*3.3. Generator*

The generator can be simplified to Equation (4) from the side of the control system. The torque command and generator speed are input from the control algorithm and the drivetrain model, respectively, to calculate the generator torque and electrical power. The electrical power is calculated using Equation (5).

$$\frac{T\_{\mathcal{S}}(s)}{T\_{\mathcal{S}}^c(s)} = \frac{1}{1 + \tau\_{\mathcal{S}}s} \tag{4}$$

$$P = T\_{\mathcal{K}} \Omega\_{\mathcal{K}} \eta\_{\mathcal{K}} \tag{5}$$

*3.4. Tower*

The tower model is expressed as the equation of motion given in Equation (6). In this model, the velocity at which the nacelle sways fore and aft because of the wind speed is added to the input wind speed and entered into the aerodynamics model.

$$k\mathbf{u}\mathbf{\dot{r}}\dot{\mathbf{x}}\_{faft} + c\mathbf{r}\dot{\mathbf{x}}\_{faft} + k\mathbf{r}\mathbf{x}\_{faft} = F\mathbf{r} \tag{6}$$

#### *3.5. Pitch Actuator*

The dynamic characteristics of the pitch actuator model are expressed by Equation (7). The pitch actuator operates within the limits of Equations (8) and (9) according to the design specification.

$$\frac{\beta(s)}{\beta^c(s)} = \frac{1}{1 + \tau\_P s} \tag{7}$$

$$
\delta - \mathsf{S}^{\diamond} \le \beta \le \theta \mathbf{0}^{\diamond} \tag{8}
$$

$$-10\ \text{(}^{\circ}\text{/s)} \le \dot{\beta} \le 10\ \text{(}^{\circ}\text{/s)}\tag{9}$$

## **4. Control Algorithms**

This section introduces PI control, LQR control, and LQR-PI hybrid control algorithms. The PI control algorithm is a control technique applied to the target wind turbines and in this study, it is presented as a reference control algorithm to compare the performance of LQR control and LQR-PI control.

#### *4.1. PI Control Algorithm*

The PI control algorithm adopted and used by the target wind turbine receives feedback on the measured pitch angle, electrical power, and generator speed, and sends pitch angle and torque commands to the pitch actuator and generator [7,10,15]. In practice, mechanical load-reduction control techniques such as tower dampers and peak shaving are usually applied, but in this study, only power control was considered.

Figure 4 shows a block diagram of the pitch PI control algorithm containing the gain schedule. The configuration consists of a pitch PI control, torque schedule, and mode switch. The torque schedule is a lookup table with an input of generator speed and an output of generator torque. It was constructed to perform an open-loop maximum power point tracking (MPPT) control to achieve the optimal tip speed ratio with the measured generator speed in region 2. The optimal values of the generator torque with respect to the input generator speeds were calculated on the basis of the aerodynamic analysis of the rotor using DNVGL-Bladed. The gain selection of pitch PI control and operation of the mode switch are explained in detail below.

**Figure 4.** Block diagram of the proportional integral (PI) control algorithm.

Figure 5 shows the block diagram of the mode switch. The mode switch determines the control mode using an internal logic with the measurement values of generator speed, electrical power, and pitch angle. A set-reset (SR) flip-flop is a logic that remembers one bit and remains in the current state until a change in the state signal (clock) is generated. If the measured generator speed or power exceeds the rated values, the mode switch outputs a signal of 1 (switched on). Also, the mode switch outputs a signal of 1 (switched on) if the measured pitch angle is greater than the fine pitch angle.

**Figure 5.** Block diagram of mode switch.

When the mode switch is on, pitch PI control is performed, and the generator torque command is fixed to be the rated value. When the mode switch is <sup>o</sup>ff, to perform open-loop MPPT control, the pitch angle command is fixed to be the fine pitch angle, and the torque control is performed through the torque schedule.

Figure 6 shows the frequency response of the pitch control loop gain. Figure 6a shows the frequency response of the open-loop transfer function given in Equation (10). The frequency response was drawn over all the wind speeds from the rated wind speed to the cutout wind speed by 0.5 m/s intervals. As shown in the figure, the frequency response varies depending on wind speed. Therefore, to have uniform pitch sensitivity, gain scheduling should be applied to maintain a constant value of the cross frequency of the pitch control loop. In this study, the cross frequency was set to 1 rad/s, taking into

account the fact that most energy components in wind speed exist at frequencies lower than 1 rad/s based on the power spectrum of wind speed [4].

**Figure 6.** Frequency response of the pitch loop. (**a**) Open-loop transfer function for the pitch input; (**b**) pitch control loop gain transfer function with gain scheduling.

Figure 6b shows the frequency response of the pitch control loop gain transfer function of Equation (11). The pitch control loop consists of gain scheduling, PI control, pitch actuator dynamics (Equation (7)) and open-loop transfer function (Equation (10)). As shown in the figure, all the frequency responses (magnitude plot) at different wind speeds had a cross frequency of 1 rad/s, and the phase margin (phase plot) of at least 30 degrees was achieved for system stability.

$$G(s) = \frac{\delta \Omega\_{\%}(s)}{\delta \beta(s)}\tag{10}$$

$$L(\mathbf{s}) = k\_{\mathbf{G}}(\boldsymbol{\beta}) \Big( \mathbf{k}\_{\mathcal{P}} + \frac{k\_{\mathbf{i}}}{\mathbf{s}} \Big) \Big( \frac{1}{\tau\_{\mathcal{P}} \mathbf{s} + 1} \Big) \Big( \frac{\delta \Omega\_{\mathcal{S}}(\mathbf{s})}{\delta \beta(\mathbf{s})} \Big) \tag{11}$$

#### *4.2. LQR Control Algorithm*

Figure 7 shows the control structure of the LQR control algorithm. The pitch control was replaced by LQR control. The LQR control received the generator speed, torque, and pitch angle as inputs. In addition, the wind speed estimator used the current generator speed, torque, and pitch angle to deliver the estimated wind speed to the designed LQR control. The design of the wind speed estimator is presented in detail in Section 5.

**Figure 7.** Block diagram of the linear quadratic regulator (LQR) control algorithm.

Linearization models are required to select the optimal gain for LQR control. In this study, a linearization model was acquired through DNVGL-Bladed. In Equation (12), the state matrix A and

input matrix B are stabilizable. The state vectors and control vectors are presented in Equation (13). In state vector *x*, the pitch angle actually reflects only one result because the target wind turbine performs collective pitch control (CPC). In order to stabilize the system Equation (12), the optimum gain *K* in Equation (17) that minimizes the quadratic cost function (Equation (15)) through the state feedback method (Equation (14)) must be selected. In Equation (17), S is the symmetric positive semidefinite solution of the Riccati Equation (16). Since the size of the matrix was large, *K* was obtained using the LQR function of MATLAB (R2014a, The MathWorks, Inc, Natick, MA, USA).

Also, the weight matrices *Q* and *R* that met the conditions in Equation (18) were chosen randomly and simulated by solving the Riccati Equation. Then *Q* and/or *R* were re-selected if transient response specifications and/or size constraints were not met. This means that the weight matrix was selected as a tuning process until a satisfactory performance was achieved. The *Q* matrix was weighted more heavily for the state variables, so that the objective was achieved in a short time, and the *R* matrix was chosen through the simulation response.

$$
\dot{\mathbf{x}} = A\mathbf{x} + Bu\tag{12}
$$

$$\mathbf{x} = \begin{bmatrix} \mathbf{x}\_{faf} \ \dot{\mathbf{x}}\_{faf} \ \mathbf{x}\_{\text{side}} \ \dot{\mathbf{x}}\_{\text{side}} \ \boldsymbol{\theta}\_{\text{g}} \ \boldsymbol{\Omega}\_{\text{g}} \ \mathbf{T}\_{\text{g}} \ \boldsymbol{\beta}\_{1} \ \dot{\boldsymbol{\beta}}\_{1} \ \boldsymbol{\beta}\_{2} \ \dot{\boldsymbol{\beta}}\_{2} \ \dot{\boldsymbol{\beta}}\_{3} \ \dot{\boldsymbol{\beta}}\_{3} \end{bmatrix}^{T}, \boldsymbol{\mu} = \begin{bmatrix} \boldsymbol{\beta}^{cmd} \end{bmatrix} \tag{13}$$

$$
\mu = -\mathbb{K}\mathfrak{x} \tag{14}
$$

$$J = \int\_0^\infty \left(\mathbf{x}^T Q \mathbf{x} + \mathbf{u}^T R \mathbf{u}\right) dt\tag{15}$$

$$A^T S + SA - SBR^{-1}B^T S + Q = 0\tag{16}$$

$$K = R^{-1}B^T S \tag{17}$$

$$Q = Q^T \ge 0,\\
R = R^T > 0\tag{18}$$

#### *4.3. LQR-PI Control Algorithm*

Figure 8 shows a block diagram of the LQR-PI control algorithm proposed in this paper. The LQR-PI control combines the control output of the PI pitch control and that of the LQR control to transmit the combined pitch command to the pitch actuator. As shown in the figure, the LQR control was applied in this structure as a feedforward term for pitch control.

**Figure 8.** Block diagram of the LQR-PI hybrid control algorithm.

Although it is simply a combination of LQR control and PI control, both controls can complement each other in the pitch control domain to improve the operating stability of the wind turbine. If the LQR control delivers the optimum pitch angle command to the pitch actuator and performs a sufficiently stable control, the contribution of the PI control will be small, because the pitch PI control will intervene in control when the measured generator speed exceeds the rated value. However, because of noise or unexpected circumstances, the LQR control may send incorrect pitch commands to the pitch actuator, thus not achieving its original purpose of maintaining the rated generator speed. In this case, PI control takes the lead in pitch control.

LQR-PI control was configured to perform a PI control also when LQR control was removed (or disconnected by a switch) from the control structure. The advantage of this feature is that when the LQR-PI control is applied to actual wind turbines, only a feed-forward loop of LQR control can be added to the pitch loop of the existing PI controller without modifying existing control algorithms.

#### **5. Wind Speed Estimator**

The nacelle wind speed is not suitable for feeding a control-scheduling logic because it is disturbed by the rotation of the rotor, which introduces a periodic decrease with multiples of the rotor frequency, as well as higher frequency disturbances due to wake turbulence [33]. Therefore, a wind speed estimator was designed for LQR control and used to calculate the rotor average wind speed.

Figure 9 shows a block diagram of the wind speed estimator. It consists of an aerodynamic torque estimator, a 3D look-up table for wind speed, and a low-pass filter. The aerodynamic torque estimator is just a Simulink representation of Equation (3), which is the two-mass drivetrain model. The aerodynamic torque was firstly estimated from the aerodynamic torque estimator using the rotor speed and the generator torque and then was supplied to the 3D lookup table as an input. Two more inputs of pitch angle and rotor speed were provided to the 3D lookup table to ge<sup>t</sup> the wind speed as an output. The wind speed from the 3D lookup table was finally low-pass filtered to remove high-frequency components [21].

**Figure 9.** Block diagram of the wind speed estimator. Low pass filter (LPF); Revolutions per minute (RPM).

The 3D lookup table in Figure 9 was created using the fminsearch function of MATLAB. The fminsearch is a function minimization algorithm based on the Nelder-Mead simplex method [34–36]. It can be applied to nonlinear functions whose derivatives are not known and is one of the most widely used function minimization algorithms for a direct search method. This method uses a simple value known as a polytope with *n* + 1 vertices (or *n* + 1 test point) in the *n* variable of the objective function. To find a value that can minimize objective function values, compare the function values at the *n* + 1 test point, avoid the test points that provide the worst function values, and repeat the reflection, contraction, and extension of the variables [36].

That is, the fminsearch finds the wind speed to minimize the error between the power calculated and the rated power, as shown in Figure 10. To calculate the electrical power, the aerodynamic torque was firstly calculated using Equation (1) with inputs of wind speed, rotor speed, and pitch angle. The wind speed, in this case, was a trial value from the fminsearch algorithm, and the others were the given inputs. The aerodynamic torque was finally multiplied by the rotor spee, and the generator loss and converted into electrical power. The error between the calculated power and the rated power was used as a cost function in the fminsearch algorithm to be minimized. At the function minimum, the wind speed could be obtained. The inputs were varied to construct a 3-D lookup table

whose output was the wind speed. The inputs in the lookup table were rotor speed, pitch angle, and aerodynamic torque.

**Figure 10.** Block diagram of a function to produce a 3D lookup table.

Figure 11 shows the results of the 3D lookup table constructed through this process. The 3D lookup table was constructed for the operation range as a function of rotor speed, pitch angle, and estimated aerodynamic torque. For example, if the estimated aerodynamic torque for a given rotor speed of 40 RPM, as in Figure 11 was 25 kNm and the pitch angle was 14◦, the wind speed would be 15 m/s.

**Figure 11.** 3D lookup table for the estimated rotor averaged wind speed.

To validate the wind speed estimator using simulation, the wind speed estimated from the wind speed estimator was compared with the rotor averaged wind speed from DNVGL-Bladed (Figure 12). A dynamic simulation at a mean wind speed of 14 m/s was performed with the target wind turbine, and the rotor averaged wind speed was obtained from DNVGL-Bladed. The generator torque, rotor speed, pitch angle from DNVGL-Bladed at a time interval of 10 ms were used as inputs to the wind speed estimator, and the turbulent wind speeds were obtained as outputs. Although a delay of less than 1 second was found, the wind speed estimator considered appeared valid because the mean and the standard deviation of the two wind speed data were almost identical.

**Figure 12.** Comparison of the wind speed estimator results with those obtained with DNVGL-Bladed.
