# Aeroservoelastic Pitch Control of Stall-Induced Flap/Lag Flutter of Wind Turbine Blade Section

Authors: Tingrui Liu
Journal: Shock and Vibration
 (2015)

---

## Abstract

The aim of this paper is to analyze aeroelastic stability, especially flutter suppression for aeroelastic instability. Effects of aeroservoelastic pitch control for flutter suppression on wind turbine blade section subjected to combined flap and lag motions are rarely studied. The work is dedicated to solving destructive flapwise and edgewise instability of stall-induced flutter of wind turbine blade by aeroservoelastic pitch control. The aeroelastic governing equations combine a flap/lag structural model and an unsteady nonlinear aerodynamic model. The nonlinear resulting equations are linearized by small perturbation about the equilibrium point. The instability characteristics of stall-induced flap/lag flutter are investigated. Pitch actuator is described by a second-order model. The aeroservoelastic control is analyzed by three types of optimal PID controllers, two types of fuzzy PID controllers, and neural network PID controllers. The fuzzy controllers are developed based on Sugeno model and intuition method with good results achieved. A single neuron PID control strategy with improved Hebb learning algorithm and a radial basic function neural network PID algorithm are applied and performed well in the range of extreme wind speeds.


---

## Body

## 1. Introduction

Stall flutter denotes unstable aeroelastic performance in stall conditions. Stall-induced nonlinear flutter of wind turbine blade, as a typical nonlinear aeroelastic instable phenomenon, is an important reason of fatigue damage for wind turbine. How to effectively avoid stall-induced aeroelastic instability has become an important research needed to be investigated. Moreover, particular concepts and technologies that can alleviate fatigue loads should be investigated. Meanwhile the simplified stall-induced flap/lag flutter investigation of typical blade section plays an important role in this area due to its simplicity and convenience, so in this study, researches on aeroelastic instability and aeroservoelastic control will be depicted based on typical blade section.In the past ten years, a number of issues related to the modeling, nonlinear vibration, and nonlinear stall flutter are investigated. Structure modeling and aeroelastic stability analysis of individual blade section subjected to combined flap/lag motion are investigated [1]. Considering a wide ranging parametric variation, Sunetra and Hester [2] investigate effects on the aeroelastic instability and overall nonlinear dynamical response of the flap/lag system. Numerical simulation and local bifurcation analysis of an airfoil section with a single degree of freedom subject to quasi-steady aerodynamic forces in the stochastic inflow are studied by Hansen [3, 4] in stall-induced vibration process.Recently most issues related to aeroservoelastic control today are focused on flexible aircraft wing or rotorcraft blade [5]. A single-input and single-output controller and a multi-input and multioutput controller are designed based on the unconstrained optimization techniques to demonstrate the effects of nonlinearity on the linear flutter control system [6]. A nonlinear control law designed to limit the control surface pitch response was found to increase the flutter speed considerably by forcing the system to undergo limit cycle oscillations instead of fluttering [7]. Coupling the structural state equations with the aerodynamic state equations of the wing and the control surface based on the reduced order models, a transonic aeroservoelastic model in state-space is constructed and the suboptimal control method based on output feedback is used to design the flutter suppressing law [8].However aeroelastic control in field of wind turbine blade is rarely studied. In this study, aeroservoelastic control for wind turbine blade is implemented by pitch angle controller. Pitch control is the most efficient and popular power control method, especially for variable-speed wind turbines. Today most of pitch control processes for wind turbine focus on power conversion, wind energy utilization, and multivariable and nonlinear control [9–12]. Reliable and powerful control strategies are needed for wind energy conversion systems to achieve maximum performance [9]. Lin et al. [10] present the design of a fuzzy sliding mode loss-minimization control for the speed of a permanent magnet synchronous generator and a high-performance on-line training radial basis function network for the turbine pitch angle control. Duong et al. [11] present a hybrid controller based on PI and fuzzy technique for the pitch angle controller which has been one of the most common methods for smoothing output power fluctuations in order to improve power quality and maintain the stable output generated from wind farm. Yilmaz and Özer [12] propose an artificial neural network pitch angle controller based on two types of artificial neural network methods for wind turbines.However pitch control for aeroelastic instability of wind turbine blade in stall situation is more rarely studied. Kallesøe [13] studies Beddoes-Leishman aerodynamic model and presents an aeroservoelastic pitch control strategy based on conventional PID controller and quadratic regulator (QR) for blade load-reducing control process. Performance of QR regulation depends on the choice of weighting matrices which have no analytic solution. So in some sense, this optimal control in fact is entirely man-made.In this paper, the emphasis of pitch control is not on power and efficiency but on improving flutter suppression effect for destructive flapwise and edgewise motions of stall-induced flutter for wind turbine blade in large angle of attack behavior. An effort is to solve the instability problem of the stall flutter of rotating blade section combined flap/lag motion. Several intelligent control techniques have been applied to improve the quality of nonlinear aeroelastic stability. Based on nonlinear flutter calculations carried out using Fourier analysis to extract the harmonics from the nonlinear ONERA aerodynamics [14] suitable for large angle of attack, a set of unsteady aerodynamic equations suitable for pure pitching motion are applied. The resulting nonlinear aeroelastic formulations are linearized to analyze aeroelastic stability. The influences of some related factors, rotating speed and wind speed, are investigated. The paper gives methods of aeroservoelastic pitch control based on optimal PID control, two types of fuzzy PID control, a single neuron PID control strategy with improved Hebb learning algorithm, and a radial basic function neural network algorithm, which can determine flutter suppression of the blade section.


## 2. Nonlinear Aeroelastic Model of the Flap/Lag System

The stall-induced flap/lag flutter investigation is simplified and based on the instability analysis of a special rotating blade section located at a distancer from the hub in Figure 1. The structure is hinged in such a way that its motions have four independent translational degrees of freedom: an edgewise direction (lag), denoted by y; one perpendicular to that direction (flap), denoted by z; pitch angle φ; and rotating angular frequency ϕ. For the specific designed section, it is assumed that the distance between the centre of gravity and the elastic axis is negligible; the distance between the aerodynamic centre and the elastic axis is negligible; the elastic twist of the section might be negligible.Figure 1
Coordinate system and aerodynamic force.The structural equations are derived using Lagrange’s equations. The kinetic energyT and potential energy V for the system are expressed as(1a)
T
=
1
2
I
C
G
φ
˙
2
+
1
2
M
r
2
ϕ
˙
2
+
1
2
ρ
b
·
y
˙
cos
⁡
φ
-
z
˙
sin
φ
-
φ
˙
y
sin
φ
+
z
cos
⁡
φ
+
r
ϕ
˙
2
+
1
2
ρ
b
z
˙
cos
⁡
φ
+
y
˙
sin
φ
+
φ
˙
(
y
cos
⁡
φ
-
z
sin
φ
)
2
(1b)
V
=
1
2
ω
y
2
ρ
b
y
2
+
ω
z
2
ρ
b
z
2
,where I
C
G is the rotational inertia of per unit spanwise length; the inertia of the rotor, gearbox, and generator is described by the mass M; the mass per length of blade section is described by ρ
b; the natural angular frequencies of motions in y and z directions are described by ω
y and ω
z, respectively.Using Lagrange’s equations, the equations of motions inz and y directions are deduced to(2a)
z
¨
+
2
ξ
z
ω
z
z
˙
+
ω
z
2
z
+
φ
¨
y
-
r
ϕ
¨
sin
φ
+
2
φ
˙
y
˙
+
φ
˙
2
z
=
ρ
b
-
1
(
L
C
cos
⁡
ψ
+
L
N
C
cos
⁡
φ
+
D
sin
ψ
)
(2b)
y
¨
+
2
ξ
y
ω
y
y
˙
+
ω
y
2
y
-
φ
¨
z
+
r
ϕ
¨
cos
⁡
φ
-
2
φ
˙
z
˙
-
φ
˙
2
y
=
ρ
b
-
1
D
cos
⁡
ψ
-
L
N
C
sin
φ
-
L
C
sin
ψ
,where ξ
y and ξ
z are the damping ratios in y and z directions, respectively. The right-hand side terms in equations are aerodynamic expressions.Based on variable parameters, such as pitch angleφ and wind velocity U, a series of Fourier analyses were performed on the lift and drag coefficients for pure pitching motion of the blade at its quarter-chord aerodynamic center with constant incoming velocity V
0 [14]. A set of extracted nonlinear ONERA aerodynamic model equations suitable for this pure pitching motion are as follows [15].Lift expressions are(3a)
L
N
C
=
1
2
ρ
S
L
(
s
L
b
V
0
φ
˙
+
l
L
b
V
0
φ
˙
+
k
L
b
2
φ
¨
)
(3b)
L
C
=
1
2
ρ
S
L
V
0
2
C
1
L
+
V
0
2
C
2
L
,where C
1
L and C
2
L are nonlinear lift coefficients; its expressions are(3c)
C
˙
1
L
+
λ
L
V
0
b
C
1
L
=
λ
L
a
0
L
V
0
b
φ
+
λ
L
σ
L
φ
˙
+
γ
L
a
0
L
φ
˙
+
γ
L
σ
L
b
V
0
φ
¨
(3d)
C
¨
2
L
+
a
L
V
0
b
C
˙
2
L
+
r
L
V
0
2
b
2
C
2
L
=
-
r
L
V
0
2
b
2
Δ
C
L
+
e
L
V
0
b
∂
Δ
C
L
∂
t
.Drag expression is(3e)
D
=
1
2
ρ
c
V
0
2
C
D
1
+
V
0
2
C
D
2
,where C
D
2 is nonlinear drag coefficient; its expression is(3f)
C
¨
D
2
+
a
D
V
0
b
C
˙
D
2
+
r
D
V
0
2
b
2
C
D
2
=
-
r
D
V
0
2
b
2
Δ
C
D
+
e
D
V
0
b
φ
˙
.Herein, nonlinear parts in nonlinear items and other aerodynamic parameters are expressed as (4)
Δ
C
L
=
a
x
l
ψ
-
φ
-
0.1396
a
L
1
-
0.3142
a
L
2
,
a
x
l
=
0
,
a
L
1
=
0
,
a
L
2
=
0
,
ψ
-
φ
≤
0.1396
a
x
l
=
a
L
1
,
a
L
1
=
6.32284
,
a
L
2
=
0
,
0.1396
<
ψ
-
φ
≤
0.3142
a
x
l
=
a
L
1
+
a
L
2
,
a
L
1
=
6.32284
,
a
L
2
=
-
0.42284
,
i
ψ
-
φ
>
0.3142
S
L
=
c
,
s
L
=
π
,
k
L
=
0.5
π
,
a
0
L
=
σ
L
=
5.9
,
λ
L
=
0.15
,
γ
L
=
0.55
,
l
L
=
0
,
e
=
c
4
,
a
L
=
0.25
+
0.4
Δ
C
L
2
,
r
L
=
0.2
+
0.23
Δ
C
L
2
2
,
e
L
=
-
2.7
Δ
C
L
2
;
V
0
=
U
2
+
r
Ω
2
,
Δ
C
D
=
-
a
D
1
α
-
a
D
2
α
2
-
a
D
3
α
3
,
a
D
=
0.32
,
a
D
1
=
0.42
,
a
D
2
=
0.1437
,
a
D
3
=
4.923
,
C
D
1
=
0.014
,
e
D
=
-
0.015
Δ
C
L
2
,
r
D
=
0.2
+
0.1
Δ
C
L
2
2
.Note that these equations have been written in terms of the force coefficients C
L, C
D rather than the circulation coefficients in original ONERA model. Each of Δ
C
L and Δ
C
D represents the deviation of the linear static curve from the quasistatic curve for the lift and drag, respectively.


## 3. System Linearization and Stability Analysis

To solve the stability analysis of the nonlinear aeroelastic system as given by (2a)–(3f), the following steps will be implemented. Inserting (3a)-(3b) and (3e) into (2a) and (2b), considering (2a) and (2b) are in conjunction with the nonlinear items of (3c)-(3d) and (3f), and assuming(5)
X
=
z
y
C
1
L
C
2
L
C
D
2
φ
ϕ
Tand ignoring the items of the first-order reciprocal product of variables according to linearizing process analysis of nonlinear aeroelastic stability in [15] result in the nonlinear equations governing the motions of the aeroelastic system. It is matrix equations with subequation structures as follows:(6)
M
M
X
¨
+
C
M
X
˙
+
K
M
X
=
Q
M
,where the 5 × 7 coefficient matrices M
M, C
M, and K
M and 5 × 1 matrix Q
M are displayed in Appendix A.Towards the goal of stability analysis, (6) is linearized. The considered numerical procedures involve the linearizing process and its postprocessing. Firstly this is done by assuming small oscillations around the equilibrium position of the equation; that is, assuming X
=
X
0
+
ε
X
1, especially  ϕ
=
ϕ
0
+
ε
ϕ
1
=
Ω
t
+
ε
ϕ
1, and inserting it into (6) give(7)
M
M
0
+
ε
M
M
1
ε
X
¨
1
+
C
M
0
+
ε
C
M
1
ε
X
˙
1
+
K
M
0
+
ε
K
M
1
(
X
0
+
ε
X
1
)
=
Q
M
0
+
ε
Q
M
1
.Balancing terms ofε
0 and ε
1 leads to the equilibrium position equations and the linear approximation equations, respectively, which are shown in (8)-(9). The static equilibrium position of (8) comprises a set of nonlinear algebraic equations for the unknown state vector X
0. The linear dynamic perturbation of (9) about a static equilibrium sate is depending on the wind velocity, pitch angle, and rotating speed, which can be illustrated as(8)
K
M
0
X
0
=
Q
M
0
(9)
M
M
0
X
¨
1
+
C
M
0
X
˙
1
+
K
M
0
X
1
=
Q
M
1
-
K
M
1
X
0
,where (10)
X
0
=
z
0
y
0
C
1
L
0
C
2
L
0
C
D
20
φ
0
Ω
T
,
X
1
=
z
1
y
1
C
1
L
1
C
2
L
1
C
D
21
φ
1
ϕ
1
T
.Herein, by matrix transformation, the item on right-hand side of (9) can be expressed as(11)
Q
M
1
-
K
M
1
X
0
=
Q
W
X
1and the 5 × 7 coefficient matrices K
M
0, M
M
0, and C
M
0 and 5 × 1 matrix Q
M
0, as well as Q
W, are displayed in Appendix B. So after matrix transformation, (9) is rewritten as(12)
M
M
0
X
¨
1
+
C
M
0
X
˙
1
+
(
K
M
0
-
Q
W
)
X
1
=
0
.As for solution of the static equilibrium position of (8), because of existence of stall, the equation is a complicated function itself. Hence in order to get the exact solution, a secant method [16] with an appropriate initial guess should be used. For example, the initial values can be assigned by solving the quadratic equation in the linear unstalled region, or in the nonlinear stalled region [14].In addition, to determine the time domain response of this dynamic perturbation of (12), a Runge-Kutta time-marching approach is applied. Using the method presented in [17], (12) will be expressed in state-space form. In general, upon defining the state vector Y
=
[
X
1
T
,
X
˙
1
T
]
T and adjoining the identity equation X
˙
1
=
X
˙
1, (12) can be converted to the first-order linear differential equation as(13)
Y
˙
t
=
A
Y
t
,where the 14 × 14 state matrix A is given by(14)
A
=
0
I
E
-
M
M
0
-
1
(
K
M
0
-
Q
W
)
-
M
M
0
-
1
C
M
0
.Herein, I
E is the unitary matrix.Generally system stability can be demonstrated by eigenvalue analysis of the dynamic perturbation in (13). However for this aeroelastic system of stall-induced flap/lag flutter, aerodynamic eigenvalues and structural eigenvalues cannot easily be distinguished from each other [15]; according to Lyapunov’s first linearization method, the equilibrium of the nonlinear system is inconclusive for LHP eigenvalues and one or more eigenvalues on the imaginary axis. Aeroelastic stability of stall-induced flutter can be confirmed by linearized time response of the dynamic perturbation by time-marching approach [15].The validity of the method of stability analysis is verified in [15]. Some cases intended to highlight the effects played by rotating speed and wind velocity on the vibration and aeroelastic stability of rotating blade will be presented. Basic structure parameters are c
=
0.05 m, b
=
c
/
2, ρ
b
=
15 kg m−1, ω
z
=
8 rad s−1, ω
y
=
12 rad s−1, ξ
z
=
0.01, ξ
y
=
0.02, and r
=
1 m.


### 3.1. Effect of Rotating Speed on Stability

Figure2 shows linearized time responses of flap and lag motions for rotating speed Ω from 30 rpm to 120 rpm at intervals of 30 rpm, with wind velocity U
=
5 m s−1. As far as the response tendency is concerned, aeroelastic stability of flap motion decreases apparently with the increase of rotating speed, but rotating speed has little influence on stability of lag motion. When Ω = 30 rpm and Ω = 60 rpm the flap motions are stable and convergent; when Ω = 90 rpm and Ω = 120 rpm the flap motions are unstable with divergent tendency. Figure 2(b) shows that there are no obvious changes in edgewise stability with the increase of Ω. Meanwhile, all the trends are converging.Figure 2
Flap and lag displacements forΩfrom 30 rpm to 120 rpm at intervals of 30 rpm.
(a)
Flap displacements forΩfrom 30 rpm to 120 rpm at intervals of 30 rpm
(b)
Lag displacements forΩfrom 30 rpm to 120 rpm at intervals of 30 rpm


### 3.2. Effect of Wind Velocity on Stability

Figures3(a)-3(b) show the linearized time responses of flap and lag displacements for wind velocity U from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1, with rotating speed Ω = 60 rpm. It is obvious from Figure 3(a) that the flap displacement for U
=
5 m s−1 is convergent and stable and for other cases of U
=
10 m s−1, U
=
15 m s−1, and U
=
20 m s−1 flap displacements are divergent. It can be seen from Figure 3(b) that lag displacements for U
=
5 m s−1, U
=
10 m s−1, and U
=
15 m s−1 are convergent and for another case of U
=
20 m s−1 lag displacement is unstable (flutter and divergence). From the overall tendency, aeroelastic stability is getting worse with wind speed increasing; furthermore when U
>
20 m s−1, the flutter and divergence intensify.Figure 3
Flap and lag displacements forU from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1 and critical flutter velocities of the section at given rotating speeds.
(a)
Flap displacements forU from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1
(b)
Lag displacements forU from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1
(c)
The critical flutter velocities versus rotating speedsΩ from 30 rpm to 120 rpm at intervals of 15 rpmIn addition, Figure3(c) shows plots of the critical flutter velocities of the blade section versus rotating speeds from 30 rpm to 120 rpm at intervals of 15 rpm. For the case of Ω = 60 rpm in Figure 3(c), the critical flutter velocity U
=
8.5 m s−1, which is precisely between values of 5 m s−1 (convergent) and 10 m s−1 (divergent) in Figure 3(a).


## 3.1. Effect of Rotating Speed on Stability

Figure2 shows linearized time responses of flap and lag motions for rotating speed Ω from 30 rpm to 120 rpm at intervals of 30 rpm, with wind velocity U
=
5 m s−1. As far as the response tendency is concerned, aeroelastic stability of flap motion decreases apparently with the increase of rotating speed, but rotating speed has little influence on stability of lag motion. When Ω = 30 rpm and Ω = 60 rpm the flap motions are stable and convergent; when Ω = 90 rpm and Ω = 120 rpm the flap motions are unstable with divergent tendency. Figure 2(b) shows that there are no obvious changes in edgewise stability with the increase of Ω. Meanwhile, all the trends are converging.Figure 2
Flap and lag displacements forΩfrom 30 rpm to 120 rpm at intervals of 30 rpm.
(a)
Flap displacements forΩfrom 30 rpm to 120 rpm at intervals of 30 rpm
(b)
Lag displacements forΩfrom 30 rpm to 120 rpm at intervals of 30 rpm


## 3.2. Effect of Wind Velocity on Stability

Figures3(a)-3(b) show the linearized time responses of flap and lag displacements for wind velocity U from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1, with rotating speed Ω = 60 rpm. It is obvious from Figure 3(a) that the flap displacement for U
=
5 m s−1 is convergent and stable and for other cases of U
=
10 m s−1, U
=
15 m s−1, and U
=
20 m s−1 flap displacements are divergent. It can be seen from Figure 3(b) that lag displacements for U
=
5 m s−1, U
=
10 m s−1, and U
=
15 m s−1 are convergent and for another case of U
=
20 m s−1 lag displacement is unstable (flutter and divergence). From the overall tendency, aeroelastic stability is getting worse with wind speed increasing; furthermore when U
>
20 m s−1, the flutter and divergence intensify.Figure 3
Flap and lag displacements forU from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1 and critical flutter velocities of the section at given rotating speeds.
(a)
Flap displacements forU from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1
(b)
Lag displacements forU from 5 m s−1 to 20 m s−1 at intervals of 5 m s−1
(c)
The critical flutter velocities versus rotating speedsΩ from 30 rpm to 120 rpm at intervals of 15 rpmIn addition, Figure3(c) shows plots of the critical flutter velocities of the blade section versus rotating speeds from 30 rpm to 120 rpm at intervals of 15 rpm. For the case of Ω = 60 rpm in Figure 3(c), the critical flutter velocity U
=
8.5 m s−1, which is precisely between values of 5 m s−1 (convergent) and 10 m s−1 (divergent) in Figure 3(a).


## 4. Aeroservoelastic Control

The rotating rotor blade has controllable blade pitch angle; a controller requests a pitch angle and a pitch actuator that can adjust the actual pitch angle. Here the pitch actuator is described by a second-order model:(15)
I
φ
φ
¨
+
d
φ
φ
˙
+
k
φ
φ
-
k
φ
φ
ref
=
0
,where I
φ
=
0.2, d
φ
=
1.13, and k
φ = 1 are the moment of inertia, damping, and stiffness fitted to give a response similar to a pitch actuator; φ
ref is the pitch angle requested by the controller. Due to the fact that blade itself used here is not a standard blade, an overdamping actuator is especially selected to test the effect of actuation under adverse conditions. Although for a mechanical actuator the actuation based on these parameters looks appropriate, it might be realized by a piezoelectric actuator.To design and evaluate controllers, the control objectives and the control cost are quantified by measuring the rotor speed error Δϕ. Most rotor blades today are controlled by PID controllers. The PID controllers implemented in this study use the approximate rotor speed error Δϕ to establish the control signal:(16)
φ
ref
=
K
P
Δ
ϕ
+
K
I
∫
0
t
Δ
ϕ
d
t
+
K
D
d
Δ
ϕ
d
t
=
K
P
ϕ
˙
-
Ω
+
K
I
∫
0
t
ϕ
˙
-
Ω
d
t
+
K
D
d
ϕ
˙
-
Ω
d
t
=
K
P
(
ε
ϕ
˙
1
)
+
K
I
∫
0
t
ε
ϕ
˙
1
d
t
+
K
D
d
ε
ϕ
˙
1
d
t
,where K
P, K
I, and K
D are proportional gain, integral gain, and derivative gain, respectively, all being tuned to optimize control performance. Let φ
=
φ
0
+
ε
φ
1, and inserting (16) into (15), the linearized second-order model can be expressed by(17)
I
φ
k
φ
φ
¨
1
-
K
D
ϕ
¨
1
+
d
φ
k
φ
φ
˙
+
k
φ
φ
-
K
P
ϕ
˙
1
+
φ
1
-
K
I
ϕ
1
=
0
.


### 4.1. The Linearized Aeroservoelastic Model and PID Control Strategy

Inserting (17) into (12), the linearized aeroservoelastic model can be rewritten as(18)
M
X
¨
1
+
C
X
˙
1
+
K
X
1
=
0
,where the M, C, and K are 6 × 7 coefficient matrices, and furthermore, (18) can be transformed as(19)
X
¨
1
=
A
M
2
X
˙
1
+
A
M
1
X
1
+
W
P
K
P
X
˙
1
+
K
I
∫
X
˙
1
d
t
+
K
D
d
X
˙
1
d
t
,where the coefficient matrices A
M
2, A
M
1, and W
P are displayed in Appendix B. The block diagram of (19) in Simulink is shown in Figure 4.Figure 4
Block diagram of PID control strategy in Simulink.The PID parameters,K
P, K
I, and K
D are optimized by optimal control. The optimal control is under certain concrete conditions, to achieve special control object and to make selected target maximal or least. From the optimal control, the PID parameters can be found during optimization. The target function used here adopts ITAE rule, defined as(20)
J
ITAE
=
∫
0
∝
t
e
(
t
)
d
t
.When time t is relatively big, to assure small target value, the steady error must be diminished, so as to make system speed in steady area. Each parameter in target transfer function may be evaluated and initialized.The Matlab functions are applied to search for the optimal parameters for the system [18]. It can be seen from the practice that the program is quite versatile in finding the optimal controllers. However, in some disturbed states, it may not find a solution due to unstable controlled plant. This can be a drawback in conventional optimization algorithms, but such problems can be avoided by multiple choices of initial values based on genetic algorithm [19].In order to improve the computational speed, the optional PI controller (K
D
=
0) or PD controller (K
I
=
0) is applied. Throughout the numerical illustrations, the aeroelastic instability case of U
=
30 m s−1 and Ω
=
60 rpm is considered with aeroservoelastic control implemented in Figure 5.Figure 5
Aeroelastic instability and aeroservoelastic control case ofU
=
30 m s−1 and Ω = 60 rpm.
(a)
Flap and lag displacements for the case ofU
=
30 m s−1 and Ω = 60 rpm
(b)
Lag displacements for the three optimal aeroservoelastic controllers
(c)
Flap displacements for the three optimal aeroservoelastic controllersFigure5(a) shows that both flap and lag motions are divergent (instable): flapwise displacement is positively divergent; the lag one is negatively divergent. Figures 5(b)-5(c) depict the aeroservoelastic control results of stall-induced flap/lag flutter for the three optimal aeroservoelastic controllers mentioned above.It demonstrates the obviously improved aeroelastic stability with rapidly convergent tendency for both flap and lag motions. With the change of time, the flutter amplitude decreases rapidly and gradually tends to be zero. The effect of flutter suppression is obvious. In contrast to flap case, effect of aeroservoelastic control on lag case is more prominent, and furthermore PID and PD controllers are better means of aeroservoelastic control. As for aeroservoelastic control in flapwise direction, PI controller is a good choice.


### 4.2. Fuzzy PID Control Strategy

Due to the fact that the controller in Figure4 is based on the full state feedback including the aerodynamic variables, the realization of this full state feedback controller is difficult in practice. So the simulation in Figure 4 is only the theoretical research and an object of reference in order to make comparison with other physically achievable feedback controllers for structural or mechanical variables such as fuzzy controllers and neural network controllers below.A self-tuning fuzzy PI (FPI) PID controller is designed, which can automatically tune the parameters (K
P and K
I) of PI controller. The advantage of the fuzzy controller here is that it is physically achievable feedback controllers for mechanical variable ϕ
˙
1, compared with optimal controllers. During real-time control, the self-tuning process is carried out by PLC instead of optimal control, which can also avoid large amount of calculation of genetic algorithm in the process of optimization in some disturbed states.Firstly rewrite (19) as(21)
X
¨
1
=
A
M
2
X
˙
1
+
A
M
1
X
1
+
W
P
1
K
P
ϕ
˙
1
+
K
I
∫
ϕ
˙
1
d
t
+
K
D
d
ϕ
˙
1
d
t
,where W
P
1
=
M
-
1
[
0
0
0
0
0
1
]
T. The block diagram of FPI PID control in Simulink is illustrated in Figure 6.Figure 6
The block diagram of FPI PID strategy.Fuzzy control depends on the fuzzy algorithm between the information of process and control input. Fuzzy controllers from their inception have demonstrated a vast range of applicability to processes where the plant transfer function is not defined but the control action can be described in terms of linguistic variables and are often being used to improve the performance of a system where the plant transfer function is known [20]. For this physically achievable feedback system, FPI controller is devised. Two input variables, error (E) and change in error (E
D), and two output variables K
P and K
I with five linguistic variables of membership function are used. The input linguistic variables aremf1 (negative medium),mf2 (negative small),mf3 (zero),mf4 (positive small), andmf5 (positive medium). Five triangular-membership function forms for E and five triangular-membership function forms forED(de/dt) are determined the same way, which are shown in Figure 7. Borders of both function sequences vary between +3 and −3. In this study, a fuzzy rule type of the designed FPI controller is developed based on Sugeno model in Simulink. So there are 25 weight values in both K
P and K
I designs. According to intuition method, list of linguistic rules is shown in Tables 1-2.Table 1
Rules forK
p.
E∖ED
m
f
1
m
f
2
m
f
3
m
f
4
m
f
5
m
f
1
1
0.85
0.65
0.5
0.5
m
f
2
0.85
0.85
0.65
0.5
0.5
m
f
3
0.85
0.65
0.5
0.35
0.15
m
f
4
0.65
0.5
0.35
0.35
0.15
m
f
5
0.65
0.5
0.35
0.15
0.15Table 2
Rules forK
I.
E∖ED
m
f
1
m
f
2
m
f
3
m
f
4
m
f
5
m
f
1
0
0.5
1
1
1.5
m
f
2
0.5
1
1
1.5
2
m
f
3
0.5
1
1.5
2
2.5
m
f
4
1
1.5
2
2
2.5
m
f
5
1.5
2
2
2.5
3Figure 7
Membership function plots for both function sequences.In addition, a self-tuning fuzzy PD (FPD) PID controller is designed. Its main advantage is that FPD controller parameters can be changed very quickly in response to changes of the system dynamics. Five triangular-membership function forms forE and E
D are determined the same way as Figure 7. The block diagram of FPD PID control strategy is illustrated in Figure 8. According to intuition method, list of linguistic fuzzy rules is shown in Table 3. Simulation is performed by substituting structure in Figure 6 for corresponding one of Figure 8.Table 3
Rule weight values of FPD PID controller.
E∖ED
m
f
1
m
f
2
m
f
3
m
f
4
m
f
5
m
f
1
−16
−12.8
−2.4
−8
0
m
f
2
−13.6
−12
−6.4
0
11.2
m
f
3
−14.4
−8
0
10.4
16
m
f
4
−9.6
0
13.6
12.8
15.2
m
f
5
0
8
12
15.2
16Figure 8
The block diagram of FPD PID strategy.In order to make qualitative analysis for the effectiveness of the fuzzy PID (including FPI and FPD) control algorithm, some fixed speeds in different ranges of wind speeds are selected.Case One. Consider U
=
5 m s−1 and U
=
10 m s−1, respectively, in low range of wind speeds for less than 17 m s−1.Case Two. Consider U
=
18 m s−1, an approximate intermediate value of wind velocity for the choice of the two types of fuzzy PID controllers.Case Three. Consider U
=
20 m s−1 and U
=
24 m s−1, respectively, in high range of wind speeds for more than 19 m s−1.Figures9(a)-9(b) depict the linearized responses of flap and lag displacements of Case one for different fuzzy controllers, which shows obvious superiority in flapwise stability of FPD PID controller; Figure 10 depicts the linearized responses of flap and lag displacements of Case two for the intermediate wind velocity of U
=
18 m s−1, which shows almost equal stability of both FPI and FPD PID controllers; Figures 11(a)-11(b) depict the linearized responses of flap and lag displacements of Case three for different fuzzy controllers, which shows obvious superiority in stability of FPI PID controller. It can be seen in the process of qualitative analysis that in low range of wind speeds FPD PID controller has the superiority compared with FPI PID controller. However simulation result indicates that the FPI PID controller is more effective and predominant than FPD controller in high range of wind speeds. Here the rotating speed Ω = 60 rpm.Figure 9
Flap and lag displacements of Case one for different fuzzy controller.
(a)
Displacements forU
=
5 m s−1
(b)
Displacements forU
=
10 m s−1Figure 10
Flap and lag displacements of Case two forU
=
18 m s−1.Figure 11
Flap and lag displacements of Case three for different fuzzy controllers.
(a)
Displacements forU
=
20 m s−1
(b)
Displacements forU
=
24 m s−1It should be stated that the intermediate value of wind velocity for the choice of the two types of fuzzy PID controller varies with the rotating speed. For example, the intermediate value of wind velocity isU
=
24.5 m s−1 for rotating speed Ω = 120 rpm.


### 4.3. Single Neuron PID Control Strategy with Improved Hebb Learning Algorithm

In general, under the condition of normal wind speed in the range ofU
<
32 m s−1, as mentioned above, the FPI PID control can completely meet the requirements. However in some special cases of extreme wind speeds, such as more than 50 m s−1, FPI algorithm might fail because of highly complicated and nonlinear instability. A single neuron (SN) PID control with improved Hebb learning algorithm can make up for the defect from the aspects of numerical simulation (due to the lack of real-time operation environment in the range of U
>
50 m s−1). With the interdisciplinary development of modern computational technologies and biological science, artificial neural networks, as typical artificial intelligence models and algorithms, have become an attractive approach for modeling highly complicated and nonlinear system without understanding the nature of phenomenon [21].Using the improved Hebb learning algorithm, update rules of the three weights can be written as [22, 23](22)
ω
1
(
k
)
=
ω
1
(
k
-
1
)
+
η
P
e
k
u
k
e
k
-
Δ
e
k
ω
2
(
k
)
=
ω
2
(
k
-
1
)
+
η
I
e
(
k
)
u
(
k
)
[
e
(
k
)
-
Δ
e
(
k
)
]
ω
3
k
=
ω
3
k
-
1
+
η
D
e
k
u
k
e
k
-
Δ
e
k
,where η
P, η
I, and η
D are adaptive proportional, integral, and differential learning rates, respectively. The three weights are the system state variables, so the control law is written as(23)
u
k
=
u
k
-
1
+
K
∑
i
=
1
3
ω
i
0
k
x
i
k
,where the normalized weights and the three parameters are(24)
ω
i
0
k
=
ω
i
k
∑
i
=
1
3
ω
3
k
,
x
1
k
=
e
k
,
x
2
k
=
Δ
e
k
=
e
k
-
e
k
-
1
,
x
3
k
=
Δ
2
e
k
=
e
k
-
2
e
k
-
1
+
e
k
-
2
.Summarizing the above algorithm, the SN PID control block diagram with algorithmS-function in Matlab and the system simulation diagram are illustrated in Figure 12.Figure 12
SN PID control block diagram and the system simulation diagram.
(a)
SN PID controller
(b)
System simulation diagramFigures13(a)–13(c) show three simulation results of flap and lag displacements of SN PID controller for U
=
30 m s−1 and Ω = 60 rpm. Figure 13(d) shows flap and lag displacements of FPI PID controller. It can be seen from Figures 13(a)–13(c) that the three simulation results including vibration amplitudes and steady-state values are quite different although system stability is the same and convergent. Among them Figure 13(c) is the most effective from the magnitude of vibration amplitudes of view. Figures 13(c)-13(d) show the difference between the third simulation of SN PID controller and simulation of FPI PID controller. It is seen that the flutter amplitudes of FPI PID controller in both flapwise and edgewise directions are smaller than those of SN PID controller.Figure 13
Flap and lag displacements of SN PID control (a)–(c) and FPI PID control (d) forU
=
30 m s−1 and Ω
=
60 rpm.
(a)
The first simulation of SN PID
(b)
The second simulation of SN PID
(c)
The third simulation of SN PID
(d)
Simulation of FPI PIDThe advantages of SN PID controller over FPI PID controller lie in its performance in the range of extreme wind speeds, such as more than 50 m s−1. Figures 14(a)-14(b) show the flap and lag displacements of SN PID control for U
=
50 m s−1 and U
=
55 m s−1, respectively. Although there might be some deviations from zero in flapwise motions, the control effects of SN PID controller are satisfactory within the range of extreme wind speeds where the simulation algorithm of FPI PID controller is not feasible because it is unable to reduce the step size without violating the minimum step size for a certain number of consecutive times during simulation. This might be caused by dependent dynamic states or by the high stiffness of the system.Figure 14
Flap and lag displacements of SN PID control forU
=
50 m s−1 and U
=
55 m s−1, respectively, with rotating speed Ω = 60 rpm.
(a)
U
=
50 m s−1
(b)
U
=
55 m s−1


### 4.4. Radial Basic Function Neural Network PID Control Strategy

In view of the variability of vibration amplitudes and steady-state values in different simulation processes of above-mentioned SN PID controller, an adaptive radial basic function neural network (RBFNN) is developed and employed to perform PID control. RBFNN has the radial basis function neurons for the hidden layer and linear neurons for the output layer [21–23]. Firstly, the nonlinear map implementation from input x to h
j(x) is performed in the hidden layer. In this process the radial basis function in the hidden layer gives a local response to the input signals. Commonly, radial basis functions are chosen as Gaussian function, expressed as(25)
h
j
x
=
exp
⁡
-
x
-
c
j
2
2
b
j
2
,
j
=
1,2
,
…
,
m
,where x is the dimensional input vector with elements x
j and c
j is the vector determining the centre of the basis function h
j. Then the linear map implementation from h
j(x) to y
m(k) is performed in the output layer:(26)
y
m
k
=
ω
1
h
1
+
ω
2
h
2
+
⋯
+
ω
m
h
m
.The performance index function of identifier is expressed as(27)
J
=
1
2
y
k
-
y
m
k
2
.Based on the gradient descent method, iterative algorithms of output weights, basis width parameters, and the centre functions are expressed as(28)
Δ
ω
j
k
=
η
y
k
-
y
m
k
h
j
,
ω
j
k
=
ω
j
(
k
-
1
)
+
Δ
ω
j
(
k
)
+
α
·
ω
j
k
-
1
-
ω
j
k
-
2
,
Δ
b
j
k
=
η
y
k
-
y
m
k
ω
j
h
j
X
-
C
2
b
j
3
,
b
j
k
=
b
j
(
k
-
1
)
+
Δ
b
j
(
k
)
+
α
·
b
j
k
-
1
-
b
j
k
-
2
,
Δ
c
j
i
k
=
η
y
k
-
y
m
k
ω
j
x
j
-
c
j
i
b
j
2
,
c
j
i
k
=
c
j
i
(
k
-
1
)
+
Δ
c
j
i
(
k
)
+
α
·
c
j
i
k
-
1
-
c
j
i
k
-
2
,where η is learning rate and α is momentum factor.Summarizing the above algorithm, PID input itemsx
c, control algorithm, and neural network tuning indicator are, respectively, expressed as(29)
x
c
1
=
e
k
-
e
k
-
1
,
x
c
(
2
)
=
e
(
k
)
x
c
(
3
)
=
e
(
k
)
-
2
e
(
k
-
1
)
+
e
(
k
-
2
)
u
k
=
u
k
-
1
+
Δ
u
k
,
Δ
u
k
=
K
P
e
k
-
e
k
-
1
+
K
I
e
k
+
K
D
e
k
-
2
e
k
-
1
+
e
k
-
2
E
k
=
1
2
e
k
2
,where the tuning of PID parameters is based on the gradient descent method, detailed as(30)
Δ
K
P
=
-
η
∂
E
∂
K
P
=
-
η
∂
E
∂
y
∂
y
∂
Δ
u
∂
Δ
u
∂
K
P
=
η
e
k
∂
y
∂
Δ
u
x
c
1
,
Δ
K
I
=
-
η
∂
E
∂
K
I
=
-
η
∂
E
∂
y
∂
y
∂
Δ
u
∂
Δ
u
∂
K
I
=
η
e
k
∂
y
∂
Δ
u
x
c
2
,
Δ
K
D
=
-
η
∂
E
∂
K
D
=
-
η
∂
E
∂
y
∂
y
∂
Δ
u
∂
Δ
u
∂
K
D
=
η
e
k
∂
y
∂
Δ
u
x
c
3
.Herein, ∂
y
/
∂
Δ
u is Jacobi information of the controlled object, obtained by the neural network identification.Figures15(a)-15(b) show the flap and lag displacements and corresponding pitch angles of RBFNN PID control for U
=
50 m s−1 and U
=
55 m s−1, respectively. Here the rotating speed Ω = 60 rpm. By comparison with the simulation results of SN PID control in Figure 14, RBFNN PID has obvious advantages from the point of view of vibration amplitudes and steady-state values. Furthermore, different from SN PID control, the simulation results of RBFNN PID are almost invariant in different simulation processes given that the relevant parameters are invariant, which might lie in the local approximation performance of RBFNN. Figures 15(c)-15(d) show the flap and lag displacements and corresponding pitch angles of RBFNN PID control for U
=
60 m s−1 and U
=
80 m s−1, respectively. The orders of magnitude of vibration amplitudes are consistent with the increasing of the wind speeds. Even under the condition of high wind speed U
=
150 m s−1 (theoretical value) in Figure 15(e), RBFNN PID control can achieve good effects. However, under the condition of higher theoretical speed of U
=
180 m s−1 in Figure 15(f), the efficacy of RBFNN PID control is lost with displacements being divergent. Simultaneously the proper pitch angle is not adjustable and available. In addition the disadvantage of RBFNN PID control is that the construction process of system hardware link is more complicated.Figure 15
Flap and lag displacements and corresponding pitch angles of RBFNN PID control forU = 50, 55, 60, 80, 150, and 180 m s−1, respectively.
(a)
U
=
50 m s−1
(b)
U
=
55 m s−1
(c)
U
=
60 m s−1
(d)
U
=
80 m s−1
(e)
U
=
150 m s−1
(f)
U
=
180 m s−1


## 4.1. The Linearized Aeroservoelastic Model and PID Control Strategy

Inserting (17) into (12), the linearized aeroservoelastic model can be rewritten as(18)
M
X
¨
1
+
C
X
˙
1
+
K
X
1
=
0
,where the M, C, and K are 6 × 7 coefficient matrices, and furthermore, (18) can be transformed as(19)
X
¨
1
=
A
M
2
X
˙
1
+
A
M
1
X
1
+
W
P
K
P
X
˙
1
+
K
I
∫
X
˙
1
d
t
+
K
D
d
X
˙
1
d
t
,where the coefficient matrices A
M
2, A
M
1, and W
P are displayed in Appendix B. The block diagram of (19) in Simulink is shown in Figure 4.Figure 4
Block diagram of PID control strategy in Simulink.The PID parameters,K
P, K
I, and K
D are optimized by optimal control. The optimal control is under certain concrete conditions, to achieve special control object and to make selected target maximal or least. From the optimal control, the PID parameters can be found during optimization. The target function used here adopts ITAE rule, defined as(20)
J
ITAE
=
∫
0
∝
t
e
(
t
)
d
t
.When time t is relatively big, to assure small target value, the steady error must be diminished, so as to make system speed in steady area. Each parameter in target transfer function may be evaluated and initialized.The Matlab functions are applied to search for the optimal parameters for the system [18]. It can be seen from the practice that the program is quite versatile in finding the optimal controllers. However, in some disturbed states, it may not find a solution due to unstable controlled plant. This can be a drawback in conventional optimization algorithms, but such problems can be avoided by multiple choices of initial values based on genetic algorithm [19].In order to improve the computational speed, the optional PI controller (K
D
=
0) or PD controller (K
I
=
0) is applied. Throughout the numerical illustrations, the aeroelastic instability case of U
=
30 m s−1 and Ω
=
60 rpm is considered with aeroservoelastic control implemented in Figure 5.Figure 5
Aeroelastic instability and aeroservoelastic control case ofU
=
30 m s−1 and Ω = 60 rpm.
(a)
Flap and lag displacements for the case ofU
=
30 m s−1 and Ω = 60 rpm
(b)
Lag displacements for the three optimal aeroservoelastic controllers
(c)
Flap displacements for the three optimal aeroservoelastic controllersFigure5(a) shows that both flap and lag motions are divergent (instable): flapwise displacement is positively divergent; the lag one is negatively divergent. Figures 5(b)-5(c) depict the aeroservoelastic control results of stall-induced flap/lag flutter for the three optimal aeroservoelastic controllers mentioned above.It demonstrates the obviously improved aeroelastic stability with rapidly convergent tendency for both flap and lag motions. With the change of time, the flutter amplitude decreases rapidly and gradually tends to be zero. The effect of flutter suppression is obvious. In contrast to flap case, effect of aeroservoelastic control on lag case is more prominent, and furthermore PID and PD controllers are better means of aeroservoelastic control. As for aeroservoelastic control in flapwise direction, PI controller is a good choice.


## 4.2. Fuzzy PID Control Strategy

Due to the fact that the controller in Figure4 is based on the full state feedback including the aerodynamic variables, the realization of this full state feedback controller is difficult in practice. So the simulation in Figure 4 is only the theoretical research and an object of reference in order to make comparison with other physically achievable feedback controllers for structural or mechanical variables such as fuzzy controllers and neural network controllers below.A self-tuning fuzzy PI (FPI) PID controller is designed, which can automatically tune the parameters (K
P and K
I) of PI controller. The advantage of the fuzzy controller here is that it is physically achievable feedback controllers for mechanical variable ϕ
˙
1, compared with optimal controllers. During real-time control, the self-tuning process is carried out by PLC instead of optimal control, which can also avoid large amount of calculation of genetic algorithm in the process of optimization in some disturbed states.Firstly rewrite (19) as(21)
X
¨
1
=
A
M
2
X
˙
1
+
A
M
1
X
1
+
W
P
1
K
P
ϕ
˙
1
+
K
I
∫
ϕ
˙
1
d
t
+
K
D
d
ϕ
˙
1
d
t
,where W
P
1
=
M
-
1
[
0
0
0
0
0
1
]
T. The block diagram of FPI PID control in Simulink is illustrated in Figure 6.Figure 6
The block diagram of FPI PID strategy.Fuzzy control depends on the fuzzy algorithm between the information of process and control input. Fuzzy controllers from their inception have demonstrated a vast range of applicability to processes where the plant transfer function is not defined but the control action can be described in terms of linguistic variables and are often being used to improve the performance of a system where the plant transfer function is known [20]. For this physically achievable feedback system, FPI controller is devised. Two input variables, error (E) and change in error (E
D), and two output variables K
P and K
I with five linguistic variables of membership function are used. The input linguistic variables aremf1 (negative medium),mf2 (negative small),mf3 (zero),mf4 (positive small), andmf5 (positive medium). Five triangular-membership function forms for E and five triangular-membership function forms forED(de/dt) are determined the same way, which are shown in Figure 7. Borders of both function sequences vary between +3 and −3. In this study, a fuzzy rule type of the designed FPI controller is developed based on Sugeno model in Simulink. So there are 25 weight values in both K
P and K
I designs. According to intuition method, list of linguistic rules is shown in Tables 1-2.Table 1
Rules forK
p.
E∖ED
m
f
1
m
f
2
m
f
3
m
f
4
m
f
5
m
f
1
1
0.85
0.65
0.5
0.5
m
f
2
0.85
0.85
0.65
0.5
0.5
m
f
3
0.85
0.65
0.5
0.35
0.15
m
f
4
0.65
0.5
0.35
0.35
0.15
m
f
5
0.65
0.5
0.35
0.15
0.15Table 2
Rules forK
I.
E∖ED
m
f
1
m
f
2
m
f
3
m
f
4
m
f
5
m
f
1
0
0.5
1
1
1.5
m
f
2
0.5
1
1
1.5
2
m
f
3
0.5
1
1.5
2
2.5
m
f
4
1
1.5
2
2
2.5
m
f
5
1.5
2
2
2.5
3Figure 7
Membership function plots for both function sequences.In addition, a self-tuning fuzzy PD (FPD) PID controller is designed. Its main advantage is that FPD controller parameters can be changed very quickly in response to changes of the system dynamics. Five triangular-membership function forms forE and E
D are determined the same way as Figure 7. The block diagram of FPD PID control strategy is illustrated in Figure 8. According to intuition method, list of linguistic fuzzy rules is shown in Table 3. Simulation is performed by substituting structure in Figure 6 for corresponding one of Figure 8.Table 3
Rule weight values of FPD PID controller.
E∖ED
m
f
1
m
f
2
m
f
3
m
f
4
m
f
5
m
f
1
−16
−12.8
−2.4
−8
0
m
f
2
−13.6
−12
−6.4
0
11.2
m
f
3
−14.4
−8
0
10.4
16
m
f
4
−9.6
0
13.6
12.8
15.2
m
f
5
0
8
12
15.2
16Figure 8
The block diagram of FPD PID strategy.In order to make qualitative analysis for the effectiveness of the fuzzy PID (including FPI and FPD) control algorithm, some fixed speeds in different ranges of wind speeds are selected.Case One. Consider U
=
5 m s−1 and U
=
10 m s−1, respectively, in low range of wind speeds for less than 17 m s−1.Case Two. Consider U
=
18 m s−1, an approximate intermediate value of wind velocity for the choice of the two types of fuzzy PID controllers.Case Three. Consider U
=
20 m s−1 and U
=
24 m s−1, respectively, in high range of wind speeds for more than 19 m s−1.Figures9(a)-9(b) depict the linearized responses of flap and lag displacements of Case one for different fuzzy controllers, which shows obvious superiority in flapwise stability of FPD PID controller; Figure 10 depicts the linearized responses of flap and lag displacements of Case two for the intermediate wind velocity of U
=
18 m s−1, which shows almost equal stability of both FPI and FPD PID controllers; Figures 11(a)-11(b) depict the linearized responses of flap and lag displacements of Case three for different fuzzy controllers, which shows obvious superiority in stability of FPI PID controller. It can be seen in the process of qualitative analysis that in low range of wind speeds FPD PID controller has the superiority compared with FPI PID controller. However simulation result indicates that the FPI PID controller is more effective and predominant than FPD controller in high range of wind speeds. Here the rotating speed Ω = 60 rpm.Figure 9
Flap and lag displacements of Case one for different fuzzy controller.
(a)
Displacements forU
=
5 m s−1
(b)
Displacements forU
=
10 m s−1Figure 10
Flap and lag displacements of Case two forU
=
18 m s−1.Figure 11
Flap and lag displacements of Case three for different fuzzy controllers.
(a)
Displacements forU
=
20 m s−1
(b)
Displacements forU
=
24 m s−1It should be stated that the intermediate value of wind velocity for the choice of the two types of fuzzy PID controller varies with the rotating speed. For example, the intermediate value of wind velocity isU
=
24.5 m s−1 for rotating speed Ω = 120 rpm.


## 4.3. Single Neuron PID Control Strategy with Improved Hebb Learning Algorithm

In general, under the condition of normal wind speed in the range ofU
<
32 m s−1, as mentioned above, the FPI PID control can completely meet the requirements. However in some special cases of extreme wind speeds, such as more than 50 m s−1, FPI algorithm might fail because of highly complicated and nonlinear instability. A single neuron (SN) PID control with improved Hebb learning algorithm can make up for the defect from the aspects of numerical simulation (due to the lack of real-time operation environment in the range of U
>
50 m s−1). With the interdisciplinary development of modern computational technologies and biological science, artificial neural networks, as typical artificial intelligence models and algorithms, have become an attractive approach for modeling highly complicated and nonlinear system without understanding the nature of phenomenon [21].Using the improved Hebb learning algorithm, update rules of the three weights can be written as [22, 23](22)
ω
1
(
k
)
=
ω
1
(
k
-
1
)
+
η
P
e
k
u
k
e
k
-
Δ
e
k
ω
2
(
k
)
=
ω
2
(
k
-
1
)
+
η
I
e
(
k
)
u
(
k
)
[
e
(
k
)
-
Δ
e
(
k
)
]
ω
3
k
=
ω
3
k
-
1
+
η
D
e
k
u
k
e
k
-
Δ
e
k
,where η
P, η
I, and η
D are adaptive proportional, integral, and differential learning rates, respectively. The three weights are the system state variables, so the control law is written as(23)
u
k
=
u
k
-
1
+
K
∑
i
=
1
3
ω
i
0
k
x
i
k
,where the normalized weights and the three parameters are(24)
ω
i
0
k
=
ω
i
k
∑
i
=
1
3
ω
3
k
,
x
1
k
=
e
k
,
x
2
k
=
Δ
e
k
=
e
k
-
e
k
-
1
,
x
3
k
=
Δ
2
e
k
=
e
k
-
2
e
k
-
1
+
e
k
-
2
.Summarizing the above algorithm, the SN PID control block diagram with algorithmS-function in Matlab and the system simulation diagram are illustrated in Figure 12.Figure 12
SN PID control block diagram and the system simulation diagram.
(a)
SN PID controller
(b)
System simulation diagramFigures13(a)–13(c) show three simulation results of flap and lag displacements of SN PID controller for U
=
30 m s−1 and Ω = 60 rpm. Figure 13(d) shows flap and lag displacements of FPI PID controller. It can be seen from Figures 13(a)–13(c) that the three simulation results including vibration amplitudes and steady-state values are quite different although system stability is the same and convergent. Among them Figure 13(c) is the most effective from the magnitude of vibration amplitudes of view. Figures 13(c)-13(d) show the difference between the third simulation of SN PID controller and simulation of FPI PID controller. It is seen that the flutter amplitudes of FPI PID controller in both flapwise and edgewise directions are smaller than those of SN PID controller.Figure 13
Flap and lag displacements of SN PID control (a)–(c) and FPI PID control (d) forU
=
30 m s−1 and Ω
=
60 rpm.
(a)
The first simulation of SN PID
(b)
The second simulation of SN PID
(c)
The third simulation of SN PID
(d)
Simulation of FPI PIDThe advantages of SN PID controller over FPI PID controller lie in its performance in the range of extreme wind speeds, such as more than 50 m s−1. Figures 14(a)-14(b) show the flap and lag displacements of SN PID control for U
=
50 m s−1 and U
=
55 m s−1, respectively. Although there might be some deviations from zero in flapwise motions, the control effects of SN PID controller are satisfactory within the range of extreme wind speeds where the simulation algorithm of FPI PID controller is not feasible because it is unable to reduce the step size without violating the minimum step size for a certain number of consecutive times during simulation. This might be caused by dependent dynamic states or by the high stiffness of the system.Figure 14
Flap and lag displacements of SN PID control forU
=
50 m s−1 and U
=
55 m s−1, respectively, with rotating speed Ω = 60 rpm.
(a)
U
=
50 m s−1
(b)
U
=
55 m s−1


## 4.4. Radial Basic Function Neural Network PID Control Strategy

In view of the variability of vibration amplitudes and steady-state values in different simulation processes of above-mentioned SN PID controller, an adaptive radial basic function neural network (RBFNN) is developed and employed to perform PID control. RBFNN has the radial basis function neurons for the hidden layer and linear neurons for the output layer [21–23]. Firstly, the nonlinear map implementation from input x to h
j(x) is performed in the hidden layer. In this process the radial basis function in the hidden layer gives a local response to the input signals. Commonly, radial basis functions are chosen as Gaussian function, expressed as(25)
h
j
x
=
exp
⁡
-
x
-
c
j
2
2
b
j
2
,
j
=
1,2
,
…
,
m
,where x is the dimensional input vector with elements x
j and c
j is the vector determining the centre of the basis function h
j. Then the linear map implementation from h
j(x) to y
m(k) is performed in the output layer:(26)
y
m
k
=
ω
1
h
1
+
ω
2
h
2
+
⋯
+
ω
m
h
m
.The performance index function of identifier is expressed as(27)
J
=
1
2
y
k
-
y
m
k
2
.Based on the gradient descent method, iterative algorithms of output weights, basis width parameters, and the centre functions are expressed as(28)
Δ
ω
j
k
=
η
y
k
-
y
m
k
h
j
,
ω
j
k
=
ω
j
(
k
-
1
)
+
Δ
ω
j
(
k
)
+
α
·
ω
j
k
-
1
-
ω
j
k
-
2
,
Δ
b
j
k
=
η
y
k
-
y
m
k
ω
j
h
j
X
-
C
2
b
j
3
,
b
j
k
=
b
j
(
k
-
1
)
+
Δ
b
j
(
k
)
+
α
·
b
j
k
-
1
-
b
j
k
-
2
,
Δ
c
j
i
k
=
η
y
k
-
y
m
k
ω
j
x
j
-
c
j
i
b
j
2
,
c
j
i
k
=
c
j
i
(
k
-
1
)
+
Δ
c
j
i
(
k
)
+
α
·
c
j
i
k
-
1
-
c
j
i
k
-
2
,where η is learning rate and α is momentum factor.Summarizing the above algorithm, PID input itemsx
c, control algorithm, and neural network tuning indicator are, respectively, expressed as(29)
x
c
1
=
e
k
-
e
k
-
1
,
x
c
(
2
)
=
e
(
k
)
x
c
(
3
)
=
e
(
k
)
-
2
e
(
k
-
1
)
+
e
(
k
-
2
)
u
k
=
u
k
-
1
+
Δ
u
k
,
Δ
u
k
=
K
P
e
k
-
e
k
-
1
+
K
I
e
k
+
K
D
e
k
-
2
e
k
-
1
+
e
k
-
2
E
k
=
1
2
e
k
2
,where the tuning of PID parameters is based on the gradient descent method, detailed as(30)
Δ
K
P
=
-
η
∂
E
∂
K
P
=
-
η
∂
E
∂
y
∂
y
∂
Δ
u
∂
Δ
u
∂
K
P
=
η
e
k
∂
y
∂
Δ
u
x
c
1
,
Δ
K
I
=
-
η
∂
E
∂
K
I
=
-
η
∂
E
∂
y
∂
y
∂
Δ
u
∂
Δ
u
∂
K
I
=
η
e
k
∂
y
∂
Δ
u
x
c
2
,
Δ
K
D
=
-
η
∂
E
∂
K
D
=
-
η
∂
E
∂
y
∂
y
∂
Δ
u
∂
Δ
u
∂
K
D
=
η
e
k
∂
y
∂
Δ
u
x
c
3
.Herein, ∂
y
/
∂
Δ
u is Jacobi information of the controlled object, obtained by the neural network identification.Figures15(a)-15(b) show the flap and lag displacements and corresponding pitch angles of RBFNN PID control for U
=
50 m s−1 and U
=
55 m s−1, respectively. Here the rotating speed Ω = 60 rpm. By comparison with the simulation results of SN PID control in Figure 14, RBFNN PID has obvious advantages from the point of view of vibration amplitudes and steady-state values. Furthermore, different from SN PID control, the simulation results of RBFNN PID are almost invariant in different simulation processes given that the relevant parameters are invariant, which might lie in the local approximation performance of RBFNN. Figures 15(c)-15(d) show the flap and lag displacements and corresponding pitch angles of RBFNN PID control for U
=
60 m s−1 and U
=
80 m s−1, respectively. The orders of magnitude of vibration amplitudes are consistent with the increasing of the wind speeds. Even under the condition of high wind speed U
=
150 m s−1 (theoretical value) in Figure 15(e), RBFNN PID control can achieve good effects. However, under the condition of higher theoretical speed of U
=
180 m s−1 in Figure 15(f), the efficacy of RBFNN PID control is lost with displacements being divergent. Simultaneously the proper pitch angle is not adjustable and available. In addition the disadvantage of RBFNN PID control is that the construction process of system hardware link is more complicated.Figure 15
Flap and lag displacements and corresponding pitch angles of RBFNN PID control forU = 50, 55, 60, 80, 150, and 180 m s−1, respectively.
(a)
U
=
50 m s−1
(b)
U
=
55 m s−1
(c)
U
=
60 m s−1
(d)
U
=
80 m s−1
(e)
U
=
150 m s−1
(f)
U
=
180 m s−1


## 5. Conclusions

In this study, aeroelastic stability of stall-induced flutter of special blade section with flap and lag motions based on unsteady dynamic stall nonlinear aerodynamic model is investigated. Moreover aeroservoelastic control for aeroelastic instability is implemented by different controllers. Some concluding remarks can be drawn from the results.(
1
) The influences of stall-induced flap/lag flutter for different parameters, wind speed and rotating speed, are investigated based on nonlinear aerodynamic model. Stability analysis methods and solution methodology are described in this study. The nonlinear resulting equations are then linearized for small perturbation about the equilibrium point. The stability characteristics and linearized aeroelastic responses are investigated by secant method and Runge-Kutta time-marching approach.(
2
) Pitch control is rarely investigated in the field of stall-induced flutter suppression for wind turbine blade. Aeroservoelastic pitch control for aeroelastic instability is implemented by PID controller with control parameters tuned by optimal control theory and other adaptive processes. For the normal range of wind speeds, FPI and FPD controllers based on fuzzy control theory are adopted and good results are achieved. Furthermore, in the range of extreme wind speeds, SN PID control strategy with improved Hebb learning algorithm and RBFNN PID control are applied. Generally, both of them are of no application limitations.(
3
) It is obviously demonstrated that the method of aeroservoelastic control is robust and the effects of flutter suppression are apparent. It should be stated that the main purpose of the control strategy mentioned here is to gain flutter suppression approaches about the interaction between rotating blade section and imposed control action. Flutter suppression for the whole blade body under coupling conditions is not involved here.



---

*Source: 692567-2015-03-18.xml*