Active blade pitch control and stabilization of a wind turbine driven PMSG for power output regulation

Abstract

A common strategy in controlling a permanent magnet synchronous generator (PMSG) driven by a wind turbine is the maximization of output power of the wind turbine itself. A control strategy must be adopted, is to deliver a desired reduced amount of power whenever it is required. In order to realize the direct control of wind turbine output power across a wide range of wind speeds, a linearized parameter varying dynamic model of the nonlinear wind turbine system including wind disturbances is developed and used in this paper. The stability of the wind turbine system is analyzed and a blade pitch controller is designed, based on the linearized, parameter-varying, model-predictive control and is validated. Thus, the wind turbine is regulated in a way that the generator delivers the demanded power output to the load. Moreover, the blade pitch control system also performs the key function of augmenting the stability of the wind turbine, for the right choice of the gains.

Introduction

The conventional method of power generation from a wind turbine has been based on the use of a doubly fed induction generator. However, there has been a growing interest in the development of small scale wind turbine power generating units which typically drive a permanent magnet synchronous generator (PMSG). There have been many recent studies related to the dynamic model and control of wind turbines driving a PMSG. Elbeji et al. (2014) have investigated the dynamic model and control of a PMSG driven by a wind turbine. They have investigated several suitable control schemes. Melício et al. (2011) have considered the application of both classical and fractional order control theory to wind turbines driving a PMSG. Hussein et al. (2013) have considered the application simple control schemes to control both the generator side and grid side of the PMSG system. There have been several applications related to the maximum power point tracking (MPPT) control. Rolan et al. (2009) reported the implementation of a MPPT system by simple speed adjustment of the wind turbine. Aliprantis et al. (2000) discussed the modeling and active control of a stall-regulated variable-speed wind turbine driving a PMSG. Camara et al. (2015) developed a MPPT system for PMSG speed control as well as active / reactive power control, management of the DC-bus voltage and battery’s power control. Tafticht et al. (2006) have discussed the estimation of power quality, control strategies for MPPT and the connection of the wind turbine to a variety of storage or grid systems. Hamatwi et al. (2017) have implemented rotor speed control for MPPT using phase lag compensation. Other important studies related to wind turbines driving a PMSG include the control of the inertial frequency response using a full-rated Voltage Source Converters (VSC), by Cheah-Mane et al. (2014) and the development of a systematic approach to model reduction by Hackl et al.(2018). The important issue of the coupled stability and control of the permanent magnet synchronous generator (PMSG) when driven by a wind turbine was considered by Hamied and Amary (2016). Thus it appears that most of the focus of the current research has been on the control of the electrical side of PMSG driven by a wind turbine and that there are, as yet, a few important unexplored issues. Finally, Zhou and Liu (2018) have presented a active blade pitch control approach for a wind turbine based on a nonlinear PI/PD method.

In this paper, a linearized parameter varying dynamic modeling of the nonlinear wind turbine system including wind disturbances is developed. The stability of the wind turbine system is first considered and a blade pitch control system is designed, based on the linearized, parameter-varying model. The basis of the controller is model predictive control which is then validated. Thus the wind turbine is regulated in a way that the generator delivers the desired power to the load. Operation at the maximum power point is also shown to be feasible. Moreover the blade pitch control system also performs the key function of augmenting the stability of the wind turbine system.

Dynamic modeling of the PMSG and the wind turbine

The dynamics of the PMSG is modeled in the well-known d-q co-ordinates (see e.g. Vepa (2013)). The currents in the d-q axes satisfy the differential equations given by,

$$\left[\begin{array}{cc}{L}_d& 0\\ 0& L_q\end{array}\right]\frac{d}{dt}\left[\begin{array}{c}{i}_d\\ i_q\end{array}\right]=\left[\begin{array}{cc}-R_a& +L_{q}p{\omega }_m\\ -L_{d}p{\omega }_m& -R_a\end{array}\right]\left[\begin{array}{c}{i}_d\\ i_q\end{array}\right]+\left[\begin{array}{c}0\\ -p{\omega }_m{\phi }_m\end{array}\right]+\left[\begin{array}{c}{u}_d\\ u_q\end{array}\right].$$

(1)

The electro-magnetic torque generated is given by

$$T_{em}=\frac{3}{2}p\left(\left(L_d-L_q\right)i_d{i}_q+{\phi }_m{i}_q\right)$$

(2)

Furthermore the electrical frequency is related to the mechanical speed by the relation ${\omega }_e=p{\omega }_m$ , where $p$ is the number of generator pole pairs. Thus given the driving torque generated by the wind turbine, dynamics of the mechanical shaft speed of the turbine is governed by,

$$J\frac{d{\omega }_m}{\mathrm{dt}}+B{\omega }_m=T_{\mathrm{dr}}-T_{\mathrm{em}}.$$

(3)

The power output of the wind turbine can be modeled in terms of the power coefficient $C_p(\lambda ,\beta )$ which is function of both the tip speed ratio and the blade collective pitch angle, $\beta$ . The wind power output is given by,

$$P_w=T_{\mathrm{dr}}{\omega }_m=(1/2)\rho {\mathrm{AV}}_w^3{C}_p(\lambda ,\beta ),$$

(4)

where, $\lambda ={\omega }_{m}R/V_w$ is the tip speed ratio, $C_p(\lambda ,\beta )$ is the output power coefficient and the torque driving to PMSG is given by,

$$T_{\mathrm{dr}}=T_{\mathrm{wt}}=(P_w/{\omega }_m)=(1/2{\omega }_m)\rho {\mathrm{AV}}_w^3{C}_p(\lambda ,\beta ).$$

(5)

In (5), the turbine hub atmospheric air density is $\rho$ , which is assumed to be at a height $H$ , $R$ is the outer radius of the rotor. An approximate expression for the power coefficient is given by,

$$C_p(\lambda ,\beta )=0.5176(\frac{116}{{\lambda }_i}-0.4\beta -5)\exp (-\frac{21}{{\lambda }_i})+0.006795{\lambda }_i,$$

(6)

where ${\lambda }_i$ is given by,

$${\lambda }_i=\frac{1}{\frac{1}{\lambda +0.08\beta }-\frac{0.035}{1+{\beta }^3}}.$$

(7)

The maximum pitch angle ${\beta }_{max}$ is assumed to be less than $20°$ .

A generic expression for the power coefficient $C_p(\lambda ,\beta )$ is,

$$C_p(\lambda ,\beta )=C_1(\frac{C_2}{{\lambda }_i}-C_3\beta -C_4{\beta }^{C_5}-C_6)e^{-C_7/{\lambda }_i}+C_8\lambda ,$$

$$\frac{1}{{\lambda }_i}=\frac{1}{\lambda +C_9\beta +C_{12}}-\frac{C_{10}}{1+C_{11}{\beta }^3}.$$

(8)

Table 1 lists the approximations to constant parameters $C_i$ , $i=1,2,3\cdots ,12$ , in the power coefficient $C_p(\lambda ,\beta )$ (are obtained from Wijewardana et al. (2016) where they are also compared with the coefficients given by Voltolini et al. (2012) and by Heier (1998). A more complete set of coefficients is also available from, González-Hernández and Salas-Cabrera (2019). The calculations performed in this paper for the selected set of coefficients were performed for all of the available coefficient sets. Plots of a selected representation of the power coefficient $C_p(\lambda ,\beta )$ and the torque coefficient $C_p(\lambda ,\beta )/\lambda$ are respectively compared in Figure 1(a) and (b), for increasing values of the blade pitch angle.

Table 1.

Approximations to the parameters $C_i$ , in the expression for $C_p(\lambda ,\beta )$ .

Model	Our model	Voltolini et al. (2012)	Heier (1998)	Constant speed operation	Variable speed operation
$C_1$	0.5176	0.5176	0.5	0.44	0.73
$C_2$	116	116.0	116.0	125.0	151.0
$C_3$	0.4	0.4	0.4	0	0.58
$C_4$	0	0	0	0	0.002
$C_5$	0	0	0	0	2.14
$C_6$	5	5	5	6.94	13.2
$C_7$	21	21	21	16.5	18.4
$C_8$	0.006795	0.006795	0	0	0
$C_9$	0.08	0.08	0.08	0	0.02
$C_{10}$	0.035	0.035	0.035	0.002	0.003
$C_{11}$	1	1	1	1	1
$C_{12}$	0	0	0	0	0

Figure 1.

(a) Plots of the power coefficient $C_p(\lambda ,\beta )$ and (b) Plots of the torque coefficient $C_T=C_p(\lambda ,\beta )/\lambda$ .

From Figure 1, it is seen that initially, as the rotor the speed of the wind turbine is lowered, the tip speed ratio increases, resulting in an increase in the power coefficient as can be observed from the power coefficient curves. The power coefficient curve is only valid above the cut-in tip speed ratio and below the cut-out tip speed ratio. However as the tip speed ratio increases beyond the point of maximum power, the aerodynamic conditions limit, and reduce the power absorbed by the wind turbine from the flow. Consequently there is fall in the power output and the power coefficient also falls. It must also be mentioned that when the flow around a blade section separates, the lift generated by it stalls and the dynamic stall conditions severely limit the power output of the wind turbine. In these situations the power coefficient models of the type given by equation (8) are no longer valid and the effects of dynamic stall must be included in the model as mentioned by Holierhoek (2013). When dynamic stall is present, the power coefficient is further influenced by the stalled flow and the dynamic inflow that may be present, the tip (and root) corrections due to having a finite number of blades, the presence of yawed inflow, the effects of the shadow of the tower as well as three-dimension effects on the lift generated and the state of the turbulent wake. In this case the dynamics of lift generation must be included as is done with models developed for helicopters (see e.g. Truong (2017)) and aircraft propellers, using a blade-element momentum approach as indicated in Vepa (2013). The analysis of the stability and control of the wind turbine, under these conditions is well beyond the scope of this paper and will be presented elsewhere. Moreover the power output of the wind turbine is generally very low under these conditions and the dynamic model of the power generated is used only for developing control laws for active control of a stall-regulated variable-speed wind turbine.

The dynamics of the blade pitch angle actuator, which essentially determines the blade pitch angle $\beta$ , is assumed to be represented by the first order model given by,

$${\tau }_b\stackrel{·}{\beta }+\beta =u_c,$$

(9)

where ${\tau }_b$ is the blade actuator time constant, and $u_c$ is the blade angle control input which may be expressed in terms error between the blade angle feedback and the demanded blade angle as, $e=\beta -{\beta }_d$ .

Having defined the dynamic models of the generator and the wind turbine and also discussed the limitations of the models for the power coefficient, it is important to briefly discuss the distribution of wind on the Earth’s surface. At any location on the Earth’s surface, the wind velocity is distributed probabilistically. The Weibull distribution best describes the probability distribution of the wind. Since the wind power is proportional to the third power of the wind speed, as a consequence of the fact that the wind is a random variable, the actual mean wind power at a particular location can be expressed in terms of the parameters of the Weibull distribution and is given by.

$${\overline{P}}_w=T_{\mathrm{dr}}{\omega }_m=\frac{1}{2}\rho A{C}_p(\lambda ,\beta ){\overline{V}}_w^3\frac{1}{{\overline{c}}^3}\Gamma (1+\frac{3}{k}),\overline{c}=\Gamma (1+\frac{1}{k}),$$

(10)

where $\overline{c}$ and $k$ are the parameters of Weibull distribution, ${\overline{V}}_w$ is mean wind speed, $\Gamma (.)$ is the Gamma function.

Finally one may specify the desired power output as:

For ${\overline{V}}_w$ below cut-off,

$${\overline{V}}_w<V_{\mathrm{cut}\_\mathrm{off}},P_d=0;\beta ={\beta }_{max}.$$

(11)

For

$${\overline{V}}_w,V_{\mathrm{cut}\_\mathrm{off}}\le {\overline{V}}_w<V_{\mathrm{rated}}$$

$$P_d={\overline{P}}_w=\frac{1}{2}\rho A{C}_p(\lambda ,\beta ){\overline{V}}_w^3\frac{1}{{\overline{c}}^3}\Gamma (1+\frac{3}{k});\beta =0.$$

(12)

For ${\overline{V}}_w$ , $V_{\mathrm{rated}}\le {\overline{V}}_w<V_{max}$ ,

$$P_d=\frac{1}{2}\rho A{C}_p(\lambda ,0)V_{\mathrm{rated}}^3\frac{1}{{\overline{c}}^3}\Gamma (1+\frac{3}{k}),{\overline{P}}_w=\frac{1}{2}\rho A{C}_p(\lambda ,\beta ){\overline{V}}_w^3\frac{1}{{\overline{c}}^3}\Gamma (1+\frac{3}{k}).$$

(13)

For ${\overline{V}}_w$ , ${\overline{V}}_w\ge V_{max}$ , $P_d=0$ ; $\beta ={\beta }_{max}$ . Hence, given the expression for the rated $C_p(\lambda ,\beta )=C_p(\lambda ,0)(V_{\mathrm{rated}}^3/{\overline{V}}_w^3)$ and the actual operational value of the power coefficient, $C_p(\lambda ,\beta )$ , the problem is often reduced to one finding the value to be set for $\beta$ to achieve the desired value of the power coefficient, $C_p=C_p(\lambda ,\beta )$ .

De-coupling the electrical and mechanical subsystems

To simplify the design of the controller for the wind turbine it is useful to decouple the electrical and mechanical subsystems. Thus to decouple the electrical and mechanical sub-systems, re-consider the electrical sub-system dynamics given by (1),

$$\left[\begin{array}{cc}{L}_d\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& L_q\hfill \end{array}\right]\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{i}_d\hfill \\ \multicolumn{1}{c}{i_q}\end{array}\right]=\left[\begin{array}{cc}-R_a\hfill & +L_{q}p{\omega }_m\hfill \\ \multicolumn{1}{c}{-L_{d}p{\omega }_m}& -R_a\hfill \end{array}\right]\left[\begin{array}{c}{i}_d\hfill \\ \multicolumn{1}{c}{i_q}\end{array}\right]+\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{-p{\omega }_m{\phi }_m}\end{array}\right]+\left[\begin{array}{c}{u}_d\hfill \\ \multicolumn{1}{c}{u_q}\end{array}\right].$$

(14)

With electrical feedback,

$$\left[\begin{array}{c}{u}_d\hfill \\ \multicolumn{1}{c}{u_q}\end{array}\right]=\left[\begin{array}{c}+L_{q}p{\omega }_m{i}_q\hfill \\ \multicolumn{1}{c}{R_a{i}_{q0}+p{\omega }_m{\phi }_m+L_{d}p{\omega }_m{i}_d}\end{array}\right],$$

(15)

one has,

$$\left[\begin{array}{cc}{L}_d\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& L_q\hfill \end{array}\right]\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{i}_d\hfill \\ \multicolumn{1}{c}{i_q}\end{array}\right]=\left[\begin{array}{cc}-R_a\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -R_a\hfill \end{array}\right]\left[\begin{array}{c}{i}_d\hfill \\ \multicolumn{1}{c}{i_q-i_0}\end{array}\right].$$

(16)

The electro-magnetic torque generated is given by

$$T_{\mathrm{em}}=\frac{3}{2}p((L_d-L_q)i_d{i}_q+{\phi }_m{i}_q).$$

(17)

Thus the electrical and mechanical subsystems are coupled and the controllers can be independently synthesized. Moreover the stability can be independently assessed.

The state space model

The dynamics of the PMSG driven by the wind turbine may now be expressed in state space domain. Define the state vector,

$$\mathbf{x}={\left[\begin{array}{cccc}{i}_d\hfill & i_q\hfill & {\omega }_m\hfill & \beta \hfill \end{array}\right]}^T.$$

(18)

$$\left[\begin{array}{cc}{L}_d\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& L_q\hfill \end{array}\right]\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]=\left[\begin{array}{cc}-R_a\hfill & +L_{q}p{x}_3\hfill \\ \multicolumn{1}{c}{-L_{d}p{x}_3}& -R_a\hfill \end{array}\right]\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]+\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{-p{x}_3{\phi }_m}\end{array}\right]+\left[\begin{array}{c}{u}_d\hfill \\ \multicolumn{1}{c}{u_q}\end{array}\right],$$

(19)

On the other hand, based on the eCoupled model,

$$\left[\begin{array}{cc}{L}_d\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& L_q\hfill \end{array}\right]\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]=\left[\begin{array}{cc}-R_a\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -R_a\hfill \end{array}\right]\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2-x_{20}}\end{array}\right],$$

(20)

$$T_{\mathrm{em}}=\frac{3}{2}p((L_d-L_q)x_1{x}_2+{\phi }_m{x}_2).$$

(21)

The mechanical dynamics in terms of the states is,

$$J{\stackrel{·}{x}}_3+B{x}_3=T_{\mathrm{dr}}-T_{\mathrm{em}},$$

(22)

$$T_{\mathrm{dr}}=T_{\mathrm{wt}}=\frac{P_w}{x_3}=\frac{1}{2x_3}\rho {\mathrm{AV}}_w^3{C}_p(\lambda ,x_4)=\frac{1}{2}\rho {\mathrm{AV}}_w^{2}R{C}_T,$$

(23)

$$C_T\equiv C_p(\lambda ,x_4)/\lambda ,\lambda =x_{3}R/V_w,$$

(24)

$${\tau }_b{\stackrel{·}{x}}_4=u_c-x_4.$$

(25)

Linearized dynamics of the wind turbine

The complete dynamics may now be linearized and the conditions for equilibrium and stability may be established. The derivative of the mechanical torque driving the PMSG is,

$$\frac{d{C}_T}{d\lambda }=\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)=\frac{1}{\lambda }(\frac{d{C}_p(\lambda ,\beta )}{d\lambda }-\frac{C_p(\lambda ,\beta )}{\lambda }).$$

(26)

Thus, linearizing the expression for $T_{\mathrm{dr}}$ ,

$$T_{\mathrm{dr}}=\frac{1}{2}\rho {\mathrm{AV}}_w^2\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)(\lambda -{\lambda }_0)+\frac{1}{2}\rho {\mathrm{AV}}_w^2\frac{d}{d\beta }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)\beta .$$

(27)

In terms of the states of the system,

$$T_{\mathrm{dr}}=\frac{1}{2}\rho {\mathrm{AV}}_w^2\frac{d}{d\lambda }\left(\frac{C_p}{\lambda }\right)\frac{R}{V_w}(x_3-x_{30})+\frac{1}{2}\rho {\mathrm{AV}}_w^2\frac{d}{d\beta }\left(\frac{C_p}{\lambda }\right)\beta ,$$

(28)

where

$$\frac{d}{d\beta }{C}_p(\lambda ,\beta )=C_1(C_{2}g(1-C_7/{\lambda }_i)-C_3(1-C_7\beta g))-C_1(C_4{\beta }^{C_5-1}(C_5-C_7\beta g)-C_7{C}_{6}g)e^{-C_7/{\lambda }_i},$$

(29)

with,

$$g=-\frac{C_9}{{(\lambda +C_9\beta +C_{12})}^2}+\frac{3{\beta }^2{C}_{11}{C}_{10}}{{(1+C_{11}{\beta }^3)}^2}.$$

(30)

Thus, the derivative of the power coefficient $C_p(\lambda ,\beta )$ with the tip speed ratio $\lambda$ is,

$$\frac{d}{d\lambda }{C}_p(\lambda ,\beta )=C_8(1+C_7\lambda \frac{d{q}_i}{d\lambda })-C_7{C}_p(\lambda ,\beta )\frac{d{q}_i}{d\lambda }+C_1{C}_2{e}^{-C_{7}q}\frac{d{q}_i}{d\lambda },$$

(31)

where,

$$q_i\equiv \frac{1}{{\lambda }_i}=\frac{1}{\lambda +C_9\beta +C_{12}}-\frac{C_{10}}{1+C_{11}{\beta }^3},$$

(32)

$$\frac{d{q}_i}{d\lambda }=-{\left(\frac{1}{\lambda +C_9\beta +C_{12}}\right)}^2\lambda ={R{\omega }_m/V}_w.$$

(33)

The nominal dynamic pressure force on the wind turbine disk is denoted as, $f_d=(1/2)\rho V_w^{2}A$ . Hence, after linearizing, $T_{\mathrm{dr}}$ is,

$$T_{\mathrm{dr}}=\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)\frac{f_{d}R}{V_w}(x_3-x_{30})+\frac{f_d}{\lambda }\frac{d{C}_p(\lambda ,\beta )}{d\beta }\beta .$$

(34)

The linearized driving torque is expressed as,

$$T_{\mathrm{dr}}=d_{\mathrm{er}1}(x_3-x_{30})+d_{\mathrm{er}2}\beta ,$$

(35)

where the stability derivatives $d_{\mathrm{er}1}$ are $d_{\mathrm{er}2}$ are defined as,

$$d_{\mathrm{er}1}=\frac{d{C}_T}{d\lambda }\frac{f_{d}R}{V_w}=\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)\frac{f_{d}R}{V_w},$$

(36)

$$d_{\mathrm{er}2}=\frac{f_d}{\lambda }\frac{d{C}_p(\lambda ,\beta )}{d\beta }.$$

(37)

Plots of the derivatives of the power coefficient $C_p(\lambda ,\beta )$ and the derivatives of the torque coefficient $C_T=C_p(\lambda ,\beta )/\lambda$ are respectively shown in Figure 2(a) and (b), for increasing blade pitch angle.

Figure 2.

(a) Stability derivative of the power coefficient and (b) Stability derivative of the torque coefficient.

It is clear from the Figures 1 and 2 that wind turbine’s power coefficient is effectively fully controllable only beyond $\lambda$ equal to about 9. Positive $C_p(\lambda ,\beta )$ values are feasible for $9\le \lambda \le 13$ . Thus if one requires to reduce the value of the output power, or the value of $C_p(\lambda ,\beta )$ , one must necessarily increase the speed of the wind turbine. Linearizing $T_{\mathrm{em}}$ ,

$$T_{\mathrm{em}}=\frac{3}{2}p((L_d-L_q)x_{10}+{\phi }_m)x_2+\frac{3}{2}p\left((L_d-L_q)x_{20}\right)x_1.$$

(38)

Equation (38) is expressed as,

$$T_{\mathrm{em}}=d_{\mathrm{er}4}{x}_2+d_{\mathrm{er}3}{x}_1,$$

(39)

where the stability derivatives $d_{\mathrm{er}3}$ and $d_{\mathrm{er}4}$ are defined as,

$$d_{\mathrm{er}3}=(3/2)p\left((L_d-L_q)x_{20}\right)\approx 0,$$

(40)

$$d_{\mathrm{er}4}=\frac{3}{2}p((L_d-L_q)x_{10}+{\phi }_m)=\frac{3}{2}p((L_d-L_q)i_{d0}+{\phi }_m).$$

(41)

Thus the complete linear model is given by,

$$\left[\begin{array}{cc}{L}_d\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& L_q\hfill \end{array}\right]\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]=\left[\begin{array}{cc}-R_a\hfill & +L_{q}p{x}_{30}\hfill \\ \multicolumn{1}{c}{-L_{d}p{x}_{30}}& -R_a\hfill \end{array}\right]\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]+\left[\begin{array}{c}{L}_{q}p{x}_{20}\hfill \\ \multicolumn{1}{c}{-L_{d}p{x}_{10}-p{\phi }_m}\end{array}\right]x_3+\left[\begin{array}{c}{u}_d\hfill \\ \multicolumn{1}{c}{u_q}\end{array}\right],$$

(42)

or in decoupled form, with the electrical feedback,

$$\left[\begin{array}{cc}{L}_d\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& L_q\hfill \end{array}\right]\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]=\left[\begin{array}{cc}-R_a\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -R_a\hfill \end{array}\right]\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2-x_{20}}\end{array}\right],$$

(43)

where,

$$\mathbf{x}={\left[\begin{array}{cccc}{i}_d\hfill & i_q\hfill & {\omega }_m\hfill & \beta \hfill \end{array}\right]}^T,$$

(44)

The electromagnetic torque is,

$$T_{\mathrm{em}}=\frac{3}{2}p((L_d-L_q)x_1{x}_2+{\phi }_m{x}_2),$$

(45)

$$T_{\mathrm{dr}}=\frac{P_w}{x_3}=\frac{1}{2x_3}\rho {\mathrm{AV}}_w^3{C}_p(\lambda ,x_4)=\frac{1}{2}\rho {\mathrm{AV}}_w^{2}R{C}_T,$$

(46)

$$C_T\equiv C_p(\lambda ,x_4)/\lambda ,\lambda =x_{3}R/V_w,$$

(47)

$$\frac{d{T}_{\mathrm{em}}}{d{i}_q}=\frac{3}{2}p((L_d-L_q)x_1+{\phi }_m),T_{\mathrm{em}}=\frac{d{T}_{\mathrm{em}}}{d{i}_q}{x}_2.$$

(48)

The mechanical dynamics is,

$$J{\stackrel{·}{x}}_3+B{x}_3=T_{\mathrm{dr}}-T_{\mathrm{em}},T_{\mathrm{em}}=d_{\mathrm{er}4}{x}_2+d_{\mathrm{er}3}{x}_1.$$

(49)

Thus, eliminating $T_{\mathrm{dr}}$ , the linear mechanical dynamics is,

$$J{\stackrel{·}{x}}_3=-d_{\mathrm{er}3}{x}_1-d_{\mathrm{er}4}{x}_2-B{x}_3+d_{\mathrm{er}1}{x}_3+d_{\mathrm{er}2}\beta -d_{\mathrm{er}1}{x}_{30}.$$

(50)

and

$${\tau }_b{\stackrel{·}{x}}_4=u_c-x_4.$$

(51)

Equilibrium and stability of the wind turbine

The steady state conditions for equilibrium point operation are used to define the steady state rotor mechanical speed ${\omega }_{\mathrm{ms}}$ and the steady state $q$ (quadrature) axis current, $i_{\mathrm{qs}}$ .

Thus, the steady state rotor mechanical speed ${\omega }_{\mathrm{ms}}$ is,

$$B{\omega }_{\mathrm{ms}}=T_{\mathrm{dr}}-T_{\mathrm{em}}.$$

(52)

Hence, it follows that,

$$T_{\mathrm{em}}=T_{\mathrm{dr}}-B{\omega }_{\mathrm{ms}}=\frac{3}{2}p((L_d-L_q)i_d+{\phi }_m)i_q.$$

(53)

It also follows that, the steady state $q$ (quadrature) axis current, $i_{\mathrm{qs}}$ is,

$$i_{\mathrm{qs}}=2(T_{\mathrm{dr}}-B{\omega }_{\mathrm{ms}})/3p((L_d-L_q)i_d+{\phi }_m)=2(T_{\mathrm{dr}}-B{\omega }_{\mathrm{ms}})/3p{\phi }_m.$$

(54)

To consider the stability, since the electrical system is decoupled and stable it is only essential to consider the mechanical subsystem. Thus,

$${\stackrel{·}{x}}_3=-\frac{(d_{\mathrm{er}3}{x}_1+d_{\mathrm{er}4}{x}_2)}{J}+\frac{(d_{\mathrm{er}1}-B)}{J}{x}_3+\frac{d_{\mathrm{er}2}}{J}{x}_4-\frac{d_{\mathrm{er}1}}{J}{x}_{30}.$$

(55)

$${\stackrel{·}{x}}_4=-(x_4/{\tau }_b)+u_c/{\tau }_b.$$

(56)

Hence the open loop stability is guaranteed when, $d_{\mathrm{er}1}-B<0$ , or in terms of the stability derivative $d_{\mathrm{er}1}$ as,

$$d_{\mathrm{er}1}=\frac{d{C}_T}{d\lambda }\frac{f_{d}R}{V_w}=\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)\frac{f_{d}R}{V_w}<B,$$

(57)

$$\frac{d{C}_T}{d\lambda }=\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)<\frac{B{V}_w}{f_{d}R}.$$

(58)

Since the expression on the right hand side of the preceding equation is generally small as $B\approx 0$ , the stability condition is approximately reduced to,

$$\frac{d{C}_T}{d\lambda }=\frac{d}{d\lambda }\left(\frac{C_p(\lambda ,\beta )}{\lambda }\right)<\approx 0.$$

(59)

Referring to Figure 3, the actual point of operation is determined by the equilibrium point which lies on the $C_T-\lambda$ performance curve to the right of the peak of the curve. For an equilibrium point on the left any disturbance causing an increase in the mechanical speed, will result in an increased torque, which in turn will drive the wind turbine to run faster. Hence all equilibrium points to the left of the peak in the $C_T-\lambda$ curve are unstable while those to the right of the peak are stable. Thus operation at the maximum power point is also seen to be feasible as it is to the right of the peak of the $C_T-\lambda$ characteristic.

Figure 3.

Stability diagram.

To analyze the closed loop stability, assume a full state control law of the form,

$$u_c=-K_1{x}_3-K_2{x}_4.$$

(60)

The last equation is,

$${\stackrel{·}{x}}_4=-((K_2+1)x_4/{\tau }_b)-K_1{x}_3/{\tau }_b.$$

(61)

The characteristic polynomial is,

$$\left|\begin{array}{cc}{d}_{\mathrm{er}1}-B-s\hfill & d_{\mathrm{er}2}\hfill \\ \multicolumn{1}{c}{-K_1}& -K_2-1-s\hfill \end{array}\right|=0.$$

(62)

It reduces to,

$$(s+K_2+1)(s+B-d_{\mathrm{er}1})+K_1{d}_{\mathrm{er}2}=0.$$

(63)

The characteristic polynomial expands to,

$$s^2+s(K_2+1+B-d_{\mathrm{er}1})+(K_2+1)B+K_1{d}_{\mathrm{er}2}-(K_2+1)d_{\mathrm{er}1}=0.$$

(64)

Thus, in principle, $K_1$ and $K_2$ may be chosen to guarantee stability. Ensuring that the coefficient of the second order characteristic polynomial 64, are positive,

$$K_2+1+B>d_{\mathrm{er}1}$$

(65)

$$(K_2+1)B+K_1{d}_{\mathrm{er}2}>(K_2+1)d_{\mathrm{er}1}.$$

(66)

The second condition (66) may be expressed as,

$$(\frac{K_1}{K_2+1}{d}_{\mathrm{er}2}+B)>d_{\mathrm{er}1}$$

(67)

provided,

$$K_2>0.$$

(68)

In particular when $B\approx 0$ , the second condition (66) is expressed entirely in terms of the stability derivatives as,

$$\frac{d_{\mathrm{er}1}}{d_{\mathrm{er}2}}<\frac{K_1}{K_2+1}.$$

(69)

Thus the closed loop stability is guaranteed by a proper choice of $K_1$ and $K_2$ , which is feasible as long as, the stability derivative $d_{\mathrm{er}2}\ne 0$ . Clearly when $d_{\mathrm{er}2}\approx 0$ , the controller synthesis must be done while also ensuring that at the equilibrium point the stability derivative $\left|d_{\mathrm{er}2}\right|>>0$ . This may require operating with the tip speed ratio $\lambda$ being sufficiently large, which is beyond the maximum power point. This indicates that the benefits of blade pitch control can be marginal, unless the operating set point is chosen appropriately.

Control law synthesis: Nonlinear model predictive control

In this section we shall briefly consider the synthesis of a control law for the nonlinear plant by locally linearizing the plant dynamics and then applying the methodology of model predictive control at each time step. The concept of MPC is explained by Rawlings (2000). The methodology considered is based on the implementation due to Vepa (2018), but altered to suit the current application. To illustrate the process of synthesis of a linear control law at each time step, a typical discrete time system is defined as,

$$\mathbf{x}(k+1)=\mathbf{A}\left(k\right)\mathbf{x}\left(k\right)+\mathbf{B}\left(k\right)\mathbf{u}\left(k\right),\mathbf{y}\left(k\right)=\mathbf{C}\left(k\right)\mathbf{x}\left(k\right).$$

(70)

An control input sequence is defined as $\mathbf{u}\left(j\right)$ , $j=0,1,2\dots N-1$ , let the sequence of control inputs be expressed as a single vector defined by,

$$\mathbf{U}={\left[\begin{array}{cccc}{\mathbf{u}}^T\left(0\right)\hfill & {\mathbf{u}}^T\left(1\right)\hfill & \cdots \hfill & {\mathbf{u}}^T(N-1)\hfill \end{array}\right]}^T.$$

(71)

The objective is to minimize a performance index which is function of the output sequence $\mathbf{y}\left(k\right)$ , the control input sequence $\mathbf{u}\left(j\right)$ , the terminal state over the horizon, $\mathbf{y}\left(N\right)$ and a terminal weighting matrix ${\mathbf{Q}}_N$ and is assumed to be given by,

$$J(\mathbf{x}\left(0\right),\mathbf{U})=\sum _{k=0}^{N-1}\{{\mathbf{y}}^T\left(k\right)q\mathbf{Qy}\left(k\right)+{\mathbf{u}}^T\left(k\right)r\mathbf{Ru}\left(k\right)\}+{\mathbf{y}}^T\left(N\right)q{\mathbf{Q}}_N\mathbf{y}\left(N\right).$$

(72)

In (72) $q$ and $r$ are scalar scaling parameters. They are primarily used to re-scale the relative contributions of the states and the control inputs to the cost function. One may also define the sequence of state vectors expressed as a single vector as,

$$\mathbf{X}=\left[\begin{array}{ccccc}{\mathbf{x}}^T\left(1\right)\hfill & {\mathbf{x}}^T\left(2\right)\hfill & \cdots \hfill & {\mathbf{x}}^T(N-1)\hfill & {\mathbf{x}}^T\left(N\right)\hfill \end{array}\right].$$

(73)

It may be noted that the state vector $\mathbf{x}\left(k\right)$ represents the state at the following instant of time. The cost function $J$ may be written as,

$$J(\mathbf{x}\left(0\right),\mathbf{U})=q{\mathbf{x}}^T\left(0\right){\mathbf{C}}^T\left(0\right)\mathbf{QC}\left(0\right)\mathbf{x}\left(0\right)+q{\mathbf{X}}^T\overline{\mathbf{Q}}\mathbf{X}+r{\mathbf{U}}^T\overline{\mathbf{R}}\mathbf{U},$$

(74)

where the block diagonal matrix $\overline{\mathbf{Q}}$ is one with the matrices ${\mathbf{C}}^T\left(k\right)\mathbf{QC}\left(k\right)$ , $k=1,2\dots N-1$ along the diagonal except the last element which is ${\mathbf{C}}^T\left(N\right){\mathbf{Q}}_N\mathbf{C}\left(N\right)$ . The block diagonal matrix $\overline{\mathbf{R}}$ is one along the diagonal with the matrix $\overline{\mathbf{R}}$ . The next step is to construct a prediction model,

$$\mathbf{X}=\mathbf{SU}+\mathbf{Tx}\left(0\right)$$

(75)

where,

$$\mathbf{S}=\left[\begin{array}{cccc}\mathbf{B}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\mathbf{A}\left(1\right)\mathbf{B}}& \mathbf{B}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\cdots }& \cdots \hfill & \cdots \hfill & \cdots \hfill \\ \multicolumn{1}{c}{\left(\underset{k=N-1}{\overset{1}{\Pi }}\mathbf{A}\left(k\right)\right)\mathbf{B}}& \left(\underset{k=N-2}{\overset{1}{\Pi }}\mathbf{A}\left(k\right)\right)\mathbf{B}\hfill & \cdots \hfill & \mathbf{B}\hfill \end{array}\right],\mathbf{T}=\left[\begin{array}{c}\mathbf{A}\left(0\right)\hfill \\ \multicolumn{1}{c}{\mathbf{A}\left(1\right)\mathbf{A}\left(0\right)}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{\underset{k=N-1}{\overset{0}{\Pi }}\mathbf{A}\left(k\right)}\end{array}\right].$$

(76)

It has been assumed for simplicity that $\mathbf{B}$ is constant but not $\mathbf{A}$ . The cost function could be expressed compactly in the form,

$$J(\mathbf{x}\left(0\right),\mathbf{U})=\frac{1}{2}\mathbf{x}{\left(0\right)}^T\mathbf{Gx}\left(0\right)+\frac{1}{2}{\mathbf{U}}^T\mathbf{HU}+{\mathbf{x}}^T\left(0\right)\mathbf{FU},$$

(77)

with,

$$\mathbf{G}=2q({\mathbf{T}}^T\overline{\mathbf{Q}}\mathbf{T}+{\mathbf{C}}^T\left(0\right)\mathbf{QC}\left(0\right)),\mathbf{H}=2r\overline{\mathbf{R}}+2q{\mathbf{S}}^T\overline{\mathbf{Q}}\mathbf{S}\mathrm{and}\mathbf{F}=q{\mathbf{T}}^T\overline{\mathbf{Q}}\mathbf{S}.$$

(78)

To obtain the optimal control sequence one may set the gradient of the cost function with respect to the sequence of control inputs, $\mathbf{U}$ to zero. By this process the cost function $J$ is minimized resulting in,

$$\mathrm{dJ}(\mathbf{x}\left(0\right),\mathbf{U})/d\mathbf{U}={\mathbf{U}}^T\mathbf{H}+{\mathbf{x}}^T\left(0\right)\mathbf{F}=0\Rightarrow \mathbf{HU}+{\mathbf{F}}^T\mathbf{x}\left(0\right)=0\Rightarrow \mathbf{U}=-{\mathbf{H}}^{-1}{\mathbf{F}}^T\mathbf{x}\left(0\right).$$

(79)

Considering only the receding horizon, the optimal control input sequence is given by,

$$\mathbf{u}\left(k\right)=-\left[\begin{array}{cccc}1\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right]{\mathbf{H}}^{-1}{\mathbf{F}}^T\mathbf{x}\left(k\right).$$

(80)

It is recursively obtained over successive prediction windows of the control sequence. The product ${\mathbf{H}}^{-1}{\mathbf{F}}^T$ is expressed as.

$${\mathbf{H}}^{-1}{\mathbf{F}}^T=\frac{q^2}{2}{(\frac{r}{q}\overline{\mathbf{R}}+q{\mathbf{S}}^T\overline{\mathbf{Q}}\mathbf{S})}^{-1}{\mathbf{S}}^T\overline{\mathbf{Q}}\mathbf{T}.$$

(81)

The parameter $q$ may be set to 1, and (81) reduces to,

$${\mathbf{H}}^{-1}{\mathbf{F}}^T=\frac{1}{2}{(r\overline{\mathbf{R}}+{\mathbf{S}}^T\overline{\mathbf{Q}}\mathbf{S})}^{-1}{\mathbf{S}}^T\overline{\mathbf{Q}}\mathbf{T},$$

(82)

where $r$ is treated as a free parameter to be chosen.

It must be emphasized that the MPC based control ensures the closed loop stability of the wind turbine as it is based on the theory of linear optimal control. Thus it is important to focus on the steady state operating point.

Determining the optimum operating set point

To establish the steady state operating point, it is important to recognize the desired output power is our primary requirement. Thus based on the demanded power and the current operating speed, the desired steady state operating point is determined. The desired power coefficient is then determined from the demanded power output and the nominal or current operating speed. Once desired power coefficient is known, for the given tip speed ratio $\lambda$ , the desired maximum value of the pitch angle command is determined from the power coefficient, $C_p=C_p(\lambda ,\beta )$ model defined in (8) relating the power output coefficient to blade collective pitch angle $\beta$ and the tip speed ratio $\lambda$ . This requires inverting the relationship $C_p=C_p(\lambda ,\beta )$ given by (8) for a fixed tip speed ratio $\lambda$ to determine the blade collective pitch angle, $\beta ={\beta }_{d,max}$ . This inversion is done by minimizing the norm of the error between the desired power output coefficient and the mathematical model relating it to the tip speed ratio $\lambda$ and blade collective pitch angle, $\beta$ . Thus this is done by employing the Gray Wolf optimization of Long and Xu (2016) as well as by the MATLAB function provided for unconstrained optimization, “fminunc.m” for obtaining the solution for desired blade collective pitch angle, $\beta ={\beta }_{d,max}$ . Both the Gray Wolf optimization algorithm of Long and Xu (2016) as well as by the MATLAB function provided for unconstrained optimization, “fminunc.m,” were extensively compared. It was found that the results obtained were not different and the minimization time taken by the Gray Wolf optimization algorithm of Long and Xu (2016) was always consistently slightly less than the MATLAB function provided for unconstrained optimization, “fminunc.m.”

It must be recognized that although the commanded pitch angle is, ${\beta }_{d,max}$ , the actual steady state pitch angle is $0\le {\beta }_{\mathrm{ss}}\le {\beta }_{d,max}$ . The control law only tracks the desired power and maintains stability, but does not track ${\beta }_{d,max}$ . This is because the mechanical speed and hence the tip speed ratio $\lambda$ is not forced to track a demanded value but is a free variable.

It must be said that the MPPT algorithm may also be implemented. It is exactly same algorithm implemented by Vepa (2011). Briefly the current drawn from the generator is increased slowly till the maximum power point is reached, so the wind turbine is operating just beyond the maximum of the $C_p-\lambda$ curve. The details can be found in Vepa (2011) and will not be repeated here.

Finally the case of active stall control is not considered as in this case the relation for the power coefficient given by (8) is no longer valid. This case will be considered independent using a set of complementary controller synthesis tools and reported separately.

Typical simulations and results

A typical three bladed wind turbine is simulated and controlled as application of the above theory. The system parameter are listed in Table 2. The nominal time step for integrating the equations of motion is $\Delta t=0.0001$ s. The prediction window to obtain the control input at each time instant is 10 time-steps. The actual nominal power output of the wind turbine is about 60 kW. The desired power output is nominally set at 50 kW but could be chosen as desired.

Table 2.

Typical parameter values for simulation.

Parameter	Value
$C_{p\_\mathrm{nom}}$	0.38
$R$	5 m
$R_{\mathrm{root}}$	9.3564 m
$J$	$5.\times {10}^{-5}\mathrm{kg}{m}^2$
$L_q$	$5\mathrm{mH}$
$L_d$	$5\mathrm{mH}$
${\phi }_m$	$4.4\mathrm{Nm}/{(\mathrm{rad}/s)}^2$
$P$	10
$B$	$0.2\times J$
$V_w$	10 m/s
${\tau }_b$	0.01
$P_{\mathrm{nom}}$	60 kW
${\lambda }_{\mathrm{nom}}$	8
$N$	3
$H$	60 m

Initially the wind turbine is operating in steady state with $i_q$ equal to 40 amps, $i_d$ equal to 0, the power output being 62.58 kW, and the mechanical speed being 23.7 rads/s. It is sought to reduce the power output to 0.7 of the initial power output. The change is initiated after 0.04 seconds and the results are shown in Figure 4. From Figure 4 it is clear that the power responds to the desired power output that has been set. This implies a reduction in the q axis current, $i_q$ . The current in the d axis continues to remain at zero. However it is observed that the mechanical speed has gone up implying an increase in the tip speed ratio, as well, which is up by almost 50%. In Figure 5 the change in the collective blade pitch angle is plotted.

Figure 4.

Current, mechanical speed, and power output response of the wind turbine in response to power reduction command.

Figure 5.

Response of the blade collective pitch angle during the reduction in power.

The steady state pitch angle is observed to be just over 1° while the commanded pitch angle (not shown on the Figure 5) was about 3.4°. Thus although the control system does not track the commanded pitch angle it does respond to the power demanded set point and delivers the demanded power in the steady state.

Discussion and conclusions

In this paper, de-coupled dynamic and equilibrium models of a PMSG driven by a variable speed wind turbine were presented and used in the analysis of the stability and the synthesis of a power output controller. A key feature of the paper is the systematic approach the problem was approached to guarantee both the mechanical stability and independent control of the mechanical dynamics of the wind turbine.

A representative $C_p=C_p(\lambda ,\beta )$ was chosen to present the results of this paper. However these calculations have been done for all of the eight models in the same class, quoted in González-Hernández and Salas-Cabrera (2019). It must be said that it is not matter of choosing between these models, as each of these models were postulated from data obtained from experiments done under different conditions, which we were not able to ascertain. Critical reviews of the various classes of models in the literature have been presented by Reyes et al. (2015) and by Sohoni et al. (2016). Ideally by introducing additional non-dimensional parameters to represent the conditions under which the experiments were done, such as cut-in and cut-off behavior, effect of the height where the wind turbine is sited and the atmospheric properties at this height, as well as other factors, the models quoted in González-Hernández and Salas-Cabrera (2019) could be integrated into a single model. But this was not possible in the present work.

The generator is modeled in the synchronous rotating d-q reference frame and de-coupled from the turbine by feedback. Based on de-coupled dynamic model, parameter bounds for the stability of the wind turbine are derived. The stabilization and the active blade pitch control of the wind turbine with the PMSG for power output regulation has been validated and has been successfully demonstrated. The issue of the instability of the wind turbine has been addressed by a linear feedback control law that stabilizes the system and also delivers the desired power output. The control law, which is designed by applying the MPC procedure, and which has the structure of a proportional-derivative stabilizing control law that ensures the demanded power output is delivered, has been validated. Moreover the power delivered is tracked by consideration of the equilibrium conditions, using a generic non-dimensional mathematical model for the power-speed characteristics. For this reason the controller can be applied to any wind turbine driving a PMSG, even when a turbine-specific matched mathematical model is not available.