Wind Turbine Control Using Nonlinear Economic Model Predictive Control over All Operating Regions

Abstract

With the gradual increase in the installed capacity of wind turbines, more and more attention has been paid to the economy of wind power. Economic model-predictive control (EMPC) has been developed as an effective advanced control strategy, which can improve the dynamic economy performance of the system. However, the variable-speed wind turbine (VSWT) system widely used is generally nonlinear and highly coupled nonaffine systems, containing multiple economic terms. Therefore, a nonlinear EMPC strategy considering power maximization and mechanical load minimization is proposed based on the comprehensive VSWT model, including the dynamics of the tower and the gearbox in this paper. Three groups of simulations verify the effectiveness and reliability/practicability of the proposed nonlinear EMPC strategy.

1. Introduction

As a consequence of energy shortages around the world, environmental protection requirements, and the higher cost of traditional power, renewable energy sources are receiving a great deal of attention worldwide nowadays. Wind energy becomes the dominant renewable energy source due to its plentifulness, cleanness, and economic advantages. By the ending of 2018, the cumulative installed capacity of wind is estimated to be more than 221 GW in China, with an increase of 19 GW and 25.9 GW in 2017 and 2018 respectively, accounting for 31.9% and 32.8% of global overall wind power capacity [1]. Therefore, enhancing the efficiency and economy of wind energy conversion systems has become an important issue to handle the growing need for energy.

Numerous control schemes have been applied in variable-speed wind turbine (VSWT) controller designs to enhance efficiency. Generally, PID controllers are widely used in the pitch control of wind turbine systems, which only focuses on decreasing the blades’ mechanical loads within region three, as shown in Figure 1 [2]. A feedback gain-scheduled PI controller was designed to maximize the captured power and mitigate the mechanical loads, where the fixed-phase and gain-margins are obtained through a frequency response analysis [3]. In [4], a gain-scheduled L1 adaptive optimal controller is proposed for the VSWT control system, where the nonconvex optimization problem is solved by using a genetic algorithm based on a single linearized model of the wind-turbine (WT) system at different operating points. To obtain the rated generator power and reduce the mechanical load, a nonlinear-PI controller with an extended state-observer was developed in order to control the blade pitch angle [5]. Although these PI controllers have led to the effective steady-state performance of VSWTs, the distortion of the transient performance according to variations in the wind speed is still a drawback. Additionally, PI/PID controllers offer good performance, but they still are not optimal [6].

However, owing to multiple constraints, multiple variables, the highly coupled nonaffine basis of the VSWT systems, and the stochastic shape of wind speed input, high-performance wind turbine control (WTC) has become a challenging problem [7]. Various nonlinear control schemes have been applied to the WTC problem. In [8], an artificial neural network technique is proposed for the pitch control in region three, where the VSWT dynamics are modeled using a back-propagation learning algorithm. In [9], a fuzzy logic pitch controller is proposed for WTC, considering the generator power and speed as control inputs. In [10], a multiple feedback pitch controller is proposed for the nacelle fore-aft speed and generator power regulation. This controller is designed by employing a trade-off between pitch angle fluctuations and output power fluctuations. In [11], an adaptive neural pitch controller is proposed for the WTC problem. An online, two-layer neural network model is proposed for estimating the unknown wind turbine aerodynamics. Recently in [12], a T-S fuzzy modeling method with nonlinear consequents is proposed for the VSWT control problem. However, the wind turbine actuator constraints are not considered in the controller designs presented in these approaches, which may cause a wind-up phenomenon and significant performance deterioration if the actuator inputs reach saturation limits.

Among all the control strategies, model predictive control (MPC) is emerging as a powerful candidate for WTC thanks to its capacity to handle all the system constraints directly, such as physical limitations and operational constraints [13]. Classical MPC schemes have been used for WTC to ensure that the maximum power is captured in region two, and the rated power is tracked in region three. In [14,15], a linear MPC strategy is proposed based on a single linear state-space model to deal with the actuator constraints. Although these approaches have good performance concerning certain operating points and system constraints handling, no special attention has been paid to either the strong nonlinearity of the system or the mechanical load mitigation to enhance the economy of the VSWT system.

Some MPC techniques have been utilized by taking the wind turbine nonlinearity into account. In [16,17,18], a multivariable MPC is proposed for VSWT pitch control in region three, using Takagi-Sugeno (T-S) fuzzy models. In [19], a nonlinear MPC (NMPC) strategy is investigated for VSWT control for both regions two and three by tuning the penalty parameters. In [20,21], multiple model predictive controllers are proposed based on multiple linearized time-invariant state-space models, guaranteeing the actuator limits in all operation regions. However, it is obvious that the controller switching results in more difficulties. The authors take the load of WT into account as one of completing penalties based on a simplified nonlinear WT model. The shortcoming of these control schemes is the difficulty of adequately tuning the tracking weights; furthermore, economic factors are not considered in the design of the controller.

To enhance the dynamic economy, the economic MPC (EMPC) was developed for the VSWT system. Generally, the EMPC uses a measure of the system performance directly as an economic cost function [22]. However, the economic cost function may not be a positive definite with regard to any tracking trajectory in order to guarantee the stability of the operating points. In [23], the closed-loop stability of the system is guaranteed via terminal constraints depending on the strong duality assumption. In addition, in [24], a Lyapunov-based EMPC is constituted using two separate operation modes to make the closed-loop system ultimately bounded within a small region. In the context of the WTC problem, the EMPC objective function directly considers VSWT energy generation maximization against the actuator fatigue minimization, as opposed to tracking some reference states and inputs. The proposed NEMPC in this paper has a unique weight in the objective function, while the weights must be chosen carefully for every input and state in the classical NMPC strategy [19]. In addition, in [25,26], an EMPC strategy has been designed for wind turbine control systems for all regions. In [26], the wind speed is produced based on the Van der Hoven spectrum [27], in which the wind speed is considered as a slowly-changing average wind speed superposed by a rapidly-changing turbulence wind speed. However, these EMPC approaches pay no special attention to the mitigation of mechanical loads on the tower in the economic objective function. Additionally, different treatment of the random wind speed can be considered, consisting of summing the mean wind speed, which follows a two-parameter Weibull distribution, and the turbulent wind speed, which follows a zero-mean Gaussian, white noise distribution [28].

In this study, a nonlinear EMPC (NEMPC) strategy is proposed for all the operating regions of the VSWT system. A nonlinear VSWT model, which considers the tower and gearbox dynamics, is established. The NEMPC optimization problem considers all the actuator constraints (pitch angle, and torque constraints with their rate of change constraints) and the hard constraints (rotor speed, generator speed, and electrical power). The economic cost function seeks the maximum VSWT generated power against the competing penalties regardless of wind input, in order to achieve the best economic operation, as well as the fatigue load mitigation of key mechanical structures, including both drive-shaft torsion and tower fore-aft motion. This control strategy can provide potential improvements in the closed-loop performance, and satisfy the economy in contrast with classical WTC strategies. Finally, in the simulations, random wind speed is considered by summing the mean wind speed and the turbulent wind speed, which follow a two-parameter Weibull distribution and a zero-mean Gaussian white noise distribution, respectively.

The paper is organized as follows. The nonlinear modeling of the VSWT system is established in Section 2. The proposed NEMPC strategy is derived in Section 3. In Section 4, the simulation results for a 5-MW VSWT system are demonstrated, with comparisons between the classical NMPC strategy and the proposed NEMPC strategy. Also, in this section, a typical wind turbine benchmark simulator FAST is used to test the validity and practicability of the proposed NEMPC strategy. Lastly, Section 5 provides conclusions.

2. Wind Turbine Modeling

Wind-energy conversion systems consist of a windmill, a gearbox, and a generator, as demonstrated in Figure 2. In this section, nonlinear VSWT modeling is investigated. To establish this nonlinear model, the aerodynamics, tower, drive train, and the generator should be considered. In the following subsections, these subsystems are discussed separately. Finally, a nonlinear state-space model is established.

2.1. Rotor Aerodynamics

The airflow on VSWT blades causes aerodynamic torque, which represents one of the causes of VSWT nonlinearities. The aerodynamic torque $T_r$ can be stated as below [29]:

$$T_r=\frac{1}{2}\rho \pi R^2\frac{v_r^3}{w_r}{C}_p(\lambda (v_r,w_r),\beta )$$

(1)

The aerodynamic power that the VSWT extracts is expressed by [29]:

$$P^A=\frac{1}{2}\rho \pi R^2{v}_r^3{C}_p(\lambda (v_r,w_r),\beta )$$

(2)

where $\rho$ is the air density, $v_r$ represents the effective wind speed on the turbine rotor [30], $R$ is the rotor radius of VSWT, $w_r$ is the VSWT rotor angular speed, $C_p(\lambda (v_r,w_r),\beta )$ is the aerodynamic power coefficient, and $\beta$ is the blade pitch angle. The tip speed ratio $\lambda$ is defined as [30]:

$$\lambda (w_r,v_r)=\frac{R{w}_r}{v_r}$$

(3)

The following approximate expression of the aerodynamic power coefficient derived from [18] is widely used:

$$C_p(\lambda ,\beta )=c_1(\frac{c_2}{\Lambda }-c_3\beta -c_4)\times e^{-\frac{c_5}{\Lambda }}+c_6\lambda$$

(4)

where $1/\Lambda =1/(\lambda +0.08\beta )-0.035/(1+{\beta }^3)$, and the coefficients $c_i(i=1,2,\cdots ,6)$ depend on the blade shape and its aerodynamic performance [18,29]. In this paper, NREL 5-MW VSWT system is discussed; the six chosen coefficient values are [18]: $c_1$ = 0.5176, $c_2$ = 116, $c_3$ = 0.4, $c_4$ = 5, $c_5$ = 21, $c_6$ = 0.0068. This nonlinear function $C_P(\lambda ,\beta )$ is indicated in Figure 3.

2.2. Tower Dynamics

The blade vibration results in the tower deflection. Thus, the thrust force on a tower top over the nacelle, which causes undesired nodding for the tower and fatigue loads on the VSWT, must be considered [31]:

$$F_t=\frac{1}{2}\rho \pi R^2{v}_r{}^2{C}_t(\lambda (v_r,w_r),\beta )$$

(5)

where $C_t(\lambda (v_r,w_r),\beta )$ is the thrust coefficient, which is the nonlinear function in terms of blade pitch angle $\beta$ and tip speed ratio $\lambda$. By adopting a polynomial fitting algorithm based on the data-driven from a typical wind turbine benchmark simulator, the thrust coefficient can be derived as:

$$\begin{array}{ll}{C}_t(\lambda ,\beta )=& 0.08698-0.003371\beta -0.053272\lambda +0.06499{\lambda }^2-0.003096\beta {\lambda }^2-0.009575{\lambda }^3\\ & +0.0002103\beta {\lambda }^3+0.0005667{\lambda }^4-5.24\times {10}^{-6}\beta {\lambda }^4-1.199\times {10}^{-5}{\lambda }^5\end{array}$$

(6)

where $C_t(\lambda ,\beta )$ represents a nonlinear function, as indicated in Figure 4.

Thrust force towards the tower top over the nacelle causes tower fore-aft motion. The dynamics of tower fore-aft motion can be stated as a simplified second-order differential equation:

$$M_t\ddot{y}+B_t\dot{y}+K_{t}y=F_t$$

(7)

where $M_t$ represents the mass model, $B_t$ and $K_t$ are the structural damping and structural stiffness coefficients of the tower model respectively, $y$ represents the tower top displacement, and $\dot{y}$ is the tower deflection rate.

Tower fore-aft motion causes a nonnegligible effect on effective wind speed. The tower deflection rate is defined as $v_t=\dot{y}$. Thus, the effective (relative) wind speed $v_r$ on the turbine rotor can be expressed with the normal wind speed $w$ and the tower deflection as follows [30]:

$$v_r=w-v_t$$

(8)

The above formula shows that the actual wind speed acting on the blade is a relative wind speed, which is the difference between the natural wind speed $w$ and the displacement speed of the tower $v_t$. It can be imagined that when the tower moves forward and backward in the horizontal direction with the wind speed, the wind speed acting on the blade is not exactly equal to the wind speed directly given by the natural environment. This formula makes the wind speed acting on the blade more accurate.

In Figure 2, $w$ represents the nature wind speed towards the blade, which cannot be replaced by the relative (effective) wind speed $v_r$ because the relative (effective) wind speed is just a virtual physical variable which is defined for the aerodynamic torque and power, not a real physical variable. Therefore, the relative (effective) wind speed $v_r$ cannot appear before the blade, replacing wind speed $w$ in Figure 2.

2.3. Drive Train

The wind mechanical energy is transformed through a drive train to the electrical generator. The drive train is rigidly coupled in the rotor side with a flexible connection in the generator side, as demonstrated in Figure 2. The shaft torsion represents the torque of the drive train flexible shaft, which is a vital factor affecting economic performance. Oscillations in the shaft torsion cause damage to the turbine components. The shaft torsion is replaced by the absolute angular position of the generator and rotor; its derivative is [17]:

$$\dot{\theta }=w_r-\frac{w_g}{N_g}$$

(9)

where $\theta$ is the shaft torsional angle, $w_g$ is the generator angular speed, and $N_g$ is the gearbox ratio.

Assuming that the low-speed shaft is one spring and one damper, the two-mass model is used to describe the drive train from a low- to a high-speed shaft through a gearbox as follows [17,26]:

$$J_r{\dot{w}}_r=T_r-K_{\theta }\theta -B_{\theta }\dot{\theta }$$

(10)

$$J_g{\dot{w}}_g=-T_g+\frac{K_{\theta }}{N_g}\theta +\frac{B_{\theta }}{N_g}\dot{\theta }$$

(11)

where $J_r$ is the inertia of the VSWT rotor and the low-speed shaft, $J_g$ is the inertia of the generator, $K_{\theta }$ and $B_{\theta }$ represent the stiffness and damping of the drive train, $T_g$ is the generator torque.

2.4. Generator

The electrical generator power $P^E$ can be derived as [17]:

$$P^E={\eta }_g{w}_g{T}_g$$

(12)

where ${\eta }_g$ represents the generator and power efficiency. Ignoring the losses, the derivation between $P^A$ (the mechanical power) and $P^E$ stems notably from the energy tentatively stored as the rotor kinetic energy [25].

2.5. The Nonlinear State-Space Model of VSWT

Based on Equations (1)–(12), the nonlinear dynamics of the VSWT system can be transformed into the following nonlinear state-space form:

$$\begin{array}{c}\dot{x}=f_c(x,u,w)\\ y=g_c(x,u,w)\end{array}$$

(13)

where $g_c(x,u,w)=\left[\begin{array}{c}{w}_g\\ {\eta }_g{w}_g{T}_g\end{array}\right]$, $f_c(x,u,w)=\left[\begin{array}{c}-\frac{B_{\theta }}{J_g{N}_g^2}{w}_g+\frac{B_{\theta }}{J_g{N}_g}{w}_r+\frac{K_{\theta }}{J_g{N}_g}\theta -\frac{1}{J_g}{T}_g\\ \frac{B_{\theta }}{J_r{N}_g}{w}_g-\frac{B_{\theta }}{J_r}{w}_r-\frac{K_{\theta }}{J_r}\theta +\frac{\pi }{2J_r}{\rho }_a{R}_r^2\frac{{(w-v_t)}^3}{w_r}{C}_p(\lambda ,\beta )\\ -\frac{w_g}{N_g}+w_r\\ v_t\\ -\frac{K_t}{M_t}y-\frac{B_t}{M_t}{v}_t+\frac{\pi }{2M_t}{\rho }_a{R}_r^2{(w-v_t)}^2{C}_t(\lambda ,\beta )\end{array}\right]$.

Define the state variables, input variables, and output variables as: $x={\left[\begin{array}{ccccc}{w}_g& w_r& \theta & y& v_t\end{array}\right]}^T,u={\left[\begin{array}{cc}{T}_g& \beta \end{array}\right]}^T,y_o={\left[\begin{array}{cc}{w}_g& P^E\end{array}\right]}^T$. Then, the discrete wind turbine system model (13) can be obtained by forth stage Runge–Kutta (RK4) method with the sampling time $T$,

$$\begin{array}{c}x(k+1)=f(x(k),u(k),w(k))=x(k)+\frac{T}{6}(k_1+2k_2+2k_3+k_4)\\ y(k)=g(x(k),u(k),w(k))\end{array}$$

(14)

where $k_1=f_c(x(k),u(k))$, $k_2=f_c(x(k)+\frac{T}{2}{k}_1,u(k))$, $k_3=f_c(x(k)+\frac{T}{2}{k}_2,u(k))$, $k_4=f_c(x(k)+T{k}_3,u(k))$, $g(x(k),u(k),w(k))=g_c(x,u,w)$.

3. NEMPC Strategy for VSWT Control

EMPC for WTC seeks the optimal operation between the maximization of the generated power and the minimization of the cost related to care and maintenance. In a real wind farm, the oscillatory transient of tower deflection and shaft torsion can create microcracks in the materials, which can lead to component failure and increase maintenance costs. Therefore, not only generated power, pitch-angle and generator-torque, but also shaft torsion angle and tower displacement are considered in the economic cost function for the reduction of fatigue loads in this paper.

3.1. Economic Cost Function

In this paper, economic indexes are defined as follows. Firstly, the generator power must be considered in the economic index $l_{e1}$ to capture the maximum power [19]:

$$l_{e1}(x,u)=-{\lambda }_1{P}^E$$

(15)

Secondly, the shaft torsional angle and tower displacement must be considered in the economic index $l_{e2}$ to reduce the fatigue of the tower structure caused by the tower deflection [19] and the gearbox load caused by the drive-shaft torsion [26]:

$$l_{e2}(x,u)={\lambda }_2{\theta }^2+{\lambda }_3{v}_t^2$$

(16)

Finally, the pitch-angle and generator-torque must be considered in the economic indexes $l_{e3}$ and $l_{e4}$, respectively, to smoothen the control performance and reduce the fluctuations of the output electric power [19]:

$$l_{e3}(x,u)={\lambda }_4\Delta {\beta }^2+{\lambda }_5{\beta }^2$$

(17)

$$l_{e4}(x,u)={\lambda }_6\Delta T_g^2$$

(18)

where ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$, ${\lambda }_4$, ${\lambda }_5$, and ${\lambda }_6$ are weighted coefficients.

In designing the NEMPC for the VSWT system, the economic cost function can be obtained as:

$$l_e(x,u)=l_{e1}+l_{e2}+l_{e3}+l_{e4}$$

(19)

which aims to decrease the load and fatigue on the wind power system and smooth the generator torque input and the output power while capturing the maximal power.

3.2. The Operational and Physical Constraints

Due to the electrical limitations of the actuator electronics and the safety requirements of the control process, a set of physical constraints must be fulfilled during operation in VSWT. Thus, the following set of constraints must be taken into consideration along the predictive/control horizon:

$$\begin{array}{ccc}{w}_{g\min }\le w_g\le w_{g\max }& w_{r\min }\le w_r\le w_{r\max }& \theta \ge \end{array}{\theta }_{\min }$$

(20)

$$\begin{array}{cc}0\le T_g\le T_g{}_{\max }& -\Delta T_g{}_{\max }\le \Delta T_g\le \Delta T_g{}_{\max }\end{array}$$

(21)

$$\begin{array}{cc}{\beta }_{\min }\le \beta \le {\beta }_{\max }& -\Delta {\beta }_{\max }\le \Delta \beta \le \Delta {\beta }_{\max }\end{array}$$

(22)

$$0\le P^E\le P^E{}_{\max }$$

(23)

where the maximum rotor speed $w_{r\max }$, generator speed $w_{g\max }$, and torque $T_g{}_{\max }$ are set slightly higher than their rated values $w_{rrated}$, $w_{grated}$, and $T_{grated}$. The definition and values of related parameters are listed in

Table A2. The values of the bounds in the 5-MW VSWT system.

Parameter	Symbol	Value
Rated generator speed, Min. generator speed	${\omega }_{g rated}, {\omega }_{g min}$	122.9096, 70.1622
Rated rotor speed, Min. rotor speed	${\omega }_{r rated}, {\omega }_{r cut-in}$	1.2671, 0.7226
Min. shaft torsion angle	${\theta }_{min}$	0
Rated torque, Min. torque	$T_{g rated}, T_{g min}$	43,093.55, 0
Max. torque rate, Min. pitch, Max. pitch	$\Delta T_{gmax}, {\beta }_{min},{\beta }_{max}$	15,000, 0, 90
Max. pitch rate, Rated electrical power	$\Delta {\beta }_{min}, P_{rated}^E$	8, 5,000,000

in Appendix A [19,32].

3.3. The Optimization Problem for the Proposed NEMPC Strategy

Based on the above nonlinear model and cost function in Section 3.1 and Section 3.2, the optimization problem for the proposed NEMPC strategy can be summed up as follows:

$$\begin{array}{l}\min J={\displaystyle \sum _{i=1}^{N_P}{l}_e(x(k+i-1|k),u(k+i-1|k))}\\ s.t.\\ x(k+i|k)=f(x(k+i-1|k),u(k+i-1|k),w(k+i-1)),i=1,\cdots ,N_P,x(k|k)=x(k)\\ y(k+i-1|k)=g(x(k+i-1|k),u(k+i-1|k),w(k+i-1))\\ \Delta u(k+i|k)=u(k+i|k)-u(k+i-1|k)\\ equation(20)-equation(23)\end{array}$$

(24)

Define $X={\left[\begin{array}{cccc}x{(k+1|k)}^T& x{(k+2|k)}^T& \cdots & x{(k+N_P|k)}^T\end{array}\right]}^T$, $U={\left[\begin{array}{cccc}u{(k|k)}^T& u{(k+1|k)}^T& \cdots & u{(k+N_P-1|k)}^T\end{array}\right]}^T$, $\mathbf{\Delta }\mathbf{U}={\left[\begin{array}{cccc}\Delta u{(k|k)}^T& \Delta u{(k+1|k)}^T& \cdots & \Delta u{(k+N_P-1|k)}^T\end{array}\right]}^T$, $P^E={\left[\begin{array}{ccc}{P}^E(k|k)& \cdots & P^E(k+N_P-1|k)\end{array}\right]}^T$, where $N_P$ is the predictive/control horizon.

Define $\Phi ={\left[\begin{array}{cccc}X& U& \mathbf{\Delta }\mathbf{U}& P^E\end{array}\right]}^T$ as the new optimal variables set. Then, the optimization problem (24) can be rewritten as:

$$\begin{array}{l}\min J\\ s.t.\\ Aeq\Phi =bu(k-1),\\ Ce{q}_i(\Phi )={\left[\begin{array}{c}x(k+i|k)-f(x(k+i-1|k),u(k+i-1|k),w(k+i-1))\\ {\widehat{P}}^E(k+i-1|k)-g(x(k+i-1|k),u(k+i-1|k),w(k+i-1))\end{array}\right]}_{i=1,\cdots ,N_P,\widehat{x}(k|k)=x(k)}=0\\ {\Phi }_{\min }\le \Phi \le {\Phi }_{\max }\end{array}$$

(25)

where $\begin{array}{l}Aeq={\left[\begin{array}{cccc}0& I& -S& 0\end{array}\right]}_{2N_P\times 10N_P},b={\left[\begin{array}{cccc}I& I& \cdots & I\end{array}\right]}_{2N_P\times 2N_P}^T\\ S=\left[\begin{array}{cccc}I& 0& \cdots & 0\\ I& I& & 0\\ ⋮& ⋮& \ddots & \\ I& I& \cdots & I\end{array}\right]\end{array}$, $\begin{array}{l}{\Phi }_{\min }=\left[\begin{array}{cccc}{X}_{\min }& U_{\min }& \Delta U_{\min }& {\widehat{P}}_{\min }^E\end{array}\right],{\Phi }_{\max }=\left[\begin{array}{cccc}{X}_{\max }& U_{\max }& \Delta U_{\max }& {\widehat{P}}_{\max }^E\end{array}\right]\\ X_{\min }=\left[\begin{array}{c}{x}_{\min }\\ x_{\min }\\ ⋮\\ x_{\min }\end{array}\right],X_{\max }=\left[\begin{array}{c}{x}_{\max }\\ x_{\max }\\ ⋮\\ x_{\max }\end{array}\right],U_{\min }=\left[\begin{array}{c}{u}_{\min }\\ u_{\min }\\ ⋮\\ u_{\min }\end{array}\right],U_{\max }=\left[\begin{array}{c}{u}_{\max }\\ u_{\max }\\ ⋮\\ u_{\max }\end{array}\right]\\ \Delta U_{\min }=\left[\begin{array}{c}\Delta u_{\min }\\ \Delta u_{\min }\\ ⋮\\ \Delta u_{\min }\end{array}\right],\Delta U_{\max }=\left[\begin{array}{c}\Delta u_{\max }\\ \Delta u_{\max }\\ ⋮\\ \Delta u_{\max }\end{array}\right],{\widehat{P}}_{\min }^E=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right],{\widehat{P}}_{\max }^E=\left[\begin{array}{c}{P}_{\max }^E\\ P_{\max }^E\\ ⋮\\ P_{\max }^E\end{array}\right]\end{array}$, $\begin{array}{l}{x}_{\min }={\left[\begin{array}{ccccc}{w}_{g\min }& w_{r\min }& {\theta }_{\min }& y_{\min }& v_{t\min }\end{array}\right]}^T,u_{\min }={\left[\begin{array}{cc}{T}_{g\min }& {\beta }_{\min }\end{array}\right]}^T,\Delta u_{\min }={\left[\begin{array}{cc}\Delta T_{g\min }& \Delta {\beta }_{\min }\end{array}\right]}^T\\ x_{\max }={\left[\begin{array}{ccccc}{w}_{g\max }& w_{r\max }& {\theta }_{\max }& y_{\max }& v_{t\max }\end{array}\right]}^T,u_{\max }={\left[\begin{array}{cc}{T}_{g\max }& {\beta }_{\max }\end{array}\right]}^T,\Delta u_{\max }={\left[\begin{array}{cc}\Delta T_{g\max }& \Delta {\beta }_{\max }\end{array}\right]}^T\end{array}$.

The optimization problem (25) can be solved by using the interior point method, which is realized by using IPOT solver in MATLAB [33]. Then, the optimal solution of the control inputs (pitch angle and the generator torque) at the current time can be extracted and implemented.

4. Simulations Results

Three simulations are performed based on a typical 5-MW NREL VSWT system to validate the proposed NEMPC strategy. The corresponding parameters of the 5-MW NREL VSWT system are listed in

Table A1. The parameters in the 5-MW VSWT system.

Parameter	Symbol	Value
Number of blades	-	3
Rotor radius	-	63 m
Hub height	-	84.3 m
Hub inertia on the low-speed shaft	$J_H$	115,926 kg·m2
Blade inertia on the low-speed shaft	$J_B$	11,776,047 kg·m2
Generator inertia on high-speed shaft	$J_g$	534.116 kg·m2
Equivalent drive shaft torsion spring	$K_{\theta }$	867,637,000 N·m/rad
Equivalent drive shaft torsion damping	$B_{\theta }$	6,215,000 N·m/(rad/s)
Gearbox ratio	$N_g$	97
tower fore-aft inertia	$M_t$	430,000 kg
tower damping	$B_t$	17,600 N·s/m
tower stiffness	$K_t$	1,770,000 N/m
Air density	${\rho }_a$	1.2231 kg/m3
Rotor radius	$R_r$	63 m
Generator efficiency	${\eta }_g$	94.4%

in the Appendix A [32]. The corresponding constraints are shown in Table A2 [19,32]. In the simulations, the weights in the cost function (15)–(18) contain: ${\lambda }_1=20,{\lambda }_2=100,{\lambda }_3=5\ast {10}^8,{\lambda }_4=1,{\lambda }_5=0,{\lambda }_6=0.001$. The maximum generator speed, rotor speed, and generator torque are all arbitrarily set at 5% above their rated values, i.e.,$w_{g\max }=1.05w_{grated}$, $w_{r\max }=1.05w_{rrated}$, $T_{g\max }=1.05T_{grated}$. For the purpose of simplicity, the maximum electric power is chosen as $P_{\max }^E=P_{rated}^E$.

For comparison purposes, the classical NMPC strategy based on [19] is also constituted, which has the same model (13) and constraints (20)–(23) as the proposed NEMPC strategy in this paper. However, the tracking objective function (not economic cost function) is used in the classical NMPC, which can be described as:

$$l_{mpc}(x,u)={\displaystyle \sum _{i=1}^{N_P}[Q_{\Omega }{\left(w_r(k+i-1|k)-w_{rref}(w)\right)}^2+}{Q}_{\beta }{(\beta (k+i-1|k)-{\beta }_{ref}(w))}^2]$$

(26)

where $Q_{\Omega }=3.4\times {10}^7$, $Q_{\beta }=4.6\times {10}^8$ for the wind speed region II, $Q_{\Omega }=3.4\times {10}^5$, $Q_{\beta }=4.6\times {10}^5$ for wind speed region III, ${\beta }_{ref}(w)$ is the pitch angle reference and $w_{rref}(w)$ is the rotor speed reference, which can be derived as:

$$w_{rref}(w)=\{\begin{array}{lll}{w}_r^{\ast }(w)\hfill & if& w_r^{\ast }(w)\in [w_{rcut-in},w_{rrated}]\\ w_{rcut-in}\hfill & if& w_r^{\ast }(w)<w_{rcut-in}\\ w_{rrated}\hfill & if& w_r^{\ast }(w)>w_{rrated}\end{array},{\beta }_{ref}(w)=\{\begin{array}{lll}{\beta }^{\ast }(w)\hfill & if& w\le w_{rated}\\ {\beta }_{AR}(w)\hfill & if& w>w_{rated}\end{array}$$

(27)

${\beta }^{\ast }$ and ${\lambda }^{\ast }$ in Equation (27) are the optimal pitch angle and tip ratio in wind speed region II as demonstrated in Figure 1, can be derived via maximizing the power coefficient $C_P$:

$${\beta }^{\ast },{\lambda }^{\ast }=\underset{\beta ,\lambda }{\max }{C}_P(\beta ,\lambda )$$

(28)

Then, the optimal rotor speed $w_r^{\ast }$ can be obtained:

$$w_r^{\ast }(w)=\frac{{\lambda }^{\ast }w}{R}$$

(29)

${\beta }_{AR}(w)$ in Equation (27) can be solved from: $\frac{1}{2}\rho A{C}_P\left({\beta }_{AR}(w),\frac{R{w}_{rrated}}{w}\right)w^3=P_{rated}$.

The predictive/control horizon is set as $N_P=14$ for both classical NMPC and the proposed NEMPC in the following three groups of simulations.

4.1. Gradient Normal Wind

As indicated in Figure 5, a gradient normal wind varying between 7.5 m/s and 15 m/s is used in the simulation. The initial condition is set as $x_0={\left[\begin{array}{ccccc}93.54& 0.96& 1.85\times {10}^{-3}& 0.21& 0\end{array}\right]}^T, and t_0=0$. All the input variables, state variables, and output variables adopting both the proposed NEMPC and the classical NMPC strategies are shown as Figure 6a–h. During the time period from 0 s to 50 s, all the dynamic trajectories are almost the same when using the proposed NEMPC and the classical NMPC strategies, because the chosen initial state is just the optimal solution for both controllers under the initial 7.5 m/s normal wind speed.

At time instant 50 s, the normal wind speed steps from 7.5 m/s to 10 m/s, which is below the rated wind speed. Obviously, the VSWT system remains working in region II by using the classical NMPC strategy, which aims to follow the maximum power point and keeps the pitch angle zero. The proposed NEMPC strategy aims to attain a trade-off between the power maximization and turbine fatigue minimization by minimizing the economic cost function (19). Therefore, the pitch-angle is not always kept at zero, while the vibrations of the shaft torsion angle, the tower displacement, and the generation torque are much slighter when using the proposed NEMPC strategy, as shown in Figure 6c–f. No constraints are violated due to the constraints in Equations (20)–(23), which are imposed in both the classical NMPC strategy and the proposed NEMPC strategy, as seen in Figure 6.

At time instant 100 s, the normal wind speed steps from 10 m/s to 15 m/s, which is above the rated wind speed. Therefore, the working region of the VSWT system changes from region II to region III. The aim of the classical NMPC changes from the maximization power capturing to rated power capturing. The generator speed, rotor speed, generator torque, and generator power are changed to the rated values. Simultaneously, the pitch angle is enabled. As shown in Figure 6c–f, the vibration of the shaft torsional angle, the pitch angle, and especially the tower displacement, is much smaller when using the proposed NEMPC strategy than when using the classical NMPC strategy, due to $\theta$, $v_t$ and $\beta$ being considered in the economic cost function (19).

At time instant 150 s, the normal wind speed decreases from 15 m/s to 10 m/s. Thus, the wind turbine works back to region II. The pitch angle is disabled, the generator speed, rotor speed, and generator torque are performed to follow the maximum power point by using the classical NMPC strategy. From 6b–h, it is obvious that the vibration of each state variable is much more severe when using the classical NMPC strategy compared with the proposed NEMPC strategy. This is because the shaft load and the tower fatigue are not considered in the tracking objective function (26) for the classical NMPC strategy, while they are emphasized in the economic cost function (19) for the proposed NEMPC strategy. Vibrations in the shaft torsional angle, the tower displacement, and the tower displacement rate are harmful to wind turbine systems. Thus, the proposed NEMPC strategy is much better than the classical NMPC strategy in this respect.

From a more intuitive point of view, numerical comparisons have been made to illustrate the difference between these two controllers. Define the root mean square (RMS) values of mechanical loads as follows:

$$RM{S}_{shaft}=\sqrt{\frac{1}{Nsim}{\displaystyle \sum _{k=1}^{Nsim}{\theta }^2(k)}}$$

(30)

$$RM{S}_{tower}=\sqrt{\frac{1}{Nsim}{\displaystyle \sum _{k=1}^{Nsim}{v}_t^2(k)}}$$

(31)

where $Nsim$ is the length of the simulation, $RM{S}_{shaft}$ is the RMS value of gearbox load on the shaft, $RM{S}_{tower}$ is the RMS value of fatigue on the tower.

As shown in

Table 1. Data analysis based on the results in Figure 6.

	7.5 m/s	10 m/s	15 m/s	10 m/s
$AVG(P_g^{NEMPC})$ (MW)	1.4575	3.1873	4.9493	3.4317
$AVG(P_g^{NMPC})$ (MW)	1.4570	3.1714	4.9343	3.3323
$RM{S}_{shaft}^{NEMPC}$($\times {10}^{-3}$ rad)	1.820	3.288	4.784	3.416
$RM{S}_{shaft}^{NMPC}$ ($\times {10}^{-3}$ rad)	1.846	3.300	4.851	3.476
$RM{S}_{tower}^{NEMPC}$ ($\times {10}^{-2}$ m/s)	0	1.643	1.10	3.31
$RM{S}_{tower}^{NMPC}$ ($\times {10}^{-2}$ m/s)	0	4.483	76.58	16.44

, the average values of generator power $AVG\left(P_g\right)$, the RMS values of the gearbox load on shaft $RM{S}_{shaft}$, and the RMS values of fatigue on tower $RM{S}_{tower}$ by adopting both the proposed NEMPC strategy and the classical NMPC strategy at different normal wind speed are listed for comparison purposes.

From Table 1, it is obvious that the proposed NEMPC strategy enhances the average generator power over the classical NMPC strategy by 0.03%, 0.5%, 0.3%, and 2.98%, when the normal wind speeds are 7.5 m/s, 10 m/s, 15 m/s, and 10 m/s, respectively. During different normal wind speed periods, the proposed NEMPC strategy reduces the gearbox load on the shaft over the classical NMPC strategy by 1.41%, 0.36%, 1.38%, and 1.73%, respectively. Furthermore, it is obvious that the proposed NEMPC strategy reduces the tower fatigue a lot compared to the classical NMPC strategy from the RMS values of the fatigue on the tower, which proves that the fatigue on the tower can’t be ignored with regards to economic factors.

4.2. Stepwise Normal Wind

In general, 11.4 m/s represents the rated wind-speed, which is used to divide the wind speed into region II and III. Due to the different control objectives in regions II and III, a switching controller is widely used in the actual control of VSWT, which is a great challenge in controller design. To achieve the various operational requirements, the objective function of the classical NMPC strategy switches frequently according to wind speed operating regions. Thus, the rated wind speed 11.4 m/s is actually a very important speed in practical controller designs for VSWT systems.

The proposed NEMPC strategy in this paper adopts a unique economic cost function (19) which does not change according to the normal wind speed between the cut-in and cut-off wind speeds. Suppose the stepwise wind speed is initially set to 10.6 m/s, then steps to 12 m/s with an increment of 0.2 m within 25 s (from region II to region III), as shown in Figure 7a. As indicated in Figure 7b, it’s obvious that the pitch angle has a quick and small overshoot, before quickly returning to zero before t = 125 s. After 125 s, the fluctuations of the pitch angle accumulate around a certain optimal value when the proposed NEMPC strategy is used. At 125 s, the wind speed changes from 11.4 m/s to 11.6 m/s. This means that the special wind speed is around 11.4 m/s due to the economic objective function (19), above which the stable pitch angle no longer remains at zero. Thus, these simulation results prove that a wind speed of around 11.4 m/s was set to divide the region II and III, which is reasonable for considering the economic operation of the VSWT system.

4.3. Randomly Varying Wind Speed Based on the FAST Simulator

The NREL FAST is an inclusive aeroelastic simulator of two- and three-bladed horizontal-axis WTs [34]. To validate the practicality and validity of the proposed NEMPC strategy, the simulation is then tested on NREL FAST. In this group of simulations, a FAST subroutine was written into the MATLAB S-function, as indicated in Figure 8 and Figure 9.

The normal wind speed $w$ is always modeled as the combination of a slowly-changing mean wind-speed $w_m$ and a rapidly-changing turbulent wind-speed $w_d$:

$$w=w_m+w_d$$

(32)

where we suppose that $w_m$ follows the two-parameter Weibull distribution through a scale parameter $s1=10$ and shape parameter $s2=14$, and $w_d$ follows zero-mean Gaussian white noise distribution with standard variation $\sigma =1.899$ [28].

Then, a group of 600 s random normal wind speed data generated by the FAST simulator as indicated in Figure 10 is adopted in Case 3. The initial condition is set as $x_0={\left[\begin{array}{ccccc}99.77& 1.03& 2.1\times {10}^{-3}& 0.25& 0\end{array}\right]}^T, t_0=0$. Figure 11a–h indicate the performances of the generator speed, rotor speed, shaft torsion angle, tower displacement, tower displacement rate, pitch angle, generator torque, and generator power via the proposed NEMPC strategy.

From a more intuitive viewpoint, it can be observed that the performances of the proposed NEMPC strategy outperform the classical NMPC strategy. More specifically, the classical NMPC strategy shows more noticeable fluctuations in the performance of generator power, shaft torsion displacement, tower displacement, generator torque, and pitch angle when the wind speed has a sudden increase around its rated wind speed, i.e., 11.4 m/s, compared with the proposed NEMPC scheme.

Meanwhile, the classical NMPC strategy shows more successive fluctuations in the performance of generator power, shaft torsion displacement, tower displacement, generator torque, and pitch angle when the normal wind speed has a sudden drop around the rated wind speed, compared with the proposed NEMPC strategy.

As shown in Figure 11d–e, the tower displacement and the tower displacement by using the classical NMPC strategy fluctuates more severely according to the wind speed variation compared with the proposed NEMPC strategy, due to the economic index $l_{e2}$ related with the tower displacement considered in the economic cost function (19) for the proposed NEMPC. It is similar in the performances of the generator power, shaft torsion displacement, generator torque, and pitch angle, as shown in Figure 10. It can be concluded that the proposed NEMPC strategy decreases the operating costs while maintaining energy utilization, compared with the classical NMPC strategy. This is especially obvious when the wind speed varies around the rated wind speed. From

Table 2. Data analysis based on the results in Figure 11.

	9.7 m/s	9.4 m/s	10.8 m/s	9 m/s
$AVG(P_g^{NEMPC})$ (MW)	3.1162	2.9598	4.1231	2.5662
$AVG(P_g^{NMPC})$ (MW)	3.095	2.9291	4.1116	2.5221
$RM{S}_{shaft}^{NEMPC}$($\times {10}^{-3}$ rad)	3.234	3.092	3.913	2.581
$RM{S}_{shaft}^{NMPC}$ ($\times {10}^{-3}$ rad)	3.365	3.215	3.939	2.761
$RM{S}_{tower}^{NEMPC}$ ($\times {10}^{-2}$ m/s)	36.38	35.80	42.59	33.40
$RM{S}_{tower}^{NMPC}$ ($\times {10}^{-2}$ m/s)	39.44	35.83	43.22	32.86

, it is obvious that the proposed NEMPC strategy enhances the average generator power over the classical NMPC strategy by 0.67%, 1.05%, 0.28%, and 1.75%, when the average wind speeds are 9.7 m/s, 9.4 m/s, 10.8 m/s, and 9 m/s, respectively.

During different average wind speed periods, the proposed NEMPC strategy reduces the gearbox load on the shaft over the classical NMPC strategy by 4.05%, 3.96%, 0.68%, and 6.96%, respectively. Furthermore, it is obvious that the proposed NEMPC strategy reduces the tower fatigue a lot compared to the classical NMPC strategy from the RMS values of the fatigue on the tower. Therefore, this group of simulation-based FAST verifies the economy of the proposed NEMPC strategy.

It is well known that the computational burden of nonlinear MPC strategies increases exponentially according to the incensement of the predictive horizon. Various predictive horizons have been considered to investigate the computational burden by using classical NMPC and the proposed NEMPC strategies. The relative computational times of these two controllers are indicated in Figure 12. From Figure 12, it is obvious that the computational time of the proposed NEMPC is much smaller than the classical NMPC with the same predictive horizon. Simultaneously, the computational time of the classical NMPC increases more rapidly according to the incensement of the predictive horizon, compared with the proposed NEMPC strategy.

5. Conclusions

In this paper, the tower displacement and shaft torsional angle are taken into consideration in the development of a nonlinear model for a 5-MW VSWT system. Based on this comprehensive model, a NEMPC strategy for the VSWT system is proposed to increase the dynamic economy, i.e., through power generation maximization and fatigue load minimization during operation. Three groups of simulations for comparisons with the classical NMPC strategy demonstrate the effectiveness of the proposed NEMPC strategy. The simulation under random wind speeds based on the FAST simulator verified the reliability practicability of the proposed NEMPC strategy.

The contributions of this study compared with other current studies can be outlined as follows: Firstly, the tower dynamics are focused on in the modeling process for the VSWT with flexible connections on the shaft, and the tower displacement is considered as a part of the economic objective function, which can reduce tower fatigue. Secondly, all the related operational constraints and physical constraints in the VSWT system with flexible connections on the shaft are considered in the proposed NEMPC strategy in detail. Finally, simulations based on the FAST simulator, which is a realistic WT simulator for testing the system’s validity and practicability, were performed to show the practical implementation of the proposed controller using a turbulence wind speed profile. A future research topic is the design of a robust EMPC strategy taking into account the uncertainty of wind speed.