# Improved Grey Particle Swarm Optimization and New Luus-Jaakola Hybrid Algorithm Optimized IMC-PID Controller for Diverse Wing Vibration Systems

Authors: Nailu Li; Hua Yang; Anle Mu
Journal: Complexity
 (2019)

---

## Abstract

The PID control plays important role in wing vibration control systems. However, how to efficiently optimize the PID parameters for different kinds of wing vibration systems is still an open issue for control designers. The problem of PID control optimization is first converted into internal mode control based PID (IMC-PID) parameters optimization problem for complex wing vibration systems. To solve this problem, a novel optimization technique, called GNPSO is proposed based on the hybridization of improved grey particle swarm optimization (GPSO) and new Luus-Jaakola algorithm (NLJ). The original GPSO is modified by using small population size/iteration number, employing new grey analysis rule and designing new updating formula of acceleration coefficients. The hybrid GNPSO benefits improved global exploration of GPSO and strong local search of new Luus-Jaakola (NLJ), so as to avoid arbitrary and inefficient search of global optimum and prevent the trap in local optimum. Diverse wing vibration systems, including linear system, nonlinear system and multiple-input-multiple-output system are considered to verify the effectiveness of proposed method. Simulation results show that GNPSO optimized method obtains improved vibration control performance, stronger robustness and wide applicability on all system cases, compared to existing evolutionary algorithm based tuning methods. Enhanced optimization convergence and computation efficiency obtained by GNPSO tuning technique are also verified by statistical analysis.


---

## Body

## 1. Introduction

Control of wing vibration systems is a demanding task. Small damped vibrations, causing over fatigue load and threaten the safe operation of the wing, require vibration suppression in multiple degrees of freedoms (DOF) by the control surface installed on the wing [1, 2]. Besides, the system nonlinearity existing in aerodynamic and structural forces brings additional challenge in effective control of wing vibration systems [3, 4]. Moreover, the complex nonlinear wing vibration system [5], involving multiple control surfaces needs the accurate controller design to properly manipulate all surfaces, overcome the effect of system nonlinearity and achieve the vibration control in all DOFs.In the wing vibration control system, the proportional-integral-derivative (PID) controller [6, 8] has been developed and most widely used to suppress the vibrations for linear wing vibration system. For nonlinear wing vibration systems, sliding mode control [9, 10], robust control [11, 12] and H
∞ control [13] have been studied to control aeroelastic vibrations under unsteady flow and nonlinear spring. To deal with multiple-input-multiple-output (MIMO) wing vibration system, a state estimation based adaptive output feedback controller [14] was designed for vibration suppression on a nonlinear wing section. Also, a continuous robust controller [15] was proposed for control of a nonlinear wing section with a leading-and trailing control surface.Although those advanced controllers have improved control performance for complex nonlinear wing vibration systems, they have quite limitations of implementation due to the high order and complicated control structure. Nevertheless, the PID controller has the strength in implementation due to its simple structure and promising application for a wide range of different wing vibration systems. However, it is always time-consuming and difficult to obtain the proper coefficients of PID, especially for complex wing vibration control system, which could involve system nonlinearity or large number of control parameters for MIMO case.To solve the problem of control parameter design, evolutionary optimization techniques are considered. Based on the defined objection function, the problem of control parameter design can be converted into a problem of the parameter optimization. Besides, in order to efficiently determine the parameters of PID controller, internal mode control (IMC) based PID structure has been applied to obtain PID parameters for diverse industrial systems [16–19]. Several evolutionary algorithms have been applied for IMC-PID parameters optimization. Multiple objective based genetic algorithm (GA) [20] was adopted for tuning the parameter of IMC-PID controller for linear system. Particle swarm intelligence (PSO) [21] based technique was applied to obtain better IMC-PID control parameter than conventional Ziegler Nichols (Z-N) tuning method for first order plus delay temperature process. Beside, hybrid of differential evolution and new Luus-Jaakola algorithm (DE–NLJ) [22] was proposed as a new IMC-PID tuning method and achieved improved control parameters than classic PID tuning methods. An NPSO hybrid algorithm [23] based method obtained improved control performance than other PID tuning methods for a MIMO low-order system. However, those algorithm based tuning methods are mainly targeted to linear system, low-order systems and first-order plus delay process. It could be much difficult to search for global optimum when it comes to complex high-order nonlinear wing vibration systems. On this condition, the global search could be arbitrary and inefficient, and the optimization search could easily be trapped into local optimum under the complexity of wing vibration systems. The vibration control performance can be degraded with unsatisfied optimized results. Besides, long optimization time due to large population size P and large iterations G are not preferred by practical control systems. Thus, the effective and low-cost optimization algorithm, which can handle control optimization problem for complex wing vibration systems, needs to be further investigated.In recent years, a highly efficient grey-based PSO (GPSO) algorithm [24] was designed to optimize the PID parameters for a linear second-order system. Improved control performance, reduced population size and smaller iteration number were achieved, compared to GA and PSO algorithm. However, it is indicated that the performance of the evolutionary algorithm is problem-dependent [25, 26]. Since nonlinear high-order wing vibration system has quite different dynamic characteristics than the linear low-order system in reference [24], more suitable optimization algorithm should be studied to solve the control optimization problem in this paper, so as to obtain ensure the control performance. Moreover, all previous studies focused on the application for a specific type of system. An applicable optimization technique for a wide range of PID control systems, including linear system, nonlinear system, and MIMO system, has rarely been studied.Based on the above discussion, there is a motivation of developing a new hybrid evolutionary algorithm, call GNPSO to tune PID parameters for expanded applications in diverse wing vibration systems, from simple linear case to nonlinear MIMO case. In this paper, the problem of PID control optimization design for different wing vibration system is first converted into the IMC-PID parameter optimization problem under control constraints. A hybrid algorithm of improved GPSO and NLJ (GNPSO) with small population size/iteration number is proposed based on new grey analysis rule and new updating formula of algorithm parameters. Then, the GNPSO algorithm is applied to optimize the IMC-PID parameter in closed-loop wing vibration control system. The fitness function is defined as a linear function regarding vibration control performance in the time domain. The effectiveness of the proposed method is examined on a wide range of different wing vibration systems. The superiority of proposed GNPSO algorithm to classic algorithms and other hybrid algorithms in previous studies is studied by optimization performance, statistical results, vibration control performance and computation cost.The contributions of this paper are presented as follows:(1)
The task of PID control design for different type of wing vibration system is converted into the IMC-PID parameter optimization problem, and the GNPSO algorithm is presented to automatically optimize the control parameters according to different system case.
(2)
In proposed hybrid GNPSO, the improved GPSO (IGPSO) is designed by introducing new grey analysis rule to avoid arbitrary and inefficient global search, and new updating formula of acceleration coefficients to fasten convergence. The local search is enhanced by NLJ algorithm to reach the global optimum. Small population size and small iteration number are obtained for low-cost and efficient tuning of PID controller.
(3)
Proposed method expands the application of PID control to a wide range of wing vibration systems, including linear system, nonlinear system and MIMO systems, while feasible PID control was mostly limited to linear wing vibration system.
(4)
The superiority of GNPSO algorithm to existing tuning algorithms (i.e., classic algorithms [20, 21], hybrid algorithms [22–24]) and recent grey-based algorithm [27] is comprehensively verified based on convergence performance, different statistical tests, vibration control performance and robustness tests for parameters optimization.
(5)
Reduced computation cost of GNPSO based tuning method is also confirmed by less computation time on all system cases, compared to the existing evolutionary algorithm based tuning methods.The rest of the paper is arranged as follows. Section2 presents the IMC-PID controller for generalized wing vibration system and the formulation of control optimization problem. Section 3 presents the design of a hybrid GNPSO algorithm, its improvements and optimization procedures. Section 4 illustrates the simulation tests, statistical tests and robustness tests of proposed GNPSO on linear, nonlinear and MIMO wing vibration control systems separately, compared to classic evolutionary algorithms and other published hybrid algorithms. Computation efficiency of GNPSO is also discussed in all cases. Section 5 gives the conclusion of the paper.


## 2. Wing Vibration Control

### 2.1. System Description and Control Specifications

The wing vibration system is normally a rigid airfoil with pitch and plunge motions, which are manipulated by the flap control surface. The pitch and plunge motions are defined by pitch angleθ and plunge displacement h separately. Assuming the variable vector q
=
h
θ
T and the control signal of flap angle u
, general dynamic equations of the wing vibration system can be presented as:(1)
M
q
q
¨
+
C
q
q
˙
+
K
q
q
=
F
q
,
q
˙
,
u
,where the left side of the equation represents the structural dynamics of the wing vibration system with the mass matrixM, the damping matrix C and the spring matrix K. The right side of the equation describes the aerodynamics with respect to the vibration motions and the control command. As indicated in references [8, 12, 15], there are mainly three different types of the wing vibration system based on Equation (1). The first is linear wing vibration system with linear spring/damping characteristics and single control surface. The second is nonlinear wing vibration system with nonlinear spring/damping/aerodynamics and single control surface. The last is nonlinear MIMO wing vibration system, which has two control surfaces of both leading edge flap and trailing edge flap. The system nonlinearity and multiple control surfaces could further increase the difficulty in design of control parameters for the wing vibration system.The control specifications for wing vibration systems are:(1)
The wing vibration in pitch motion should be well suppressed.
(2)
The wing vibration in plunge motion should be well suppressed.
(3)
The satisfactory control performance should be obtained under control constraints.
(4)
The vibration control performance must be ensured under system nonlinearity.
(5)
The tear-and-wear of actuators should be minimized.
(6)
A sufficient closed-loop stability should be shown.In order to meet all the specific control requirements, the objective functions are usually given based on the vibration control performance, indicated by Integral of Absolute Error (IAE), Integral of Squared Error (ISE), Integral of Time multiplied Absolute Error (ITAE) and Integral of Time multiplied Square Error (ITSE). In this paper, IAE based objective function is taken as it emphasizes the instant response and the damping behavior, which are appropriate for the evaluation of the control performance for wing vibration systems.


### 2.2. Controller Structure and Problem Formulation

#### 2.2.1. Internal Model Control (IMC) Based PID Structure

In wing vibration control system, PID controllers are mostly used to suppress the vibrations in plunge and pitch motions. However, the shortcomings of using the conventional PID controller include slow closed-loop responses, unsatisfied robustness to system change, inefficient utilization of the control surface and great control design effort. The reason is the difficulty in design of proper PID parameters, especially when complex system is involved, such as the nonlinear system and the MIMO system. An internal model control (IMC) based PID controller improves the settling time, overshoot in the transient response and the robustness of the classic PID controller, since the IMC can handle the model mismatch due to the nonlinearity or system change by the filter design and the process model [29, 30]. The filter in IMC-PID is designed to ensure the dynamic performance and improve the robustness of the PID controller under varied system parameters [31]. Meanwhile, three PID parameters can be automatically computed and updated by adjusting the filter parameter in the IMC-PID structure.To design the IMC controller for wing vibration systems, the process modelG
m
s of the system should first be obtained. In general, the process model of wing vibration systems can be approximated by the reduced order model G
R
O
M
sas:(2)
G
m
s
=
G
R
O
M
s
=
k
0
−
β
s
+
1
a
1
s
2
+
a
0
s
+
1
,wherek
0 is the system gain, β
,
a
1
,
a
0 are system parameters of the reduced order model. a
1
>
0
,
a
0
>
0. In this paper, the least-square method is adopted to identify the reduced order model of wing vibration systems.The control schematic of the IMC is shown in Figure1(a), where the feedback IMC controller is composed of the process model G
m
s and the controller C
s. The controller C
s is designed based on the filter f
s as [32]:(3)
C
s
=
f
s
G
m
−
s
,
G
m
s
=
G
m
−
s
G
m
+
s
,whereG
m
−
s
=
k
0
/
a
1
s
2
+
a
0
s
+
1is the minimum phase part and G
m
+
s
=
−
β
s
+
1 is the nonminimum phase part of the process model. The low-order filter f
s has the filter parameter λ, which can be adjusted to improve the dynamic control performance under system nonlinearity or system changes. In this paper, second-order filter is used as f
s
=
1
/
λ
s
+
1
2
.As indicated in Figure1, the feedback IMC controller G
c
s is presented as:(4)
G
c
s
=
C
s
1
−
G
m
s
C
s
.Substituting Equations (2) and (3) into Equation (4), the feedback IMC controller can be rewritten as:(5)
G
c
s
=
a
1
s
2
+
a
0
s
+
1
k
0
λ
2
s
2
+
2
λ
+
β
λ
s
.The equivalent IMC-PID controllerG
I
M
C
-
P
I
D
s can be derived from G
c
s as:(6)
G
I
M
C
-
P
I
D
s
=
G
c
s
≈
k
p
+
k
i
s
+
k
d
s
,
k
p
λ
=
−
1
2
k
0
λ
+
k
0
β
a
0
−
λ
2
2
λ
+
β
,
k
i
λ
=
−
1
2
k
0
λ
+
k
0
β
,
k
d
λ
=
−
1
2
k
0
λ
+
k
0
β
a
1
−
λ
2
2
λ
+
β
a
0
−
λ
2
2
λ
+
β
.For the design of IMC-PID controller for different kinds of wing vibration systems, the control parameters of the IMC-PID are obtained as followings.(1)
Linear wing vibration system: the process model of the original fourth order system is taken as a second-order model in Equation (2) by model identification. Then, the control parameters of the IMC-PID are derived based on Equations (3–6). The low-order process model is favored to avoid the high order of the controller.
(2)
Nonlinear wing vibration system: the process model of the nonlinear system is approximated by the linear reduced order model at operating points. The process model in Equation (2) is identified by the least-square method. Then, the IMC-PID control parameters are calculated according to Equations (3)–(6).
(3)
MIMO wing vibration system:The MIMO system usually has two system outputs and two control inputs. For MIMO system, two SISO control loops are formed, where one IMC-PID controller is designed independently for each control loop. For the controller design in each SISO control loop, the process model of the corresponding SISO system is identified as the reduced order model in Equation (2). Then, the parameters of the identified process model are used to calculate the parameters of the IMC-PID controller using Equations (3)–(6).It should be noted that the identified process model could vary according to different wing vibration systems. On this condition, the parameters of the corresponding IMC-PID controller could vary from case to case.As mentioned above, the implemented structure of IMC- PID based wing vibration control system is shown in Figure1(b). The vibration deflection, such as pitch angle is chosen as the system output y. The control input u is defined as the flap angle. The identified process model G
m
s and the designed filter f
λ are used to formulate the feedback IMC vibration controller G
c
s, which is further implemented in the form of PID controller. The control constraints, such as the saturation and rate limit are also considered in the control system. The filter is optimally designed by proposed GNPSO algorithm according to defined objective function J. The optimal design problem of the filter is presented in next section. The design of proposed GNPSO algorithm is detailed in Section 3.Figure 1
(a) Control schematic of Internal model control (IMC) system, (b) implemented structure of IMC-PID based wing vibration control system.
(a)
(b)


#### 2.2.2. Formulation of the Control Optimization Problem

The performance of the IMC-PID controller heavily depends on the design of the filter, especially for nonlinear system and MIMO system. The suitable filter parameter could have a fast response, less overshoot in the response and good robustness [33, 34]. The filter in the IMC-PID controller has the form as follows:(7)
f
s
,
λ
=
1
λ
s
+
1
n
,wheref
s is the low-pass filter with the filter parameter λ and the order n. The choice of order n should ensure the realization of the C
s as well as avoiding high order controller structure for the implementation. Since wing vibration systems can be approximated by second-order model [26], the order n is taken as 2 in this paper, as indicated in Section 2.2.1. More importantly, the filter parameter λ needs to be optimally designed, so as to obtain the satisfied dynamic performance for complex wing vibration control systems, including system nonlinearity, control constraints or multiple control surfaces. With the optimum of the filter parameter, the suitable PID parameters can be determined correspondingly.The optimal design problem of the filter parameter in IMC-PID is formulated as a constrained optimization problem by minimizing the vibration motions and the control effort of the control surface. The overall control performance of IMC-PID controller for wing vibration systems is evaluated by the objective functionJ, which is defined as the sum of normalized IAE based vibration responses and IAE based control signal response. IAE is chosen over the time multiplied integral criteria, such as ITAE and ITSE because they emphasize on the transient response in the later phase instead of the large error in the beginning. Then, the formulation of control optimization problem is presented as follows:(8)
min
J
=
I
A
E
=
∫
0
t
f
i
n
a
l
h
h
m
a
x
+
θ
θ
m
a
x
+
u
u
m
a
x
d
t
,
s
.
t
.
u
≤
u
m
a
x
u
˙
≤
u
r
a
t
e
λ
m
i
n
<
λ
<
λ
m
a
x
,whereh
m
a
x, θ
m
a
x are the maximum of the absolute of plunge displacement and pitch angle separately, u
m
a
x is the maximum of the control signal, t
f
i
n
a
l is the simulation time, u
m
a
x is the saturation of the control surface, u
r
a
t
e is the rate limit of the control surface, λ
m
a
x
,
λ
m
i
n are upper and lower limit of the filter parameter, respectively. The constraints for the optimization problem are handled as the death penalty [35], where the search is with feasible solutions.


## 2.1. System Description and Control Specifications

The wing vibration system is normally a rigid airfoil with pitch and plunge motions, which are manipulated by the flap control surface. The pitch and plunge motions are defined by pitch angleθ and plunge displacement h separately. Assuming the variable vector q
=
h
θ
T and the control signal of flap angle u
, general dynamic equations of the wing vibration system can be presented as:(1)
M
q
q
¨
+
C
q
q
˙
+
K
q
q
=
F
q
,
q
˙
,
u
,where the left side of the equation represents the structural dynamics of the wing vibration system with the mass matrixM, the damping matrix C and the spring matrix K. The right side of the equation describes the aerodynamics with respect to the vibration motions and the control command. As indicated in references [8, 12, 15], there are mainly three different types of the wing vibration system based on Equation (1). The first is linear wing vibration system with linear spring/damping characteristics and single control surface. The second is nonlinear wing vibration system with nonlinear spring/damping/aerodynamics and single control surface. The last is nonlinear MIMO wing vibration system, which has two control surfaces of both leading edge flap and trailing edge flap. The system nonlinearity and multiple control surfaces could further increase the difficulty in design of control parameters for the wing vibration system.The control specifications for wing vibration systems are:(1)
The wing vibration in pitch motion should be well suppressed.
(2)
The wing vibration in plunge motion should be well suppressed.
(3)
The satisfactory control performance should be obtained under control constraints.
(4)
The vibration control performance must be ensured under system nonlinearity.
(5)
The tear-and-wear of actuators should be minimized.
(6)
A sufficient closed-loop stability should be shown.In order to meet all the specific control requirements, the objective functions are usually given based on the vibration control performance, indicated by Integral of Absolute Error (IAE), Integral of Squared Error (ISE), Integral of Time multiplied Absolute Error (ITAE) and Integral of Time multiplied Square Error (ITSE). In this paper, IAE based objective function is taken as it emphasizes the instant response and the damping behavior, which are appropriate for the evaluation of the control performance for wing vibration systems.


## 2.2. Controller Structure and Problem Formulation

### 2.2.1. Internal Model Control (IMC) Based PID Structure

In wing vibration control system, PID controllers are mostly used to suppress the vibrations in plunge and pitch motions. However, the shortcomings of using the conventional PID controller include slow closed-loop responses, unsatisfied robustness to system change, inefficient utilization of the control surface and great control design effort. The reason is the difficulty in design of proper PID parameters, especially when complex system is involved, such as the nonlinear system and the MIMO system. An internal model control (IMC) based PID controller improves the settling time, overshoot in the transient response and the robustness of the classic PID controller, since the IMC can handle the model mismatch due to the nonlinearity or system change by the filter design and the process model [29, 30]. The filter in IMC-PID is designed to ensure the dynamic performance and improve the robustness of the PID controller under varied system parameters [31]. Meanwhile, three PID parameters can be automatically computed and updated by adjusting the filter parameter in the IMC-PID structure.To design the IMC controller for wing vibration systems, the process modelG
m
s of the system should first be obtained. In general, the process model of wing vibration systems can be approximated by the reduced order model G
R
O
M
sas:(2)
G
m
s
=
G
R
O
M
s
=
k
0
−
β
s
+
1
a
1
s
2
+
a
0
s
+
1
,wherek
0 is the system gain, β
,
a
1
,
a
0 are system parameters of the reduced order model. a
1
>
0
,
a
0
>
0. In this paper, the least-square method is adopted to identify the reduced order model of wing vibration systems.The control schematic of the IMC is shown in Figure1(a), where the feedback IMC controller is composed of the process model G
m
s and the controller C
s. The controller C
s is designed based on the filter f
s as [32]:(3)
C
s
=
f
s
G
m
−
s
,
G
m
s
=
G
m
−
s
G
m
+
s
,whereG
m
−
s
=
k
0
/
a
1
s
2
+
a
0
s
+
1is the minimum phase part and G
m
+
s
=
−
β
s
+
1 is the nonminimum phase part of the process model. The low-order filter f
s has the filter parameter λ, which can be adjusted to improve the dynamic control performance under system nonlinearity or system changes. In this paper, second-order filter is used as f
s
=
1
/
λ
s
+
1
2
.As indicated in Figure1, the feedback IMC controller G
c
s is presented as:(4)
G
c
s
=
C
s
1
−
G
m
s
C
s
.Substituting Equations (2) and (3) into Equation (4), the feedback IMC controller can be rewritten as:(5)
G
c
s
=
a
1
s
2
+
a
0
s
+
1
k
0
λ
2
s
2
+
2
λ
+
β
λ
s
.The equivalent IMC-PID controllerG
I
M
C
-
P
I
D
s can be derived from G
c
s as:(6)
G
I
M
C
-
P
I
D
s
=
G
c
s
≈
k
p
+
k
i
s
+
k
d
s
,
k
p
λ
=
−
1
2
k
0
λ
+
k
0
β
a
0
−
λ
2
2
λ
+
β
,
k
i
λ
=
−
1
2
k
0
λ
+
k
0
β
,
k
d
λ
=
−
1
2
k
0
λ
+
k
0
β
a
1
−
λ
2
2
λ
+
β
a
0
−
λ
2
2
λ
+
β
.For the design of IMC-PID controller for different kinds of wing vibration systems, the control parameters of the IMC-PID are obtained as followings.(1)
Linear wing vibration system: the process model of the original fourth order system is taken as a second-order model in Equation (2) by model identification. Then, the control parameters of the IMC-PID are derived based on Equations (3–6). The low-order process model is favored to avoid the high order of the controller.
(2)
Nonlinear wing vibration system: the process model of the nonlinear system is approximated by the linear reduced order model at operating points. The process model in Equation (2) is identified by the least-square method. Then, the IMC-PID control parameters are calculated according to Equations (3)–(6).
(3)
MIMO wing vibration system:The MIMO system usually has two system outputs and two control inputs. For MIMO system, two SISO control loops are formed, where one IMC-PID controller is designed independently for each control loop. For the controller design in each SISO control loop, the process model of the corresponding SISO system is identified as the reduced order model in Equation (2). Then, the parameters of the identified process model are used to calculate the parameters of the IMC-PID controller using Equations (3)–(6).It should be noted that the identified process model could vary according to different wing vibration systems. On this condition, the parameters of the corresponding IMC-PID controller could vary from case to case.As mentioned above, the implemented structure of IMC- PID based wing vibration control system is shown in Figure1(b). The vibration deflection, such as pitch angle is chosen as the system output y. The control input u is defined as the flap angle. The identified process model G
m
s and the designed filter f
λ are used to formulate the feedback IMC vibration controller G
c
s, which is further implemented in the form of PID controller. The control constraints, such as the saturation and rate limit are also considered in the control system. The filter is optimally designed by proposed GNPSO algorithm according to defined objective function J. The optimal design problem of the filter is presented in next section. The design of proposed GNPSO algorithm is detailed in Section 3.Figure 1
(a) Control schematic of Internal model control (IMC) system, (b) implemented structure of IMC-PID based wing vibration control system.
(a)
(b)


### 2.2.2. Formulation of the Control Optimization Problem

The performance of the IMC-PID controller heavily depends on the design of the filter, especially for nonlinear system and MIMO system. The suitable filter parameter could have a fast response, less overshoot in the response and good robustness [33, 34]. The filter in the IMC-PID controller has the form as follows:(7)
f
s
,
λ
=
1
λ
s
+
1
n
,wheref
s is the low-pass filter with the filter parameter λ and the order n. The choice of order n should ensure the realization of the C
s as well as avoiding high order controller structure for the implementation. Since wing vibration systems can be approximated by second-order model [26], the order n is taken as 2 in this paper, as indicated in Section 2.2.1. More importantly, the filter parameter λ needs to be optimally designed, so as to obtain the satisfied dynamic performance for complex wing vibration control systems, including system nonlinearity, control constraints or multiple control surfaces. With the optimum of the filter parameter, the suitable PID parameters can be determined correspondingly.The optimal design problem of the filter parameter in IMC-PID is formulated as a constrained optimization problem by minimizing the vibration motions and the control effort of the control surface. The overall control performance of IMC-PID controller for wing vibration systems is evaluated by the objective functionJ, which is defined as the sum of normalized IAE based vibration responses and IAE based control signal response. IAE is chosen over the time multiplied integral criteria, such as ITAE and ITSE because they emphasize on the transient response in the later phase instead of the large error in the beginning. Then, the formulation of control optimization problem is presented as follows:(8)
min
J
=
I
A
E
=
∫
0
t
f
i
n
a
l
h
h
m
a
x
+
θ
θ
m
a
x
+
u
u
m
a
x
d
t
,
s
.
t
.
u
≤
u
m
a
x
u
˙
≤
u
r
a
t
e
λ
m
i
n
<
λ
<
λ
m
a
x
,whereh
m
a
x, θ
m
a
x are the maximum of the absolute of plunge displacement and pitch angle separately, u
m
a
x is the maximum of the control signal, t
f
i
n
a
l is the simulation time, u
m
a
x is the saturation of the control surface, u
r
a
t
e is the rate limit of the control surface, λ
m
a
x
,
λ
m
i
n are upper and lower limit of the filter parameter, respectively. The constraints for the optimization problem are handled as the death penalty [35], where the search is with feasible solutions.


## 2.2.1. Internal Model Control (IMC) Based PID Structure

In wing vibration control system, PID controllers are mostly used to suppress the vibrations in plunge and pitch motions. However, the shortcomings of using the conventional PID controller include slow closed-loop responses, unsatisfied robustness to system change, inefficient utilization of the control surface and great control design effort. The reason is the difficulty in design of proper PID parameters, especially when complex system is involved, such as the nonlinear system and the MIMO system. An internal model control (IMC) based PID controller improves the settling time, overshoot in the transient response and the robustness of the classic PID controller, since the IMC can handle the model mismatch due to the nonlinearity or system change by the filter design and the process model [29, 30]. The filter in IMC-PID is designed to ensure the dynamic performance and improve the robustness of the PID controller under varied system parameters [31]. Meanwhile, three PID parameters can be automatically computed and updated by adjusting the filter parameter in the IMC-PID structure.To design the IMC controller for wing vibration systems, the process modelG
m
s of the system should first be obtained. In general, the process model of wing vibration systems can be approximated by the reduced order model G
R
O
M
sas:(2)
G
m
s
=
G
R
O
M
s
=
k
0
−
β
s
+
1
a
1
s
2
+
a
0
s
+
1
,wherek
0 is the system gain, β
,
a
1
,
a
0 are system parameters of the reduced order model. a
1
>
0
,
a
0
>
0. In this paper, the least-square method is adopted to identify the reduced order model of wing vibration systems.The control schematic of the IMC is shown in Figure1(a), where the feedback IMC controller is composed of the process model G
m
s and the controller C
s. The controller C
s is designed based on the filter f
s as [32]:(3)
C
s
=
f
s
G
m
−
s
,
G
m
s
=
G
m
−
s
G
m
+
s
,whereG
m
−
s
=
k
0
/
a
1
s
2
+
a
0
s
+
1is the minimum phase part and G
m
+
s
=
−
β
s
+
1 is the nonminimum phase part of the process model. The low-order filter f
s has the filter parameter λ, which can be adjusted to improve the dynamic control performance under system nonlinearity or system changes. In this paper, second-order filter is used as f
s
=
1
/
λ
s
+
1
2
.As indicated in Figure1, the feedback IMC controller G
c
s is presented as:(4)
G
c
s
=
C
s
1
−
G
m
s
C
s
.Substituting Equations (2) and (3) into Equation (4), the feedback IMC controller can be rewritten as:(5)
G
c
s
=
a
1
s
2
+
a
0
s
+
1
k
0
λ
2
s
2
+
2
λ
+
β
λ
s
.The equivalent IMC-PID controllerG
I
M
C
-
P
I
D
s can be derived from G
c
s as:(6)
G
I
M
C
-
P
I
D
s
=
G
c
s
≈
k
p
+
k
i
s
+
k
d
s
,
k
p
λ
=
−
1
2
k
0
λ
+
k
0
β
a
0
−
λ
2
2
λ
+
β
,
k
i
λ
=
−
1
2
k
0
λ
+
k
0
β
,
k
d
λ
=
−
1
2
k
0
λ
+
k
0
β
a
1
−
λ
2
2
λ
+
β
a
0
−
λ
2
2
λ
+
β
.For the design of IMC-PID controller for different kinds of wing vibration systems, the control parameters of the IMC-PID are obtained as followings.(1)
Linear wing vibration system: the process model of the original fourth order system is taken as a second-order model in Equation (2) by model identification. Then, the control parameters of the IMC-PID are derived based on Equations (3–6). The low-order process model is favored to avoid the high order of the controller.
(2)
Nonlinear wing vibration system: the process model of the nonlinear system is approximated by the linear reduced order model at operating points. The process model in Equation (2) is identified by the least-square method. Then, the IMC-PID control parameters are calculated according to Equations (3)–(6).
(3)
MIMO wing vibration system:The MIMO system usually has two system outputs and two control inputs. For MIMO system, two SISO control loops are formed, where one IMC-PID controller is designed independently for each control loop. For the controller design in each SISO control loop, the process model of the corresponding SISO system is identified as the reduced order model in Equation (2). Then, the parameters of the identified process model are used to calculate the parameters of the IMC-PID controller using Equations (3)–(6).It should be noted that the identified process model could vary according to different wing vibration systems. On this condition, the parameters of the corresponding IMC-PID controller could vary from case to case.As mentioned above, the implemented structure of IMC- PID based wing vibration control system is shown in Figure1(b). The vibration deflection, such as pitch angle is chosen as the system output y. The control input u is defined as the flap angle. The identified process model G
m
s and the designed filter f
λ are used to formulate the feedback IMC vibration controller G
c
s, which is further implemented in the form of PID controller. The control constraints, such as the saturation and rate limit are also considered in the control system. The filter is optimally designed by proposed GNPSO algorithm according to defined objective function J. The optimal design problem of the filter is presented in next section. The design of proposed GNPSO algorithm is detailed in Section 3.Figure 1
(a) Control schematic of Internal model control (IMC) system, (b) implemented structure of IMC-PID based wing vibration control system.
(a)
(b)


## 2.2.2. Formulation of the Control Optimization Problem

The performance of the IMC-PID controller heavily depends on the design of the filter, especially for nonlinear system and MIMO system. The suitable filter parameter could have a fast response, less overshoot in the response and good robustness [33, 34]. The filter in the IMC-PID controller has the form as follows:(7)
f
s
,
λ
=
1
λ
s
+
1
n
,wheref
s is the low-pass filter with the filter parameter λ and the order n. The choice of order n should ensure the realization of the C
s as well as avoiding high order controller structure for the implementation. Since wing vibration systems can be approximated by second-order model [26], the order n is taken as 2 in this paper, as indicated in Section 2.2.1. More importantly, the filter parameter λ needs to be optimally designed, so as to obtain the satisfied dynamic performance for complex wing vibration control systems, including system nonlinearity, control constraints or multiple control surfaces. With the optimum of the filter parameter, the suitable PID parameters can be determined correspondingly.The optimal design problem of the filter parameter in IMC-PID is formulated as a constrained optimization problem by minimizing the vibration motions and the control effort of the control surface. The overall control performance of IMC-PID controller for wing vibration systems is evaluated by the objective functionJ, which is defined as the sum of normalized IAE based vibration responses and IAE based control signal response. IAE is chosen over the time multiplied integral criteria, such as ITAE and ITSE because they emphasize on the transient response in the later phase instead of the large error in the beginning. Then, the formulation of control optimization problem is presented as follows:(8)
min
J
=
I
A
E
=
∫
0
t
f
i
n
a
l
h
h
m
a
x
+
θ
θ
m
a
x
+
u
u
m
a
x
d
t
,
s
.
t
.
u
≤
u
m
a
x
u
˙
≤
u
r
a
t
e
λ
m
i
n
<
λ
<
λ
m
a
x
,whereh
m
a
x, θ
m
a
x are the maximum of the absolute of plunge displacement and pitch angle separately, u
m
a
x is the maximum of the control signal, t
f
i
n
a
l is the simulation time, u
m
a
x is the saturation of the control surface, u
r
a
t
e is the rate limit of the control surface, λ
m
a
x
,
λ
m
i
n are upper and lower limit of the filter parameter, respectively. The constraints for the optimization problem are handled as the death penalty [35], where the search is with feasible solutions.


## 3. Proposed Improved Grey-PSO and NLJ Hybrid Optimization Technique

In this section, a hybrid algorithm of improved Grey-PSO and NLJ, called GNPS is proposed to solve the IMC-PID control optimization problem in Equation (8). The initiatives of proposed hybrid algorithm are to search for satisfying optimal IMC-PID parameters at low computation cost and to have wide applications in diverse wing vibration systems, from simple linear systems to complex nonlinear systems. The preliminaries and the main GNPSO algorithm are presented in the following subsections.


### 3.1. Classic Grey-PSO Algorithm

Particle swarm optimization (PSO) proposed by Eberhart and Kennedy is a population cooperation based stochastic searching algorithm [36]. It is inspired by the behavior of bird flocking to guide the particles to search for the global optimum by evolving velocity. The trajectory of each particle is directed by dynamically adjusting the velocity of the particle, with respect to its own flying experience (P
best) and the fly experience of the swarm (S
best). The merits of the PSO include strong global exploration ability and simple evolutionary process without the crossover operation and the mutation operation [37]. The PSO has successful and wide applications in different areas, such as system identification [38–40], automotive system [41, 42] and robot systems [43–45].For the past decade, several works about modified PSO have been investigated for a different optimization problem. In the study of optimized PID control, a novel modified PSO algorithm based on grey relational analysis theory (called GPSO) has recently been proposed by Professor Yeh [24]. Compared to conventional algorithms, such as GA and PSO, the GPSO optimization technique can obtain better optimal parameters of PID controller and have low computation burden with reduced population size. For the implementation of the GPSO, the “similarity” (the relationship) between the particle and the fittest particle is first analyzed. Based on the closeness of the relationship, algorithm parameters, such as acceleration coefficients are properly adjusted to improve the guidance of particles’ movement.In classic GPSO algorithm, the relationship betweenith particle x
i
d and the fittest particle S
d
b
e
s
t is indicated by grey relational grade g
i according to classic grey relational analysis theory:(9)
g
i
=
1
D
∑
d
=
1
D
γ
S
d
b
e
s
t
,
x
i
d
,
γ
S
d
b
e
s
t
,
x
i
d
=
Δ
m
i
n
+
ζ
Δ
m
a
x
Δ
i
d
+
ζ
Δ
m
a
x
,whereΔ
i
d
=
S
d
b
e
s
t
−
x
i
d
,
Δ
m
i
n
=
m
i
n
i
m
i
n
d
Δ
i
d
,
Δ
m
a
x
=
m
a
x
i
m
a
x
d
Δ
i
d and ζ is the scalar between 0 and 1, d is the dimension of the solution space.Based on grey relational gradeg
i, the updating law of acceleration coefficients for each particle is presented as:(10)
c
1
i
=
0.5
cos
f
g
i
π
+
2
,
c
2
i
=
4
−
c
1
i
.wheref
g
i
=
1.5
1
+
ζ
g
i
+
1.5
1
−
ζ, and the sum of c
1
i and c
2
i is fixed for each particle in every iteration during the evolutionary process.The updating rule for the velocity of the particleV
i
d is given as:(11)
V
i
d
=
ω
i
V
i
d
+
c
1
i
r
a
n
d
1
P
i
d
b
e
s
t
−
x
i
d
+
c
2
i
r
a
n
d
2
S
d
b
e
s
t
−
x
i
d
,whereω
i is the inertia weight, and r
a
n
d
1, r
a
n
d
2 are the random numbers.


### 3.2. Hybrid GNPSO Algorithm

#### 3.2.1. Improved Grey-PSO (IGPSO) Optimization Algorithm

For complex nonlinear wing vibration systems, the control optimization could cost much more computation time than that for simple linear low-order system using GPSO in [24]. On this condition, compared to original GPSO algorithm, further reduced population size and iteration number are preferred to avoid large computation time. Moreover, the modifications in grey analysis and updating formula of acceleration coefficients are also considered for the control optimization problem in this paper. Therefore, three improvements are added in the original GPSO to enhance its effectiveness, efficiency, and applicability for optimal IMC-PID control issue of wing vibration systems. Those improvements are detailed as follows.(1)Small population size and small iteration number: In a classic heuristic optimization algorithm, such as GA, PSO and other hybrid PSO algorithms, there are usually large population size and large iteration number to ensure the optimization performance. Large computation time is inevitable on those conditions, especially for complex wing vibration control systems. However, fast tuning of control parameters is favored by most control systems since the uncontrolled system could be damaged. The large time delay of the control signal could greatly degrade the control performance and even fail the stabilization of the system in practice. In this paper, compared to those for GPSO optimized PID control of simple linear system, smaller population size P and smaller iteration number G
m
a
x are given to avoid large computation time for the IMC-PID tuning for complex wing vibration systems.
(2)New grey analysis rule: In the original GPSO, the “similarity” (the relationship) between the particle and the fittest particle is actually analyzed based on the space distance between them using classic grey relational analysis theory [46]. However, the parameters of PID controller have to be in a suitable proportional relationship, so as to obtain satisfied control performance. The closeness in the space cannot guarantee the proper relationship between PID parameters. In the proposed improved GPSO, the geometric similarity instead of space distance is chosen as the criterion for the analysis of particles’ relationship. Then, the accuracy of grey analysis for IMC-PID control optimization can be improved. The new grey analysis rule based on the absolute grey relational analysis (AGRA), is proposed to find the geometric similarity between the particle x
i and the fittest particle S
best. The result of the AGRA for ith particle is indicated by the grey relational grade g
i:(12)
g
i
=
1
+
S
best
,
0
+
x
i
0
1
+
S
best
,
0
+
x
i
0
+
x
i
0
−
S
best
,
0
,
whereS
best
,
0 and x
i
0 are normalized S
best and x
i with zero start point. S
best
,
0
,
x
0 and x
i
0
−
S
best
,
0 are specified as:(13)
S
best
,
0
=
∑
d
=
2
D
−
1
g
B
e
s
t
d
0
+
1
2
S
D
b
e
s
t
,
0
,
x
0
=
∑
d
=
2
D
−
1
x
d
0
+
1
2
x
D
0
,
x
i
0
−
S
best
,
0
=
∑
d
=
2
D
−
1
x
d
0
−
S
d
b
e
s
t
,
0
+
1
2
x
D
0
−
S
D
b
e
s
t
,
0
.
(3)New updating formula of acceleration coefficients: In the original GPSO algorithm, the acceleration coefficient c
2
i changes directly by c
1
i based on the fixed sum, which gives the limited space to adjust the coefficients. Besides, the updating formula of c
1
i moves in a sinusoid way. In proposed algorithm, the updating law of acceleration coefficients is improved by efficient computation of linear functions and flexible adjusting space for enhanced diversity of the coefficients. Specifically, c
1
i and c
2
i are updated independently with un-fixed sum, which could expand the range of the coefficients to fit different applications. The new updating formula of acceleration coefficients are given as:(14)
c
1
i
=
−
g
i
c
1
m
−
c
1
n
+
c
1
m
−
c
1
n
,
(15)
c
2
i
=
g
i
c
2
m
−
c
2
n
2
+
2
c
2
m
−
c
2
n
2
,
(16)
4
≥
c
1
m
>
c
1
n
≥
0
,
3
≥
c
2
m
>
c
2
n
≥
0
,
wherec
1
m
,
c
1
n
,
c
2
m
,
c
2
n are scalars to determine the range of c
1
i, c
2
i and their sum. The altering of c
2
i is not limited by the choice of c
1
i. They are updated based on the closeness (geometric similarity) of the particle and the fittest particle using AGRA: the small g
i means the evolution process of the particle in the phase of the exploration, when increased c
1
i and decreased c
2
i are given to enhance the search of local optimum and maintain the swarm diversity. Otherwise, large g
i indicates the particle is much close to the global best position. Smaller c
1
i and larger c
2
i are assigned to enhance the local search and exploit the global optimum. Based on Equations (14–16), the sum of acceleration coefficients is flexible with varied range. In this way, to search for the global optimum for different types of wing vibration control systems, the suitable coefficients could be adjusted in different zones.


#### 3.2.2. NLJ Algorithm

The new Luus-Jaakola (NLJ) algorithm an improved Luus-Jaakola algorithm, proposed by Pan [47]. It has fast convergence capability for the local search [48]. The advantage of the NLJ lies in the adaptation of the convergence coefficient, which leads to the rapid reduction of the search horizon and significantly saves the search time. The NLJ algorithm begins with an initial vector a
0 and initial search radius r
0. The search range at the generation k is determined by the optimization result at previous generation as r
k
=
c
r
a
k
−
1
v
k
−
1, where ν is determined by the convergence coefficient φ
k as v
k
=
φ
k
v
k
−
1
,
if
a
k
/
a
k
−
1
<
η
,
k
≥
2, otherwise v
k
=
1. Then, the jth vector of the population at current generation is updated based on the improved search radius as a
j
k
=
a
j
k
−
1
+
r
a
n
d
j
r
k. Based on the evaluation of performance index, the optimal result at current generation is obtained as a
k. Keep changing the search range after taking current optimum as the initial value of the next generation until the maximum generation number is reached. The coefficients of the NLJ algorithm are detailed in Refs. [22, 23].


#### 3.2.3. Hybridization of IGPSO and NLJ Algorithm

In the proposed GNPSO hybrid algorithm, the improved Grey-PSO (IGPSO) algorithm is first applied to optimize the parameters of IMC-PID controller. Then, the NLJ is further used to enhance the exploitation and ensure the convergence to global optimum. The Grey-PSO algorithm has the strength of global search while NLJ is a local search approach. Then, a hybrid optimization method based on those techniques could be promising and have benefits from their merits. The hybrid GNPSO technique is implemented with small population size and small maximum generation number, compared to other algorithms for PID optimization in references [20–24].In this way, proposed hybrid GNPSO algorithm could obtain: (1) improved global search by IGPSO to avoid arbitrary and inefficient search under the complexity of nonlinear wing vibration systems; (2) enhanced local search by NLJ to fasten convergence and avoid sub-optimum; (3) good adaptation to different wing vibration control systems by improving dynamic balance between the exploration and the exploitation; and (4) small computation cost due to small population size and small iteration number.


### 3.3. Procedures of Hybrid GNPSO Based Control Optimization

Figure2 shows control optimization procedures by proposed GNPSO hybrid algorithm. The GNPSO optimization process includes the following steps:Figure 2
(a) Optimization procedure by proposed GNPSO algorithm. (b) Procedures of calculating distance by new grey analysis rule. (c) Procedures of updating acceleration coefficients.
(a)
(b)
(c)Step (1) Initialization of IGPSO: Set the algorithm parameters, such as population size and iteration number of IGPSO and Initialize the population of particles.
Step (2) Evaluation of objective function:Evaluate the fitness of the particle in the population and find the fittest particle in the current generation. The evaluation is made by assigning the value of particle to control parameters, running the control system and calculating the objective function J in Equation (8). The fittest particle is obtained with minimized J.
Step (3) Calculate the distance between particle and fittest one: Calculate the grey relational grade g
i, which indicates the distance between the particle and the fittest particle. Based on the analysis of geometric similarity between the particle and fittest particle, g
i is obtained by new grey analysis rule using Equations (12) and (13).
Step (4) Update acceleration coefficients:Set scalars in Equation (16) to define the range and the sum of acceleration coefficients c
1
i
,
c
2
i. Then, those acceleration coefficients are adjusted by the new updating formula using Equations (14) and (15).
Step (5) Update velocity and position of particle: Based on the adjusted acceleration coefficients in step (4), update the velocity and the position of the particle.
Step (6) Output optimum of IGPSO search: If the maximum iteration number is not reached, go back to (2). Otherwise, output the best result obtained by IGPSO search in step (1)–step (5).
Step (7) Initialization of NLJ: Define the optimum of IGPSO search in step (6) as initial points of the NLJ optimization. The population of NLJ is initialized with initial search space.
Step (8) NLJ search: The search space is narrowed to improve the local search of the GNPSO. The new population is generated with narrowed search space. Then, the fittest one among the population is found with minimized objective function.
Step (9) Output global optimum: if the condition, such as iteration number is not reached, go back to Otherwise, output global optimum by NLJ search in GNPSO algorithm. The global optimum is the final optimized control parameter of IMC-PID controller for the wing vibration system.


## 3.1. Classic Grey-PSO Algorithm

Particle swarm optimization (PSO) proposed by Eberhart and Kennedy is a population cooperation based stochastic searching algorithm [36]. It is inspired by the behavior of bird flocking to guide the particles to search for the global optimum by evolving velocity. The trajectory of each particle is directed by dynamically adjusting the velocity of the particle, with respect to its own flying experience (P
best) and the fly experience of the swarm (S
best). The merits of the PSO include strong global exploration ability and simple evolutionary process without the crossover operation and the mutation operation [37]. The PSO has successful and wide applications in different areas, such as system identification [38–40], automotive system [41, 42] and robot systems [43–45].For the past decade, several works about modified PSO have been investigated for a different optimization problem. In the study of optimized PID control, a novel modified PSO algorithm based on grey relational analysis theory (called GPSO) has recently been proposed by Professor Yeh [24]. Compared to conventional algorithms, such as GA and PSO, the GPSO optimization technique can obtain better optimal parameters of PID controller and have low computation burden with reduced population size. For the implementation of the GPSO, the “similarity” (the relationship) between the particle and the fittest particle is first analyzed. Based on the closeness of the relationship, algorithm parameters, such as acceleration coefficients are properly adjusted to improve the guidance of particles’ movement.In classic GPSO algorithm, the relationship betweenith particle x
i
d and the fittest particle S
d
b
e
s
t is indicated by grey relational grade g
i according to classic grey relational analysis theory:(9)
g
i
=
1
D
∑
d
=
1
D
γ
S
d
b
e
s
t
,
x
i
d
,
γ
S
d
b
e
s
t
,
x
i
d
=
Δ
m
i
n
+
ζ
Δ
m
a
x
Δ
i
d
+
ζ
Δ
m
a
x
,whereΔ
i
d
=
S
d
b
e
s
t
−
x
i
d
,
Δ
m
i
n
=
m
i
n
i
m
i
n
d
Δ
i
d
,
Δ
m
a
x
=
m
a
x
i
m
a
x
d
Δ
i
d and ζ is the scalar between 0 and 1, d is the dimension of the solution space.Based on grey relational gradeg
i, the updating law of acceleration coefficients for each particle is presented as:(10)
c
1
i
=
0.5
cos
f
g
i
π
+
2
,
c
2
i
=
4
−
c
1
i
.wheref
g
i
=
1.5
1
+
ζ
g
i
+
1.5
1
−
ζ, and the sum of c
1
i and c
2
i is fixed for each particle in every iteration during the evolutionary process.The updating rule for the velocity of the particleV
i
d is given as:(11)
V
i
d
=
ω
i
V
i
d
+
c
1
i
r
a
n
d
1
P
i
d
b
e
s
t
−
x
i
d
+
c
2
i
r
a
n
d
2
S
d
b
e
s
t
−
x
i
d
,whereω
i is the inertia weight, and r
a
n
d
1, r
a
n
d
2 are the random numbers.


## 3.2. Hybrid GNPSO Algorithm

### 3.2.1. Improved Grey-PSO (IGPSO) Optimization Algorithm

For complex nonlinear wing vibration systems, the control optimization could cost much more computation time than that for simple linear low-order system using GPSO in [24]. On this condition, compared to original GPSO algorithm, further reduced population size and iteration number are preferred to avoid large computation time. Moreover, the modifications in grey analysis and updating formula of acceleration coefficients are also considered for the control optimization problem in this paper. Therefore, three improvements are added in the original GPSO to enhance its effectiveness, efficiency, and applicability for optimal IMC-PID control issue of wing vibration systems. Those improvements are detailed as follows.(1)Small population size and small iteration number: In a classic heuristic optimization algorithm, such as GA, PSO and other hybrid PSO algorithms, there are usually large population size and large iteration number to ensure the optimization performance. Large computation time is inevitable on those conditions, especially for complex wing vibration control systems. However, fast tuning of control parameters is favored by most control systems since the uncontrolled system could be damaged. The large time delay of the control signal could greatly degrade the control performance and even fail the stabilization of the system in practice. In this paper, compared to those for GPSO optimized PID control of simple linear system, smaller population size P and smaller iteration number G
m
a
x are given to avoid large computation time for the IMC-PID tuning for complex wing vibration systems.
(2)New grey analysis rule: In the original GPSO, the “similarity” (the relationship) between the particle and the fittest particle is actually analyzed based on the space distance between them using classic grey relational analysis theory [46]. However, the parameters of PID controller have to be in a suitable proportional relationship, so as to obtain satisfied control performance. The closeness in the space cannot guarantee the proper relationship between PID parameters. In the proposed improved GPSO, the geometric similarity instead of space distance is chosen as the criterion for the analysis of particles’ relationship. Then, the accuracy of grey analysis for IMC-PID control optimization can be improved. The new grey analysis rule based on the absolute grey relational analysis (AGRA), is proposed to find the geometric similarity between the particle x
i and the fittest particle S
best. The result of the AGRA for ith particle is indicated by the grey relational grade g
i:(12)
g
i
=
1
+
S
best
,
0
+
x
i
0
1
+
S
best
,
0
+
x
i
0
+
x
i
0
−
S
best
,
0
,
whereS
best
,
0 and x
i
0 are normalized S
best and x
i with zero start point. S
best
,
0
,
x
0 and x
i
0
−
S
best
,
0 are specified as:(13)
S
best
,
0
=
∑
d
=
2
D
−
1
g
B
e
s
t
d
0
+
1
2
S
D
b
e
s
t
,
0
,
x
0
=
∑
d
=
2
D
−
1
x
d
0
+
1
2
x
D
0
,
x
i
0
−
S
best
,
0
=
∑
d
=
2
D
−
1
x
d
0
−
S
d
b
e
s
t
,
0
+
1
2
x
D
0
−
S
D
b
e
s
t
,
0
.
(3)New updating formula of acceleration coefficients: In the original GPSO algorithm, the acceleration coefficient c
2
i changes directly by c
1
i based on the fixed sum, which gives the limited space to adjust the coefficients. Besides, the updating formula of c
1
i moves in a sinusoid way. In proposed algorithm, the updating law of acceleration coefficients is improved by efficient computation of linear functions and flexible adjusting space for enhanced diversity of the coefficients. Specifically, c
1
i and c
2
i are updated independently with un-fixed sum, which could expand the range of the coefficients to fit different applications. The new updating formula of acceleration coefficients are given as:(14)
c
1
i
=
−
g
i
c
1
m
−
c
1
n
+
c
1
m
−
c
1
n
,
(15)
c
2
i
=
g
i
c
2
m
−
c
2
n
2
+
2
c
2
m
−
c
2
n
2
,
(16)
4
≥
c
1
m
>
c
1
n
≥
0
,
3
≥
c
2
m
>
c
2
n
≥
0
,
wherec
1
m
,
c
1
n
,
c
2
m
,
c
2
n are scalars to determine the range of c
1
i, c
2
i and their sum. The altering of c
2
i is not limited by the choice of c
1
i. They are updated based on the closeness (geometric similarity) of the particle and the fittest particle using AGRA: the small g
i means the evolution process of the particle in the phase of the exploration, when increased c
1
i and decreased c
2
i are given to enhance the search of local optimum and maintain the swarm diversity. Otherwise, large g
i indicates the particle is much close to the global best position. Smaller c
1
i and larger c
2
i are assigned to enhance the local search and exploit the global optimum. Based on Equations (14–16), the sum of acceleration coefficients is flexible with varied range. In this way, to search for the global optimum for different types of wing vibration control systems, the suitable coefficients could be adjusted in different zones.


### 3.2.2. NLJ Algorithm

The new Luus-Jaakola (NLJ) algorithm an improved Luus-Jaakola algorithm, proposed by Pan [47]. It has fast convergence capability for the local search [48]. The advantage of the NLJ lies in the adaptation of the convergence coefficient, which leads to the rapid reduction of the search horizon and significantly saves the search time. The NLJ algorithm begins with an initial vector a
0 and initial search radius r
0. The search range at the generation k is determined by the optimization result at previous generation as r
k
=
c
r
a
k
−
1
v
k
−
1, where ν is determined by the convergence coefficient φ
k as v
k
=
φ
k
v
k
−
1
,
if
a
k
/
a
k
−
1
<
η
,
k
≥
2, otherwise v
k
=
1. Then, the jth vector of the population at current generation is updated based on the improved search radius as a
j
k
=
a
j
k
−
1
+
r
a
n
d
j
r
k. Based on the evaluation of performance index, the optimal result at current generation is obtained as a
k. Keep changing the search range after taking current optimum as the initial value of the next generation until the maximum generation number is reached. The coefficients of the NLJ algorithm are detailed in Refs. [22, 23].


### 3.2.3. Hybridization of IGPSO and NLJ Algorithm

In the proposed GNPSO hybrid algorithm, the improved Grey-PSO (IGPSO) algorithm is first applied to optimize the parameters of IMC-PID controller. Then, the NLJ is further used to enhance the exploitation and ensure the convergence to global optimum. The Grey-PSO algorithm has the strength of global search while NLJ is a local search approach. Then, a hybrid optimization method based on those techniques could be promising and have benefits from their merits. The hybrid GNPSO technique is implemented with small population size and small maximum generation number, compared to other algorithms for PID optimization in references [20–24].In this way, proposed hybrid GNPSO algorithm could obtain: (1) improved global search by IGPSO to avoid arbitrary and inefficient search under the complexity of nonlinear wing vibration systems; (2) enhanced local search by NLJ to fasten convergence and avoid sub-optimum; (3) good adaptation to different wing vibration control systems by improving dynamic balance between the exploration and the exploitation; and (4) small computation cost due to small population size and small iteration number.


## 3.2.1. Improved Grey-PSO (IGPSO) Optimization Algorithm

For complex nonlinear wing vibration systems, the control optimization could cost much more computation time than that for simple linear low-order system using GPSO in [24]. On this condition, compared to original GPSO algorithm, further reduced population size and iteration number are preferred to avoid large computation time. Moreover, the modifications in grey analysis and updating formula of acceleration coefficients are also considered for the control optimization problem in this paper. Therefore, three improvements are added in the original GPSO to enhance its effectiveness, efficiency, and applicability for optimal IMC-PID control issue of wing vibration systems. Those improvements are detailed as follows.(1)Small population size and small iteration number: In a classic heuristic optimization algorithm, such as GA, PSO and other hybrid PSO algorithms, there are usually large population size and large iteration number to ensure the optimization performance. Large computation time is inevitable on those conditions, especially for complex wing vibration control systems. However, fast tuning of control parameters is favored by most control systems since the uncontrolled system could be damaged. The large time delay of the control signal could greatly degrade the control performance and even fail the stabilization of the system in practice. In this paper, compared to those for GPSO optimized PID control of simple linear system, smaller population size P and smaller iteration number G
m
a
x are given to avoid large computation time for the IMC-PID tuning for complex wing vibration systems.
(2)New grey analysis rule: In the original GPSO, the “similarity” (the relationship) between the particle and the fittest particle is actually analyzed based on the space distance between them using classic grey relational analysis theory [46]. However, the parameters of PID controller have to be in a suitable proportional relationship, so as to obtain satisfied control performance. The closeness in the space cannot guarantee the proper relationship between PID parameters. In the proposed improved GPSO, the geometric similarity instead of space distance is chosen as the criterion for the analysis of particles’ relationship. Then, the accuracy of grey analysis for IMC-PID control optimization can be improved. The new grey analysis rule based on the absolute grey relational analysis (AGRA), is proposed to find the geometric similarity between the particle x
i and the fittest particle S
best. The result of the AGRA for ith particle is indicated by the grey relational grade g
i:(12)
g
i
=
1
+
S
best
,
0
+
x
i
0
1
+
S
best
,
0
+
x
i
0
+
x
i
0
−
S
best
,
0
,
whereS
best
,
0 and x
i
0 are normalized S
best and x
i with zero start point. S
best
,
0
,
x
0 and x
i
0
−
S
best
,
0 are specified as:(13)
S
best
,
0
=
∑
d
=
2
D
−
1
g
B
e
s
t
d
0
+
1
2
S
D
b
e
s
t
,
0
,
x
0
=
∑
d
=
2
D
−
1
x
d
0
+
1
2
x
D
0
,
x
i
0
−
S
best
,
0
=
∑
d
=
2
D
−
1
x
d
0
−
S
d
b
e
s
t
,
0
+
1
2
x
D
0
−
S
D
b
e
s
t
,
0
.
(3)New updating formula of acceleration coefficients: In the original GPSO algorithm, the acceleration coefficient c
2
i changes directly by c
1
i based on the fixed sum, which gives the limited space to adjust the coefficients. Besides, the updating formula of c
1
i moves in a sinusoid way. In proposed algorithm, the updating law of acceleration coefficients is improved by efficient computation of linear functions and flexible adjusting space for enhanced diversity of the coefficients. Specifically, c
1
i and c
2
i are updated independently with un-fixed sum, which could expand the range of the coefficients to fit different applications. The new updating formula of acceleration coefficients are given as:(14)
c
1
i
=
−
g
i
c
1
m
−
c
1
n
+
c
1
m
−
c
1
n
,
(15)
c
2
i
=
g
i
c
2
m
−
c
2
n
2
+
2
c
2
m
−
c
2
n
2
,
(16)
4
≥
c
1
m
>
c
1
n
≥
0
,
3
≥
c
2
m
>
c
2
n
≥
0
,
wherec
1
m
,
c
1
n
,
c
2
m
,
c
2
n are scalars to determine the range of c
1
i, c
2
i and their sum. The altering of c
2
i is not limited by the choice of c
1
i. They are updated based on the closeness (geometric similarity) of the particle and the fittest particle using AGRA: the small g
i means the evolution process of the particle in the phase of the exploration, when increased c
1
i and decreased c
2
i are given to enhance the search of local optimum and maintain the swarm diversity. Otherwise, large g
i indicates the particle is much close to the global best position. Smaller c
1
i and larger c
2
i are assigned to enhance the local search and exploit the global optimum. Based on Equations (14–16), the sum of acceleration coefficients is flexible with varied range. In this way, to search for the global optimum for different types of wing vibration control systems, the suitable coefficients could be adjusted in different zones.


## 3.2.2. NLJ Algorithm

The new Luus-Jaakola (NLJ) algorithm an improved Luus-Jaakola algorithm, proposed by Pan [47]. It has fast convergence capability for the local search [48]. The advantage of the NLJ lies in the adaptation of the convergence coefficient, which leads to the rapid reduction of the search horizon and significantly saves the search time. The NLJ algorithm begins with an initial vector a
0 and initial search radius r
0. The search range at the generation k is determined by the optimization result at previous generation as r
k
=
c
r
a
k
−
1
v
k
−
1, where ν is determined by the convergence coefficient φ
k as v
k
=
φ
k
v
k
−
1
,
if
a
k
/
a
k
−
1
<
η
,
k
≥
2, otherwise v
k
=
1. Then, the jth vector of the population at current generation is updated based on the improved search radius as a
j
k
=
a
j
k
−
1
+
r
a
n
d
j
r
k. Based on the evaluation of performance index, the optimal result at current generation is obtained as a
k. Keep changing the search range after taking current optimum as the initial value of the next generation until the maximum generation number is reached. The coefficients of the NLJ algorithm are detailed in Refs. [22, 23].


## 3.2.3. Hybridization of IGPSO and NLJ Algorithm

In the proposed GNPSO hybrid algorithm, the improved Grey-PSO (IGPSO) algorithm is first applied to optimize the parameters of IMC-PID controller. Then, the NLJ is further used to enhance the exploitation and ensure the convergence to global optimum. The Grey-PSO algorithm has the strength of global search while NLJ is a local search approach. Then, a hybrid optimization method based on those techniques could be promising and have benefits from their merits. The hybrid GNPSO technique is implemented with small population size and small maximum generation number, compared to other algorithms for PID optimization in references [20–24].In this way, proposed hybrid GNPSO algorithm could obtain: (1) improved global search by IGPSO to avoid arbitrary and inefficient search under the complexity of nonlinear wing vibration systems; (2) enhanced local search by NLJ to fasten convergence and avoid sub-optimum; (3) good adaptation to different wing vibration control systems by improving dynamic balance between the exploration and the exploitation; and (4) small computation cost due to small population size and small iteration number.


## 3.3. Procedures of Hybrid GNPSO Based Control Optimization

Figure2 shows control optimization procedures by proposed GNPSO hybrid algorithm. The GNPSO optimization process includes the following steps:Figure 2
(a) Optimization procedure by proposed GNPSO algorithm. (b) Procedures of calculating distance by new grey analysis rule. (c) Procedures of updating acceleration coefficients.
(a)
(b)
(c)Step (1) Initialization of IGPSO: Set the algorithm parameters, such as population size and iteration number of IGPSO and Initialize the population of particles.
Step (2) Evaluation of objective function:Evaluate the fitness of the particle in the population and find the fittest particle in the current generation. The evaluation is made by assigning the value of particle to control parameters, running the control system and calculating the objective function J in Equation (8). The fittest particle is obtained with minimized J.
Step (3) Calculate the distance between particle and fittest one: Calculate the grey relational grade g
i, which indicates the distance between the particle and the fittest particle. Based on the analysis of geometric similarity between the particle and fittest particle, g
i is obtained by new grey analysis rule using Equations (12) and (13).
Step (4) Update acceleration coefficients:Set scalars in Equation (16) to define the range and the sum of acceleration coefficients c
1
i
,
c
2
i. Then, those acceleration coefficients are adjusted by the new updating formula using Equations (14) and (15).
Step (5) Update velocity and position of particle: Based on the adjusted acceleration coefficients in step (4), update the velocity and the position of the particle.
Step (6) Output optimum of IGPSO search: If the maximum iteration number is not reached, go back to (2). Otherwise, output the best result obtained by IGPSO search in step (1)–step (5).
Step (7) Initialization of NLJ: Define the optimum of IGPSO search in step (6) as initial points of the NLJ optimization. The population of NLJ is initialized with initial search space.
Step (8) NLJ search: The search space is narrowed to improve the local search of the GNPSO. The new population is generated with narrowed search space. Then, the fittest one among the population is found with minimized objective function.
Step (9) Output global optimum: if the condition, such as iteration number is not reached, go back to Otherwise, output global optimum by NLJ search in GNPSO algorithm. The global optimum is the final optimized control parameter of IMC-PID controller for the wing vibration system.


## 4. Simulation Results

To demonstrate the effectiveness and the applicability of proposed GNPSO optimized IMC-PID method, simulation results for three benchmarks of different wing vibration systems, including linear system, nonlinear system and multiple-input-multiple-output (MIMO) system. For comparison, we choose the competitive algorithms, including recent Grey-PSO algorithm (GPSO) for PID design [24], recently published grey enhanced algorithm (i.e., grey adaptive DE algorithm (GADE) [27, 28]), reported NLJ hybrid algorithms for the IMC-PID design (i.e. NLJ-PSO hybrid algorithm (NPSO) [23], DE and NLJ hybrid algorithm (DE–NLJ) [22]) and classic evolutionary algorithm of GA [20] and PSO [21] for PID design. All simulations are implemented in Matlab 2014 environment using Intel Core i7-7500U CPU, 2.9 GHz, 8G RAM memory.Based on the consideration of the algorithm parameters in [20, 24, 27, 28] and fair comparison, the control parameters setting of GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for optimal IMC-PID control of wing vibration systems from Case 1 to Case 3 are given in Table 1. Further reduced population size and iteration number are given, according to the settings in GPSO optimized PID control [24]. The 25 independent runs are implemented for all algorithms.Table 1
Parameters setting of GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for optimal IMC-PID control of the linear wing vibration system.
Algorithm
Parameters setting
GA [20]
P
=
10
,
G
m
a
x
=
20, the crossover probability p
c
=
0.6, the mutation probability p
m
=
0.1
−
0.01
∗
n
/
P, where n
=
1
,
2
,
…
,
P
PSO [21]
P
=
10
,
G
m
a
x
=
20, the upper limit and the lower limit of inertia weight w
m
a
x
=
0.9
,
w
m
i
n
=
0.1, the upper limit and the lower limit of the velocity V
m
a
x
=
1
,
V
m
i
n
=
−
1, acceleration coefficients c
1
=
1
,
c
2
=
1
DE–NLJ [22]
P
=
10
,
G
m
a
x
=
20, the scaling factor F
=
1.2, the crossover rate C
R
=
0.3, the permit change rate η
=
1.2
NPSO [23]
P
=
10
,
G
m
a
x
=
20, w
m
a
x
=
0.9
,
w
m
i
n
=
0.1, V
m
a
x
=
1
,
V
m
i
n
=
−
1, acceleration coefficients c
1
=
1
,
c
2
=
1, the permit change rate η
=
1.2
GPSO [24]
P
=
10
,
G
m
a
x
=
20, w
m
a
x
=
0.9
,
w
m
i
n
=
0.1, V
m
a
x
=
1
,
V
m
i
n
=
−
1, grey relational grade g
i by Equation (6), the adjustable acceleration coefficients by Equation (7)
GADE [27]
P
=
10
,
G
m
a
x
=
20, the scaling factor F
=
2
g
i
F
0
−
F
f
+
2
F
0
−
F
f, the crossover rate C
R
=
2
g
i
CR
f
−
CR
0
+
2
CR
0
−
CR
f, grey relational grade g
i by Equation (6)
GNPSO
P
=
10
,
G
m
a
x
=
20, w
m
a
x
=
0.9
,
w
m
i
n
=
0.1, V
m
a
x
=
1
,
V
m
i
n
=
−
1, grey relational grade g
i by Equation (9), the adjustable acceleration coefficients by Equations (11)–(13), the permit change rate η
=
1.2


### 4.1. Case 1: Linear Wing Vibration System

For the linear wing vibration system with a single control surface, the dynamic equation in is specified based on Equation (1) as [2, 8]:(17)
M
=
m
t
m
w
x
θ
b
m
w
x
θ
b
I
θ
,
C
=
c
h
0
0
c
θ
,
K
=
k
h
0
0
k
θ
,
F
=
ρ
U
2
b
c
l
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
c
l
,
β
β
ρ
U
2
b
2
c
m
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
2
c
m
,
β
β
,
u
=
β
,wherem
w
,
m
t are mass coefficients, x
θ
b is the distance between the mass centre and the elastic axis, I
θ is the mass moment of inertia, c
h
,
c
θ are linear damping coefficients, k
h
,
k
θ are linear structural stiffness, ρ is the air density, U is the wind velocity, b is the semi-chord of the wing section, c
l
,
θ
,
c
m
,
θ
,
c
l
,
β
,
c
m
,
β are aerodynamic coefficients. The linear wing vibration system is controlled by trailing edge flap control surface. The system outputs are pitch angle θ and plunge displacement h. The input of the controlled system is flap angle β. The pitch angle is chosen as the sensor signal for the feedback control.Here, the search space of the filter parameter of the IMC-PID controller is in the range of (0, 10). Then, the corresponding PID parameters can be obtained using the Equation (6). The control goal is to regulate the plunge displacement and the pitch angle by single flap control surface.Table2 reveals the statistical performance of all algorithms, such as the minimum (F
m
i
n), the maximum (F
m
a
x), average (F
a
v
e
r
.), median (F
m
e
d
i
.) and standard deviation (F
s
d) values of the IAE based control performance J
b
e
s
t for 25 independent runs. Apparently, the GNPSO performs better than other algorithms with the best statistical results. Also, several different statistical tests are used to further analyze the performance of NPSO, GPSO, DE–NLJ, GADE and proposed GNPSO for the linear wing vibration system with single trailing edge flap. First, analysis of the variance (ANOVA) test of the control performance J
b
e
s
t for different algorithms are shown in Figure 3. Then, the results of Friedman, Friedman aligned and Quade tests are presented in Table 3 to show the ranks, the statistic value and associated probability p-value, where the GNPSO algorithm is best among all algorithms. Besides, in order to compare GNPSO with other six algorithms, Table 4 reveals the results of Wilcoxon signed ranks test for the control performance J
b
e
s
t according to different algorithms. In Table 4, R
+ represents the sum of ranks for the problems where the first algorithm is better than the second one, and R
− is the sum of ranks for the opposite. GNPSO shows a significant improvement over GA, PSO and GPSO, with a level of significance α
=
0.01, over DE–NLJ and NPSO with α
=
0.05, and over GADE with α
>
0.2. It is found that GA, DE–NLJ, and GADE have larger abnormal points and much worse F
worst than PSO based algorithms, although there are several good runs by GADE algorithm during all independent runs. In practice, it is important to avoid the bad control optimization results since the large vibration responses could be disaster to the wing vibration system. The PSO hybrid algorithms, such as NPSO, GPSO, and GNPSO can efficiently explore the optimal result with relatively stable performance, which ensures the final IMC-PID control performance of the linear wing vibration system.Table 2
Statistical results of the control performance using different algorithm for linear wing vibration system.
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
F
b
e
s
t
80.9772
80.8875
80.8869
80.8641
80.8757
80.8623
80.8623
F
w
o
r
s
t
547.4927
137.1411
99.5311
96.4319
241.7747
270.6472
83.3646
F
a
v
e
r
.
226.6412
87.6250
83.5841
86.1046
90.9597
98.2388
81.2842
F
m
e
d
i
.
219.0997
83.3390
81.6630
84.6301
81.3358
80.9670
80.9572
F
s
d
116.6897
11.4575
4.4233
4.6443
32.4401
48.7330
0.7137
The bold values mean the best values among all the algorithms.Table 3
Ranks, statistics and relatedp-values achieved by the Friedman, Friedman aligned, and Quade tests for control performance of linear wing vibration control system obtained by different algorithms.
Algorithms
Friedman
Friedman aligned
Quade
GA
6.4800
159.3200
6.8062
PSO
4.4400
79.0000
4.8615
DE–NLJ
3.7600
78.9600
3.4554
GADE
3.0400
81.4400
3.0492
NPSO
3.5200
70.2800
3.2277
GPSO
4.3600
80.8800
4.2031
GNPSO
2.4000
66.1200
2.3969
Statistics
54.87
38.7916
13.251
p-value
4.9145e−10
7.8635e−07
6.4356e−12Table 4
Wilcoxon signed ranks test results for the control performance of linear wing vibration control system obtained by different algorithms.
Algorithms
R
+
R
−
p-value
GNPSO versus GA
325
0
1.2290e-05
GNPSO versus PSO
325
0
1.2290e-05
GNPSO versus DE–NLJ
255
70
0.0128
GNPSO versus GADE
203
122
0.2758
GNPSO versus NPSO
274
51
0.0027
GNPSO versus GPSO
306
19
1.1286e-04Figure 3
ANOVA tests of the control performance obtained by different algorithms for linear wing vibration system.Figure4 illustrates the convergence behavior of the best run F
b
e
s
t in Table 2, obtained by all algorithms for linear wing vibration system with single control surface. It reveals that GNPSO reaches a better fitness than GA, PSO after 10 iterations. GNPSO has faster convergence speed than NPSO, GPSO, DE–NLJ and GADE. Thus, GNPSO has the best convergence characteristics among all algorithms. To demonstrate the effectiveness of those algorithms during all runs, Table 5 gives the optimal IMC-PID control parameters obtained by different algorithms according to the worst run F
w
o
r
s
t.Figure 4
Comparison of convergence process of the best fitnessF
m
i
n obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for linear wing vibration system.Table 5
Optimization results of IMC-PID controller parameters obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO associated withF
m
a
x for linear wing vibration system.
Control parameters
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
λ
1.6366
0.2811
0.1754
0.1165
0.5713
0.6723
0.1262
k
p
λ
0.2607
0.2516
0.2451
0.2365
0.2571
0.2579
0.2384
k
i
λ
−0.3203
−1.8508
−2.9498
−4.4083
−0.9150
−0.7780
−4.0764
k
d
λ
−0.2152
−0.0479
−0.0416
−0.0440
−0.0795
−0.0918
−0.0430Table6 compares control performance indices for each algorithm, involving the overall control performance, the ITAE index of system outputs and control effort defined as I
T
A
E
h
,
I
T
A
E
θ and I
T
A
E
β, the overshoot of system outputs defined as m
a
x
h and m
a
x
θ, the settling time of system outputs and control signal defined as t
s
_
h, t
s
_
θ and t
s
_
β. The system responses of the plunging, the pitching and the control surface associated with the F
w
o
r
s
t are given in Figure 5. The GNPSO performs best in terms of most indices. Although GA and NPSO have minimum value in certain index, GA has the largest settling time of t
s
_
h
,
t
s
_
θ
,
t
s
_
β. Besides, I
T
A
E
h
,
I
T
A
E
β
,
m
a
x
h
,
t
s
_
h, t
s
_
θ and t
s
_
β obtained by NPSO are worse than those by GNPSO. Hence, the GNPSO has better control performance, compared to GA, PSO, DE–NLJ, GADE, NPSO and GPSO for IMC-PID control of linear wing vibration system with single control surface.Table 6
Comparison of the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for linear wing vibration system.
Algorithms
J
b
e
s
t
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
β
m
a
x
h
m
a
x
θ
t
s
_
h
t
s
_
θ
t
s
_
β
GA
547.4927
6.6374
266.6325
279.099
0.0249
0
10.00
>10.00
>10.00
PSO
137.1411
2.1006
11.5850
10.4237
0.0335
0.0211
5.0000
6.9000
3.5200
NPSO
99.5311
1.4641
7.3532
4.8223
0.0281
0.0712
4.3800
4.7800
2.7000
GPSO
96.4319
1.3808
17.3114
5.3215
0.0281
0.1865
5.0600
5.6800
3.8200
DE–NLJ
241.7747
3.7608
36.8931
38.4890
0.0360
0.0014
6.2200
10.0000
5.6400
GADE
270.6472
3.6972
50.1241
52.3179
0.0358
0.0011
6.1800
10.0000
6.0600
GNPSO
83.3646
0.9969
9.9141
3.7646
0.0281
0.1345
3.9800
5.0800
2.5800
The bold values mean the best values among all the algorithms.Figure 5
Comparison of (a) plunging response, (b) pitch response and (c) control effort obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO based IMC-PID control of linear wing vibration system.
(a)
(b)
(c)


### 4.2. Case 2: Nonlinear Wing Vibration System

To further demonstrate the effectiveness of proposed GNPSO based IMC-PID method, the nonlinear wing vibration system with single control surface in reference [4, 10, 14] is tested. The nonlinear wing vibration system with nonlinear torsinal stiffness K
θ is represented based on Equation (1):(18)
K
θ
=
k
h
0
0
k
θ
0
+
k
θ
1
θ
+
k
θ
2
θ
2
+
k
θ
3
θ
3
+
k
θ
4
θ
4
+
⋯
,
M
=
m
t
m
w
x
θ
b
m
w
x
θ
b
I
θ
,
C
=
c
h
0
0
c
θ
,
F
=
ρ
U
2
b
c
l
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
c
l
,
β
β
ρ
U
2
b
2
c
m
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
2
c
m
,
β
β
,
u
=
β
,where the nonlinear function of the torsional stiffnessk
θ has coefficients of k
θ
1
,
k
θ
2
,
k
θ
3
,
k
θ
4. The nonlinearity in the system brings more difficulty in the design of the IMC-PID parameters. The system outputs are pitch angle θ and plunge displacement h. The input of the controlled system is flap angle β. The pitch angle is chosen as the sensor signal for the feedback control. The lower and upper limit of the filter parameter is taken as 0 and 10. The reference for system outputs are set as zero.Table7 shows the statistical results of all algorithms for 25 independent runs. For Case 2, GNPSO performs best on all performance indices of F
m
i
n
,
F
m
a
x
,
F
a
v
e
r
.
,
F
m
e
d
i
. and F
s
d. In addition, several statistical tests are carried out to compare the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system with single control surface. The ANOVA test results in Figure 6 demonstrate the best control performance of GNPSO. Table 8 reveals the ranks, the statistic value, and associated probability p-value, obtained by Friedman, Friedman aligned and Quade tests. Besides, Wilcoxon signed ranks test results are shown in Table 9, to compare the GNPSO with other six algorithms. Clearly, the GNPSO shows the improvement over GA, PSO, DE–NLJ, and GPSO, with α
=
0.01, and over GADE and NPSO with α
=
0.05 for optimal IMC-PID control of nonlinear wing vibration system with a single control surface.Table 7
Statistical results of the control performance using different algorithm for nonlinear wing vibration system.
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
F
b
e
s
t
104.9663
99.1740
99.1408
99.1501
99.1464
99.1403
99.1399
F
w
o
r
s
t
391.7290
107.4104
104.2291
104.3953
109.1813
345.3578
101.0503
F
a
v
e
r
.
172.9153
100.1079
100.3184
100.3567
101.3655
114.4408
99.5329
F
m
e
d
i
.
144.7146
99.5356
99.7550
99.7277
100.1793
99.4662
99.2658
F
s
d
75.5766
1.6712
1.4137
1.3830
2.9320
49.2240
0.4935
The bold values mean the best values among all the algorithms.Figure 6
ANOVA tests of the control performance obtained by different algorithms for nonlinear wing vibration system.Table 8
Ranks, statistics and relatedp-values achieved by the Friedman, Friedman aligned, and Quade tests for control performance of nonlinear wing vibration system obtained by different algorithms.
Algorithms
Friedman
Friedman aligned
Quade
GA
6.8400
157.6000
6.8677
PSO
3.6800
74.7200
4.2892
DE–NLJ
3.9600
75.2400
4.3077
GADE
3.7200
94.0400
3.4800
NPSO
3.5600
74.5200
3.0123
GPSO
3.6000
72.2800
3.7446
GNPSO
2.6400
67.6000
2.2985
Statistics
55.9886
38.6520
12.605
p-value
2.9265e−10
8.3744e−07
2.1316e−11Table 9
Wilcoxon signed ranks test results for the control performance of nonlinear wing vibration system obtained by different algorithms.
Algorithms
R
+
R
−
p-value
GNPSO versus GA
325
0
1.2290e-05
GNPSO versus PSO
320
5
2.2568e-05
GNPSO versus DE–NLJ
279
46
0.0017
GNPSO versus GADE
239
86
0.0396
GNPSO versus NPSO
240
85
0.0370
GNPSO versus GPSO
260
65
0.0087Figure7 shows the convergence behavior of the best run F
b
e
s
t obtained by each algorithm for nonlinear wing vibration system with single control surface. It is clear that GNPSO obtains a better fitness than GA, PSO, GPSO and DE–NLJ after 10 iterations and it also has faster convergence than NPSO and GADE, which indicates the better convergence characteristics of GNPSO than other algorithms. Table 10 presents the optimal control parameters obtained by each algorithm with respect to the F
worst for nonlinear wing vibration system with single control surface. The corresponding control performance indices are summarized in Table 11, involving the IAE control performance of system outputs and control signal, denoted by I
T
A
E
h
,
I
T
A
E
θ and I
T
A
E
β, the overshoot of system outputs, defined as m
a
x
h and m
a
x
θ and the settling time of system outputs and control signal, denoted as t
s
_
h
,
t
s
_
θ and t
s
_
β. It is obvious that more than half of control performance indices obtained by GNPSO are better than those by other six algorithms. It is struggling for GA, PSO, NPSO and DE–NLJ to handle the good control performance in different vibration motion as well as the reduction of control effort, although they have best result in terms of one performance index. Seven out of nine performance indices obtained by GPSO are worse than those by GNPSO. Besides, GA and GADE lead to relatively large IAE values and settling time among all algorithms. The results reveal that it is very difficult to obtain the optimum of control parameters on Case 2 since there is only one control surface to manipulate the different vibration motions and deal with system nonlinearity.Figure 7
Comparison of convergence process of theF
m
i
n obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system.Table 10
Optimization results of IMC-PID controller parameters obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO associated withF
m
a
x for nonlinear wing vibration system.
Control parameters
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
λ
1.5401
0.2374
0.2859
0.2909
0.3209
1.3378
0.2543
k
p
λ
0.7604
0.7782
0.7745
0.7742
0.7726
0.7608
0.7767
k
i
λ
−0.9855
−6.4746
−5.3626
−5.2693
−4.7713
−1.1350
−6.0383
k
d
λ
−0.5895
−0.1113
−0.1267
−0.1283
−0.1383
−0.5133
−0.1165Table 11
Comparison of the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system.
Algorithms
J
b
e
s
t
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
β
m
a
x
h
m
a
x
θ
t
s
_
h
t
s
_
θ
t
s
_
β
GA
391.7290
1.2609
78.5972
128.1772
0.0062
0.0239
>10.00
>10.00
>10.00
PSO
107.4104
0.6742
4.3397
5.7009
0.0077
0.0521
6.1200
4.9000
2.7600
NPSO
104.2291
0.2464
5.1856
7.1951
0.0059
0.0496
4.7800
4.8000
2.8600
GPSO
104.3953
0.2474
5.1528
7.2151
0.0058
0.0472
4.8000
4.8400
2.8800
DE–NLJ
109.1813
0.2813
5.3502
7.7510
0.0052
0.0381
5.2200
4.8600
2.6800
GADE
345.3578
1.0240
60.5494
98.1137
0.0063
0.0241
>10.00
>10.00
>10.00
GNPSO
101.0503
0.4216
4.2941
5.8438
0.0072
0.0377
4.7000
4.6600
2.5600
The bold values mean the best values among all the algorithms.The system responses and control effort are shown in Figure8. The steady-state error of the control system and the control effort obtained by GA and GADE are the worst. Compared to PSO, NPSO, GPSO, and DE–NLJ, GNPSO obtains improved overshoot of the pitching response, smaller settling time of system outputs and smoother control signal. With the difficulty of the control problem and multiple control requirements in Section 2.1, the control optimization could be easily trapped in local optimum and the search could be guided to meet certain control task rather than all tasks. Those reasons explain the sub-optimum convergence of GA, PSO, GPSO and DE–NLJ in Figure 7 and the one best performance index of GA, PSO, NPSO, and DE–NLJ in Table 11. Overall, GNPSO obtains the best control performance for IMC-PID control of nonlinear wing vibration system with single control surface.Figure 8
Comparison of (a) plunging response, (b) pitch response, and (c) control effort obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO based IMC-PID control of nonlinear wing vibration system.
(a)
(b)
(c)


### 4.3. Case 3: MIMO Nonlinear Wing Vibration System

In order to further demonstrate the effectiveness of GNPSO based method on the multiple input, multiple output (MIMO) system, this subsection represents the optimal IMC-PID control of nonlinear wing vibration system with two control surfaces using GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO. The nonlinear MIMO wing vibration system is a four degree of freedom system with two control surfaces as indicated in reference [11, 15, 49]. The dynamic equation of the nonlinear MIMO system based on Equation (1) is presented as:(19)
F
=
ρ
U
2
b
s
c
l
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
s
c
l
,
β
β
+
ρ
U
2
b
s
c
l
,
γ
γ
ρ
U
2
b
2
s
c
m
,
θ
−
e
f
f
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
2
s
c
m
,
β
−
e
f
f
β
+
ρ
U
2
b
2
s
c
m
,
γ
−
e
f
f
γ
,
u
=
β
γ
,
K
θ
=
k
h
0
0
k
θ
0
+
k
θ
1
θ
+
k
θ
2
θ
2
+
k
θ
3
θ
3
+
k
θ
4
θ
4
+
⋯
M
=
m
t
m
w
x
θ
b
m
w
x
θ
b
I
θ
,
C
=
c
h
0
0
c
θ
,where the vibration motions are controlled by both leading edge flap angleγ and trailing edge flap angle β. The effective aerodynamic coefficients are taken as c
m
,
θ
−
e
f
f
=
0.5
+
a
c
l
,
θ
+
2
c
m
,
θ
,
c
m
,
β
−
e
f
f
=
0.5
+
a
c
l
,
β
+
2
c
m
,
β
,
c
m
,
γ
−
e
f
f
=
0.5
+
a
c
l
,
γ
+
2
c
m
,
γ
. The system outputs are denoted as the plunge displacement h and the pitch angle θ. The control signals are taken as γ and β.For this MIMO system, two feedback SISO control loops are designed with one IMC-PID controller for each loop. One control loop is fromβ to θ and the other control loop is from γ to h. For each SISO control loop, a different process model is identified to design the corresponding IMC-PID controller. The search space of the filter parameter is in the range of (0, 10). The reference for system outputs are zero. To emphasize on the vibration suppression on Case 3, the objective function is modified as J
=
∫
0
t
final
h
/
h
m
a
x
+
θ
/
θ
m
a
x
d
t. Here, the control effort is not considered for the control optimization.Table12 reveals the statistical results of each algorithms in terms of performance indices F
b
e
s
t
,
F
worst
,
F
a
v
e
r
.
,
F
m
e
d
i
. and F
sd for 25 independent runs. Those statistical indices obtained by GNPSO are better than other algorithms in Table 12. Also, statistical tests are implemented to compare the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system with two control surface. Figure 9 shows the ANOVA test results, which indicates the best performance of GNPSO.Table 12
Statistical results of the control performance using different algorithm for MIMO nonlinear wing vibration system.
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
F
m
i
n
125.0941
114.6958
112.7629
112.6886
112.6937
112.6758
112.6719
F
m
a
x
412.5511
360.1039
156.2333
248.2655
153.5866
424.4457
131.0180
F
a
v
e
r
.
204.3738
194.7840
119.8541
151.7625
118.2837
156.3286
115.5874
F
m
e
d
i
.
179.8023
192.9319
115.1734
131.1979
114.5902
117.8508
114.3619
F
s
d
83.0342
64.5084
10.5619
47.2295
9.4717
88.8636
4.8418
The bold values mean the best values among all the algorithms.Figure 9
ANOVA tests of the control performance obtained by different algorithms for MIMO nonlinear wing vibration system.Table13 presents the ranks, the statistic value and associated probability p-value, obtained by Friedman, Friedman aligned, and Quade tests. Table 14 shows Wilcoxon signed ranks test results by comparing the GNPSO with other six algorithms. It is obvious that GNPSO shows the improvement over GA, PSO and GPSO with α
=
0.01, over GADE with α
=
0.05, over NPSO with α
=
0.2 and over DE–NLJ with α
>
0.2 for optimal IMC-PID control of nonlinear MIMO wing vibration system.Table 13
Ranks, statistics and relatedp-values achieved by the Friedman, Friedman aligned, and Quade tests for control performance of MIMO nonlinear wing vibration system obtained by different algorithms.
Algorithms
Friedman
Friedman aligned
Quade
GA
6.2000
135.8800
6.0800
PSO
5.4800
131.9600
5.8215
DE–NLJ
2.8400
57.0800
2.5446
GADE
3.7200
85.5200
4.0400
NPSO
2.9200
59.3200
2.7846
GPSO
4.6000
94.5200
4.7231
GNPSO
2.2400
51.7200
2.0062
Statistics
70.0629
39.0801
18.095
p-value
3.9693e−13
6.9034e−07
1.4433e−15Table 14
Wilcoxon signed ranks test results for the best control performance of MIMO nonlinear wing vibration system with two control surfaces obtained by different algorithms.
Algorithms
R
+
R
−
p-value
GNPSO versus GA
325
0
1.2290e-05
GNPSO versus PSO
318
7
2.8639e-05
GNPSO versus DE–NLJ
201
124
0.3002
GNPSO versus GADE
280
45
0.0016
GNPSO versus NPSO
222
103
0.1094
GNPSO versus GPSO
316
9
3.6243e-05Figure10 illustrates the convergence characteristics of the best run F
m
i
n obtained by each algorithm for nonlinear MIMO wing vibration system. Clearly, GNPSO reaches a better fitness with faster convergence than GA, PSO, GPSO, DE–NLJ, NPSO and GADE, which indicates the best convergence characteristics of GNPSO. Table 15 shows optimized control parameters by all algorithms with respect to the F
m
a
x for nonlinear MIMO wing vibration system. Table 16 gives the control performance indices, including ITAE values of I
T
A
E
h
,
I
T
A
E
θ, the overshoot of m
a
x
h and m
a
x
θ and settling time of t
s
_
h
,
t
s
_
θ. Obviously, seven out of nine control performance indices obtained by GNPSO are best among those obtained by all algorithms. GA, PSO and GPSO have undamped oscillations. DE–NLJ and GADE have large settling time over 10 seconds. Beside, eight out of nine control performance indices obtained by NPSO are worse than those by GNPSO. Figure 11 shows the system responses. Clearly, only GNPSO and NPSO have totally suppressed vibration responses in 10 seconds while other algorithms have limit cycle oscillations or even unstable oscillation in system responses. Compared to NPSO, GNPSO improves the oscillations in both plunging and pitching responses. Therefore, GNPSO obtains the better control performance than other six algorithms for nonlinear MIMO wing vibration system with two control surface.Figure 10
Comparison of convergence process of the best fitnessF
m
i
n obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for MIMO nonlinear wing vibration system.Table 15
Optimization results of IMC-PID controller parameters obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO associated withF
m
a
x for MIMO nonlinear wing vibration system.
Control parameters
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
λ
1
1.0539
0.1696
0.2379
0.2382
0.1348
2.1651
0.1820
λ
2
1.8535
3.0133
1.3843
2.6545
0.1616
2.2436
0.3495
k
p
1
λ
1
12.77
12.01
12.27
12.27
11.79
12.85
12.07
k
i
1
λ
1
−24.39
−147.00
−105.8
−105.75
−183.26
−11.91
−137.32
k
d
1
λ
1
−6.72
−1.13
−1.53
−1.53
−0.95
−13.89
−1.20
k
p
2
λ
2
−495.53
−497.33
−493.95
−496.94
−446.81
−496.34
−475.46
k
i
2
λ
2
539.53
331.92
722.31
376.77
6163.2
445.75
2856.5
k
d
2
λ
2
461.83
750.83
345.43
661.33
67.69
558.92
97.63Table 16
Comparison of the best control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for MIMO nonlinear wing vibration system.
Algorithms
F
b
e
s
t
I
T
A
E
h
I
T
A
E
θ
I
A
E
h
I
A
E
θ
m
a
x
h
m
a
x
θ
t
s
_
h
t
s
_
θ
GA
125.0941
3.6018
123.3021
0.7933
25.6712
0.0006
0.0499
Unstable
Unstable
PSO
114.6958
5.0979
64.7627
0.9644
13.9298
0.0035
0.0406
Oscillation
Oscillation
NPSO
112.7629
0.6801
22.1730
0.3285
8.8805
0.0013
0.0113
6.5200
6.6000
GPSO
112.6886
2.6314
41.3592
0.6063
11.5521
0.0015
0.0204
oscillation
oscillation
DE–NLJ
112.6937
1.4080
19.2991
0.3802
6.8511
0.0012
0.0085
>10.00
>10.00
GADE
112.6758
4.1785
136.7696
0.8526
26.4504
0.0004
0
>10.00
>10.00
GNPSO
112.6719
0.3548
17.5068
0.2082
7.3019
0.0003
0
7.320
6.380
The bold values mean the best values among all the algorithms.Figure 11
Comparison of (a) plunging response, (b) pitch response obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO based IMC-PID control of MIMO nonlinear wing vibration system.
(a)
(b)


### 4.4. Robustness Test

In order to further demonstrate the robustness of GNPSO, the comparison of system responses under varied system parameters obtained by GA, PSO, GPSO, NPSO, DE–NLJ GADE and GNPSO for nonlinear MIMO wing vibration system are presented in this subsection. The variations of system parameters are given in Table17. The robust performance of all algorithms is evaluated by the indices of I
T
A
E
h and I
T
A
E
θ. Table 17 presents the robust performance of nonlinear MIMO wing vibration control system under varied system parameters, obtained by different algorithms. Clearly, GNPSO performs best among seven algorithms since the values of I
T
A
E
h and I
T
A
E
θ vary in relatively small range according to variations of system parameters. For GA, PSO, GPSO and GADE, they cannot maintain good robustness with significant change of I
T
A
E
θ from case to case. DE–NLJ and NPSO have doubled the value of I
T
A
E
h when it comes to the parameter change of c
l
,
β
,
c
m
,
β. Hence, GNPSO obtains stronger robustness than other six algorithms due to the smaller variations in performance indices under varied system parameters.Table 17
Robust performance of GA, PSO, NPSO, GPSO, DE–NLJ, GADE, and GNPSO for MIMO nonlinear wing vibration control system.
Parameter variation
U
C
h
,
C
θ
c
l
,
θ
,
c
m
,
θ
c
l
,
β
,
c
m
,
β
% change
−10%
+10%
−10%
+10%
Algorithm
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
h
I
T
A
E
θ
GA
3.3648
113.0802
3.4791
121.0771
3.5214
118.4546
8.0568
242.6463
PSO
2.2940
25.4754
2.4915
28.4621
3.5588
41.9787
6.3555
89.2480
NPSO
0.4015
20.1432
0.5364
21.4972
0.4227
19.6518
0.9005
24.0710
GPSO
0.9014
19.0715
2.3205
37.1065
2.0433
30.8812
2.9261
48.9147
DE–NLJ
0.3434
17.7616
0.5704
19.8848
0.4441
18.3618
1.0579
20.2744
GADE
3.7324
123.2064
4.1695
142.4313
4.4529
144.3820
15.5231
463.8153
GNPSO
0.2510
15.2886
0.2846
17.5032
0.2519
16.2986
0.3190
18.6249


### 4.5. Computational Efficiency

In order to further analyze the computational efficiency of GNPSO, the statistical results of the computation time obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO on 25 independent runs are summarized in Table18, Table 19 and Table 20 for Case 1, Case 2, and Case 3 respectively. The statistical analysis of computation time for each case is implemented based on performance indices, including the minimum (t
m
i
n), the maximum (t
m
a
x), average (t
a
v
e
r
.), median (t
m
e
d
i
.) and standard deviation (t
sd). In practice, the fast tuning of the IMC-PID controller is favored to avoid the large time delay of the control signal and the uncontrolled system in possible damage. For the previous study of GPSO based PID control issue [24], the computation complexity is nearly O(PGn) with the population size of P, max. iteration number of G and the number of control parameters n, while those settings could cost large computation time of more than 200 seconds for complex wing vibration control systems. In this paper, to ensure the fast tuning of the IMC-PID controller, the computational cost is approximated as O(0.3
P
0.4
G
0.3
n) with reduced population size/iteration number and less control parameter. Tables 18–20 reveal that computation time obtained by the algorithms for complex wing vibration control systems are relatively small, due to reduced computation complexity. Clearly, among all algorithms, GNPSO performs best in terms of all performance indices for Cases 2 and 3, and three out of five indices for Case 1. Thus, GNPSO has better computational efficiency than other six algorithms on diverse wing vibration systems.Table 18
Statistical results of computation time (sec) using different algorithms for Case 1: linear wing vibration system with single control surface.
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
t
m
i
n
38.0083
36.4092
36.9464
41.9241
46.0005
35.4421
t
m
a
x
45.8949
38.5420
46.0858
49.6778
51.8526
46.4936
t
a
v
e
r
.
43.0086
37.4733
37.9925
44.1037
47.8175
36.6189
t
m
e
d
i
.
43.3580
37.4422
37.6460
43.6939
47.2339
36.3542
t
s
d
1.8242
0.4833
1.7357
1.7207
1.8155
2.1087
The bold values mean the best values among all the algorithms.Table 19
Statistical results of computation time (sec) using different algorithms for Case 2: nonlinear wing vibration system with single control surface.
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
t
m
i
n
41.5313
36.8535
36.8671
44.4655
53.3908
35.4934
t
m
a
x
44.2525
39.8687
40.1606
48.1763
60.1680
38.3784
t
a
v
e
r
.
42.4075
37.9693
38.4001
46.1848
55.0292
37.2044
t
m
e
d
i
.
42.3507
37.7893
38.4244
46.1072
54.5242
37.2493
t
s
d
0.7156
0.7506
0.7647
1.1402
1.5666
0.6607
The bold values mean the best values among all the algorithms.Table 20
Statistical results of computation time (sec) using different algorithms for Case 3: MIMO nonlinear wing vibration system with two control surfaces.
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
t
m
i
n
51.9130
46.5234
48.7196
55.4146
59.2800
45.8710
t
m
a
x
73.7114
61.4269
62.4517
64.0192
83.6449
51.5678
t
a
v
e
r
.
62.6266
52.2032
51.5982
58.5969
64.6545
47.6273
t
m
e
d
i
.
61.2508
49.2194
50.9832
58.2329
62.6389
47.3272
t
s
d
6.0806
5.0044
2.7682
1.9535
5.8313
1.0530
The bold values mean the best values among all the algorithms.


## 4.1. Case 1: Linear Wing Vibration System

For the linear wing vibration system with a single control surface, the dynamic equation in is specified based on Equation (1) as [2, 8]:(17)
M
=
m
t
m
w
x
θ
b
m
w
x
θ
b
I
θ
,
C
=
c
h
0
0
c
θ
,
K
=
k
h
0
0
k
θ
,
F
=
ρ
U
2
b
c
l
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
c
l
,
β
β
ρ
U
2
b
2
c
m
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
2
c
m
,
β
β
,
u
=
β
,wherem
w
,
m
t are mass coefficients, x
θ
b is the distance between the mass centre and the elastic axis, I
θ is the mass moment of inertia, c
h
,
c
θ are linear damping coefficients, k
h
,
k
θ are linear structural stiffness, ρ is the air density, U is the wind velocity, b is the semi-chord of the wing section, c
l
,
θ
,
c
m
,
θ
,
c
l
,
β
,
c
m
,
β are aerodynamic coefficients. The linear wing vibration system is controlled by trailing edge flap control surface. The system outputs are pitch angle θ and plunge displacement h. The input of the controlled system is flap angle β. The pitch angle is chosen as the sensor signal for the feedback control.Here, the search space of the filter parameter of the IMC-PID controller is in the range of (0, 10). Then, the corresponding PID parameters can be obtained using the Equation (6). The control goal is to regulate the plunge displacement and the pitch angle by single flap control surface.Table2 reveals the statistical performance of all algorithms, such as the minimum (F
m
i
n), the maximum (F
m
a
x), average (F
a
v
e
r
.), median (F
m
e
d
i
.) and standard deviation (F
s
d) values of the IAE based control performance J
b
e
s
t for 25 independent runs. Apparently, the GNPSO performs better than other algorithms with the best statistical results. Also, several different statistical tests are used to further analyze the performance of NPSO, GPSO, DE–NLJ, GADE and proposed GNPSO for the linear wing vibration system with single trailing edge flap. First, analysis of the variance (ANOVA) test of the control performance J
b
e
s
t for different algorithms are shown in Figure 3. Then, the results of Friedman, Friedman aligned and Quade tests are presented in Table 3 to show the ranks, the statistic value and associated probability p-value, where the GNPSO algorithm is best among all algorithms. Besides, in order to compare GNPSO with other six algorithms, Table 4 reveals the results of Wilcoxon signed ranks test for the control performance J
b
e
s
t according to different algorithms. In Table 4, R
+ represents the sum of ranks for the problems where the first algorithm is better than the second one, and R
− is the sum of ranks for the opposite. GNPSO shows a significant improvement over GA, PSO and GPSO, with a level of significance α
=
0.01, over DE–NLJ and NPSO with α
=
0.05, and over GADE with α
>
0.2. It is found that GA, DE–NLJ, and GADE have larger abnormal points and much worse F
worst than PSO based algorithms, although there are several good runs by GADE algorithm during all independent runs. In practice, it is important to avoid the bad control optimization results since the large vibration responses could be disaster to the wing vibration system. The PSO hybrid algorithms, such as NPSO, GPSO, and GNPSO can efficiently explore the optimal result with relatively stable performance, which ensures the final IMC-PID control performance of the linear wing vibration system.Table 2
Statistical results of the control performance using different algorithm for linear wing vibration system.
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
F
b
e
s
t
80.9772
80.8875
80.8869
80.8641
80.8757
80.8623
80.8623
F
w
o
r
s
t
547.4927
137.1411
99.5311
96.4319
241.7747
270.6472
83.3646
F
a
v
e
r
.
226.6412
87.6250
83.5841
86.1046
90.9597
98.2388
81.2842
F
m
e
d
i
.
219.0997
83.3390
81.6630
84.6301
81.3358
80.9670
80.9572
F
s
d
116.6897
11.4575
4.4233
4.6443
32.4401
48.7330
0.7137
The bold values mean the best values among all the algorithms.Table 3
Ranks, statistics and relatedp-values achieved by the Friedman, Friedman aligned, and Quade tests for control performance of linear wing vibration control system obtained by different algorithms.
Algorithms
Friedman
Friedman aligned
Quade
GA
6.4800
159.3200
6.8062
PSO
4.4400
79.0000
4.8615
DE–NLJ
3.7600
78.9600
3.4554
GADE
3.0400
81.4400
3.0492
NPSO
3.5200
70.2800
3.2277
GPSO
4.3600
80.8800
4.2031
GNPSO
2.4000
66.1200
2.3969
Statistics
54.87
38.7916
13.251
p-value
4.9145e−10
7.8635e−07
6.4356e−12Table 4
Wilcoxon signed ranks test results for the control performance of linear wing vibration control system obtained by different algorithms.
Algorithms
R
+
R
−
p-value
GNPSO versus GA
325
0
1.2290e-05
GNPSO versus PSO
325
0
1.2290e-05
GNPSO versus DE–NLJ
255
70
0.0128
GNPSO versus GADE
203
122
0.2758
GNPSO versus NPSO
274
51
0.0027
GNPSO versus GPSO
306
19
1.1286e-04Figure 3
ANOVA tests of the control performance obtained by different algorithms for linear wing vibration system.Figure4 illustrates the convergence behavior of the best run F
b
e
s
t in Table 2, obtained by all algorithms for linear wing vibration system with single control surface. It reveals that GNPSO reaches a better fitness than GA, PSO after 10 iterations. GNPSO has faster convergence speed than NPSO, GPSO, DE–NLJ and GADE. Thus, GNPSO has the best convergence characteristics among all algorithms. To demonstrate the effectiveness of those algorithms during all runs, Table 5 gives the optimal IMC-PID control parameters obtained by different algorithms according to the worst run F
w
o
r
s
t.Figure 4
Comparison of convergence process of the best fitnessF
m
i
n obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for linear wing vibration system.Table 5
Optimization results of IMC-PID controller parameters obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO associated withF
m
a
x for linear wing vibration system.
Control parameters
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
λ
1.6366
0.2811
0.1754
0.1165
0.5713
0.6723
0.1262
k
p
λ
0.2607
0.2516
0.2451
0.2365
0.2571
0.2579
0.2384
k
i
λ
−0.3203
−1.8508
−2.9498
−4.4083
−0.9150
−0.7780
−4.0764
k
d
λ
−0.2152
−0.0479
−0.0416
−0.0440
−0.0795
−0.0918
−0.0430Table6 compares control performance indices for each algorithm, involving the overall control performance, the ITAE index of system outputs and control effort defined as I
T
A
E
h
,
I
T
A
E
θ and I
T
A
E
β, the overshoot of system outputs defined as m
a
x
h and m
a
x
θ, the settling time of system outputs and control signal defined as t
s
_
h, t
s
_
θ and t
s
_
β. The system responses of the plunging, the pitching and the control surface associated with the F
w
o
r
s
t are given in Figure 5. The GNPSO performs best in terms of most indices. Although GA and NPSO have minimum value in certain index, GA has the largest settling time of t
s
_
h
,
t
s
_
θ
,
t
s
_
β. Besides, I
T
A
E
h
,
I
T
A
E
β
,
m
a
x
h
,
t
s
_
h, t
s
_
θ and t
s
_
β obtained by NPSO are worse than those by GNPSO. Hence, the GNPSO has better control performance, compared to GA, PSO, DE–NLJ, GADE, NPSO and GPSO for IMC-PID control of linear wing vibration system with single control surface.Table 6
Comparison of the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for linear wing vibration system.
Algorithms
J
b
e
s
t
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
β
m
a
x
h
m
a
x
θ
t
s
_
h
t
s
_
θ
t
s
_
β
GA
547.4927
6.6374
266.6325
279.099
0.0249
0
10.00
>10.00
>10.00
PSO
137.1411
2.1006
11.5850
10.4237
0.0335
0.0211
5.0000
6.9000
3.5200
NPSO
99.5311
1.4641
7.3532
4.8223
0.0281
0.0712
4.3800
4.7800
2.7000
GPSO
96.4319
1.3808
17.3114
5.3215
0.0281
0.1865
5.0600
5.6800
3.8200
DE–NLJ
241.7747
3.7608
36.8931
38.4890
0.0360
0.0014
6.2200
10.0000
5.6400
GADE
270.6472
3.6972
50.1241
52.3179
0.0358
0.0011
6.1800
10.0000
6.0600
GNPSO
83.3646
0.9969
9.9141
3.7646
0.0281
0.1345
3.9800
5.0800
2.5800
The bold values mean the best values among all the algorithms.Figure 5
Comparison of (a) plunging response, (b) pitch response and (c) control effort obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO based IMC-PID control of linear wing vibration system.
(a)
(b)
(c)


## 4.2. Case 2: Nonlinear Wing Vibration System

To further demonstrate the effectiveness of proposed GNPSO based IMC-PID method, the nonlinear wing vibration system with single control surface in reference [4, 10, 14] is tested. The nonlinear wing vibration system with nonlinear torsinal stiffness K
θ is represented based on Equation (1):(18)
K
θ
=
k
h
0
0
k
θ
0
+
k
θ
1
θ
+
k
θ
2
θ
2
+
k
θ
3
θ
3
+
k
θ
4
θ
4
+
⋯
,
M
=
m
t
m
w
x
θ
b
m
w
x
θ
b
I
θ
,
C
=
c
h
0
0
c
θ
,
F
=
ρ
U
2
b
c
l
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
c
l
,
β
β
ρ
U
2
b
2
c
m
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
2
c
m
,
β
β
,
u
=
β
,where the nonlinear function of the torsional stiffnessk
θ has coefficients of k
θ
1
,
k
θ
2
,
k
θ
3
,
k
θ
4. The nonlinearity in the system brings more difficulty in the design of the IMC-PID parameters. The system outputs are pitch angle θ and plunge displacement h. The input of the controlled system is flap angle β. The pitch angle is chosen as the sensor signal for the feedback control. The lower and upper limit of the filter parameter is taken as 0 and 10. The reference for system outputs are set as zero.Table7 shows the statistical results of all algorithms for 25 independent runs. For Case 2, GNPSO performs best on all performance indices of F
m
i
n
,
F
m
a
x
,
F
a
v
e
r
.
,
F
m
e
d
i
. and F
s
d. In addition, several statistical tests are carried out to compare the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system with single control surface. The ANOVA test results in Figure 6 demonstrate the best control performance of GNPSO. Table 8 reveals the ranks, the statistic value, and associated probability p-value, obtained by Friedman, Friedman aligned and Quade tests. Besides, Wilcoxon signed ranks test results are shown in Table 9, to compare the GNPSO with other six algorithms. Clearly, the GNPSO shows the improvement over GA, PSO, DE–NLJ, and GPSO, with α
=
0.01, and over GADE and NPSO with α
=
0.05 for optimal IMC-PID control of nonlinear wing vibration system with a single control surface.Table 7
Statistical results of the control performance using different algorithm for nonlinear wing vibration system.
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
F
b
e
s
t
104.9663
99.1740
99.1408
99.1501
99.1464
99.1403
99.1399
F
w
o
r
s
t
391.7290
107.4104
104.2291
104.3953
109.1813
345.3578
101.0503
F
a
v
e
r
.
172.9153
100.1079
100.3184
100.3567
101.3655
114.4408
99.5329
F
m
e
d
i
.
144.7146
99.5356
99.7550
99.7277
100.1793
99.4662
99.2658
F
s
d
75.5766
1.6712
1.4137
1.3830
2.9320
49.2240
0.4935
The bold values mean the best values among all the algorithms.Figure 6
ANOVA tests of the control performance obtained by different algorithms for nonlinear wing vibration system.Table 8
Ranks, statistics and relatedp-values achieved by the Friedman, Friedman aligned, and Quade tests for control performance of nonlinear wing vibration system obtained by different algorithms.
Algorithms
Friedman
Friedman aligned
Quade
GA
6.8400
157.6000
6.8677
PSO
3.6800
74.7200
4.2892
DE–NLJ
3.9600
75.2400
4.3077
GADE
3.7200
94.0400
3.4800
NPSO
3.5600
74.5200
3.0123
GPSO
3.6000
72.2800
3.7446
GNPSO
2.6400
67.6000
2.2985
Statistics
55.9886
38.6520
12.605
p-value
2.9265e−10
8.3744e−07
2.1316e−11Table 9
Wilcoxon signed ranks test results for the control performance of nonlinear wing vibration system obtained by different algorithms.
Algorithms
R
+
R
−
p-value
GNPSO versus GA
325
0
1.2290e-05
GNPSO versus PSO
320
5
2.2568e-05
GNPSO versus DE–NLJ
279
46
0.0017
GNPSO versus GADE
239
86
0.0396
GNPSO versus NPSO
240
85
0.0370
GNPSO versus GPSO
260
65
0.0087Figure7 shows the convergence behavior of the best run F
b
e
s
t obtained by each algorithm for nonlinear wing vibration system with single control surface. It is clear that GNPSO obtains a better fitness than GA, PSO, GPSO and DE–NLJ after 10 iterations and it also has faster convergence than NPSO and GADE, which indicates the better convergence characteristics of GNPSO than other algorithms. Table 10 presents the optimal control parameters obtained by each algorithm with respect to the F
worst for nonlinear wing vibration system with single control surface. The corresponding control performance indices are summarized in Table 11, involving the IAE control performance of system outputs and control signal, denoted by I
T
A
E
h
,
I
T
A
E
θ and I
T
A
E
β, the overshoot of system outputs, defined as m
a
x
h and m
a
x
θ and the settling time of system outputs and control signal, denoted as t
s
_
h
,
t
s
_
θ and t
s
_
β. It is obvious that more than half of control performance indices obtained by GNPSO are better than those by other six algorithms. It is struggling for GA, PSO, NPSO and DE–NLJ to handle the good control performance in different vibration motion as well as the reduction of control effort, although they have best result in terms of one performance index. Seven out of nine performance indices obtained by GPSO are worse than those by GNPSO. Besides, GA and GADE lead to relatively large IAE values and settling time among all algorithms. The results reveal that it is very difficult to obtain the optimum of control parameters on Case 2 since there is only one control surface to manipulate the different vibration motions and deal with system nonlinearity.Figure 7
Comparison of convergence process of theF
m
i
n obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system.Table 10
Optimization results of IMC-PID controller parameters obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO associated withF
m
a
x for nonlinear wing vibration system.
Control parameters
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
λ
1.5401
0.2374
0.2859
0.2909
0.3209
1.3378
0.2543
k
p
λ
0.7604
0.7782
0.7745
0.7742
0.7726
0.7608
0.7767
k
i
λ
−0.9855
−6.4746
−5.3626
−5.2693
−4.7713
−1.1350
−6.0383
k
d
λ
−0.5895
−0.1113
−0.1267
−0.1283
−0.1383
−0.5133
−0.1165Table 11
Comparison of the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system.
Algorithms
J
b
e
s
t
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
β
m
a
x
h
m
a
x
θ
t
s
_
h
t
s
_
θ
t
s
_
β
GA
391.7290
1.2609
78.5972
128.1772
0.0062
0.0239
>10.00
>10.00
>10.00
PSO
107.4104
0.6742
4.3397
5.7009
0.0077
0.0521
6.1200
4.9000
2.7600
NPSO
104.2291
0.2464
5.1856
7.1951
0.0059
0.0496
4.7800
4.8000
2.8600
GPSO
104.3953
0.2474
5.1528
7.2151
0.0058
0.0472
4.8000
4.8400
2.8800
DE–NLJ
109.1813
0.2813
5.3502
7.7510
0.0052
0.0381
5.2200
4.8600
2.6800
GADE
345.3578
1.0240
60.5494
98.1137
0.0063
0.0241
>10.00
>10.00
>10.00
GNPSO
101.0503
0.4216
4.2941
5.8438
0.0072
0.0377
4.7000
4.6600
2.5600
The bold values mean the best values among all the algorithms.The system responses and control effort are shown in Figure8. The steady-state error of the control system and the control effort obtained by GA and GADE are the worst. Compared to PSO, NPSO, GPSO, and DE–NLJ, GNPSO obtains improved overshoot of the pitching response, smaller settling time of system outputs and smoother control signal. With the difficulty of the control problem and multiple control requirements in Section 2.1, the control optimization could be easily trapped in local optimum and the search could be guided to meet certain control task rather than all tasks. Those reasons explain the sub-optimum convergence of GA, PSO, GPSO and DE–NLJ in Figure 7 and the one best performance index of GA, PSO, NPSO, and DE–NLJ in Table 11. Overall, GNPSO obtains the best control performance for IMC-PID control of nonlinear wing vibration system with single control surface.Figure 8
Comparison of (a) plunging response, (b) pitch response, and (c) control effort obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO based IMC-PID control of nonlinear wing vibration system.
(a)
(b)
(c)


## 4.3. Case 3: MIMO Nonlinear Wing Vibration System

In order to further demonstrate the effectiveness of GNPSO based method on the multiple input, multiple output (MIMO) system, this subsection represents the optimal IMC-PID control of nonlinear wing vibration system with two control surfaces using GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO. The nonlinear MIMO wing vibration system is a four degree of freedom system with two control surfaces as indicated in reference [11, 15, 49]. The dynamic equation of the nonlinear MIMO system based on Equation (1) is presented as:(19)
F
=
ρ
U
2
b
s
c
l
,
θ
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
s
c
l
,
β
β
+
ρ
U
2
b
s
c
l
,
γ
γ
ρ
U
2
b
2
s
c
m
,
θ
−
e
f
f
θ
+
h
˙
U
+
0.5
−
a
b
θ
˙
U
+
ρ
U
2
b
2
s
c
m
,
β
−
e
f
f
β
+
ρ
U
2
b
2
s
c
m
,
γ
−
e
f
f
γ
,
u
=
β
γ
,
K
θ
=
k
h
0
0
k
θ
0
+
k
θ
1
θ
+
k
θ
2
θ
2
+
k
θ
3
θ
3
+
k
θ
4
θ
4
+
⋯
M
=
m
t
m
w
x
θ
b
m
w
x
θ
b
I
θ
,
C
=
c
h
0
0
c
θ
,where the vibration motions are controlled by both leading edge flap angleγ and trailing edge flap angle β. The effective aerodynamic coefficients are taken as c
m
,
θ
−
e
f
f
=
0.5
+
a
c
l
,
θ
+
2
c
m
,
θ
,
c
m
,
β
−
e
f
f
=
0.5
+
a
c
l
,
β
+
2
c
m
,
β
,
c
m
,
γ
−
e
f
f
=
0.5
+
a
c
l
,
γ
+
2
c
m
,
γ
. The system outputs are denoted as the plunge displacement h and the pitch angle θ. The control signals are taken as γ and β.For this MIMO system, two feedback SISO control loops are designed with one IMC-PID controller for each loop. One control loop is fromβ to θ and the other control loop is from γ to h. For each SISO control loop, a different process model is identified to design the corresponding IMC-PID controller. The search space of the filter parameter is in the range of (0, 10). The reference for system outputs are zero. To emphasize on the vibration suppression on Case 3, the objective function is modified as J
=
∫
0
t
final
h
/
h
m
a
x
+
θ
/
θ
m
a
x
d
t. Here, the control effort is not considered for the control optimization.Table12 reveals the statistical results of each algorithms in terms of performance indices F
b
e
s
t
,
F
worst
,
F
a
v
e
r
.
,
F
m
e
d
i
. and F
sd for 25 independent runs. Those statistical indices obtained by GNPSO are better than other algorithms in Table 12. Also, statistical tests are implemented to compare the control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for nonlinear wing vibration system with two control surface. Figure 9 shows the ANOVA test results, which indicates the best performance of GNPSO.Table 12
Statistical results of the control performance using different algorithm for MIMO nonlinear wing vibration system.
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
F
m
i
n
125.0941
114.6958
112.7629
112.6886
112.6937
112.6758
112.6719
F
m
a
x
412.5511
360.1039
156.2333
248.2655
153.5866
424.4457
131.0180
F
a
v
e
r
.
204.3738
194.7840
119.8541
151.7625
118.2837
156.3286
115.5874
F
m
e
d
i
.
179.8023
192.9319
115.1734
131.1979
114.5902
117.8508
114.3619
F
s
d
83.0342
64.5084
10.5619
47.2295
9.4717
88.8636
4.8418
The bold values mean the best values among all the algorithms.Figure 9
ANOVA tests of the control performance obtained by different algorithms for MIMO nonlinear wing vibration system.Table13 presents the ranks, the statistic value and associated probability p-value, obtained by Friedman, Friedman aligned, and Quade tests. Table 14 shows Wilcoxon signed ranks test results by comparing the GNPSO with other six algorithms. It is obvious that GNPSO shows the improvement over GA, PSO and GPSO with α
=
0.01, over GADE with α
=
0.05, over NPSO with α
=
0.2 and over DE–NLJ with α
>
0.2 for optimal IMC-PID control of nonlinear MIMO wing vibration system.Table 13
Ranks, statistics and relatedp-values achieved by the Friedman, Friedman aligned, and Quade tests for control performance of MIMO nonlinear wing vibration system obtained by different algorithms.
Algorithms
Friedman
Friedman aligned
Quade
GA
6.2000
135.8800
6.0800
PSO
5.4800
131.9600
5.8215
DE–NLJ
2.8400
57.0800
2.5446
GADE
3.7200
85.5200
4.0400
NPSO
2.9200
59.3200
2.7846
GPSO
4.6000
94.5200
4.7231
GNPSO
2.2400
51.7200
2.0062
Statistics
70.0629
39.0801
18.095
p-value
3.9693e−13
6.9034e−07
1.4433e−15Table 14
Wilcoxon signed ranks test results for the best control performance of MIMO nonlinear wing vibration system with two control surfaces obtained by different algorithms.
Algorithms
R
+
R
−
p-value
GNPSO versus GA
325
0
1.2290e-05
GNPSO versus PSO
318
7
2.8639e-05
GNPSO versus DE–NLJ
201
124
0.3002
GNPSO versus GADE
280
45
0.0016
GNPSO versus NPSO
222
103
0.1094
GNPSO versus GPSO
316
9
3.6243e-05Figure10 illustrates the convergence characteristics of the best run F
m
i
n obtained by each algorithm for nonlinear MIMO wing vibration system. Clearly, GNPSO reaches a better fitness with faster convergence than GA, PSO, GPSO, DE–NLJ, NPSO and GADE, which indicates the best convergence characteristics of GNPSO. Table 15 shows optimized control parameters by all algorithms with respect to the F
m
a
x for nonlinear MIMO wing vibration system. Table 16 gives the control performance indices, including ITAE values of I
T
A
E
h
,
I
T
A
E
θ, the overshoot of m
a
x
h and m
a
x
θ and settling time of t
s
_
h
,
t
s
_
θ. Obviously, seven out of nine control performance indices obtained by GNPSO are best among those obtained by all algorithms. GA, PSO and GPSO have undamped oscillations. DE–NLJ and GADE have large settling time over 10 seconds. Beside, eight out of nine control performance indices obtained by NPSO are worse than those by GNPSO. Figure 11 shows the system responses. Clearly, only GNPSO and NPSO have totally suppressed vibration responses in 10 seconds while other algorithms have limit cycle oscillations or even unstable oscillation in system responses. Compared to NPSO, GNPSO improves the oscillations in both plunging and pitching responses. Therefore, GNPSO obtains the better control performance than other six algorithms for nonlinear MIMO wing vibration system with two control surface.Figure 10
Comparison of convergence process of the best fitnessF
m
i
n obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for MIMO nonlinear wing vibration system.Table 15
Optimization results of IMC-PID controller parameters obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO associated withF
m
a
x for MIMO nonlinear wing vibration system.
Control parameters
GA
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
λ
1
1.0539
0.1696
0.2379
0.2382
0.1348
2.1651
0.1820
λ
2
1.8535
3.0133
1.3843
2.6545
0.1616
2.2436
0.3495
k
p
1
λ
1
12.77
12.01
12.27
12.27
11.79
12.85
12.07
k
i
1
λ
1
−24.39
−147.00
−105.8
−105.75
−183.26
−11.91
−137.32
k
d
1
λ
1
−6.72
−1.13
−1.53
−1.53
−0.95
−13.89
−1.20
k
p
2
λ
2
−495.53
−497.33
−493.95
−496.94
−446.81
−496.34
−475.46
k
i
2
λ
2
539.53
331.92
722.31
376.77
6163.2
445.75
2856.5
k
d
2
λ
2
461.83
750.83
345.43
661.33
67.69
558.92
97.63Table 16
Comparison of the best control performance obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO for MIMO nonlinear wing vibration system.
Algorithms
F
b
e
s
t
I
T
A
E
h
I
T
A
E
θ
I
A
E
h
I
A
E
θ
m
a
x
h
m
a
x
θ
t
s
_
h
t
s
_
θ
GA
125.0941
3.6018
123.3021
0.7933
25.6712
0.0006
0.0499
Unstable
Unstable
PSO
114.6958
5.0979
64.7627
0.9644
13.9298
0.0035
0.0406
Oscillation
Oscillation
NPSO
112.7629
0.6801
22.1730
0.3285
8.8805
0.0013
0.0113
6.5200
6.6000
GPSO
112.6886
2.6314
41.3592
0.6063
11.5521
0.0015
0.0204
oscillation
oscillation
DE–NLJ
112.6937
1.4080
19.2991
0.3802
6.8511
0.0012
0.0085
>10.00
>10.00
GADE
112.6758
4.1785
136.7696
0.8526
26.4504
0.0004
0
>10.00
>10.00
GNPSO
112.6719
0.3548
17.5068
0.2082
7.3019
0.0003
0
7.320
6.380
The bold values mean the best values among all the algorithms.Figure 11
Comparison of (a) plunging response, (b) pitch response obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO based IMC-PID control of MIMO nonlinear wing vibration system.
(a)
(b)


## 4.4. Robustness Test

In order to further demonstrate the robustness of GNPSO, the comparison of system responses under varied system parameters obtained by GA, PSO, GPSO, NPSO, DE–NLJ GADE and GNPSO for nonlinear MIMO wing vibration system are presented in this subsection. The variations of system parameters are given in Table17. The robust performance of all algorithms is evaluated by the indices of I
T
A
E
h and I
T
A
E
θ. Table 17 presents the robust performance of nonlinear MIMO wing vibration control system under varied system parameters, obtained by different algorithms. Clearly, GNPSO performs best among seven algorithms since the values of I
T
A
E
h and I
T
A
E
θ vary in relatively small range according to variations of system parameters. For GA, PSO, GPSO and GADE, they cannot maintain good robustness with significant change of I
T
A
E
θ from case to case. DE–NLJ and NPSO have doubled the value of I
T
A
E
h when it comes to the parameter change of c
l
,
β
,
c
m
,
β. Hence, GNPSO obtains stronger robustness than other six algorithms due to the smaller variations in performance indices under varied system parameters.Table 17
Robust performance of GA, PSO, NPSO, GPSO, DE–NLJ, GADE, and GNPSO for MIMO nonlinear wing vibration control system.
Parameter variation
U
C
h
,
C
θ
c
l
,
θ
,
c
m
,
θ
c
l
,
β
,
c
m
,
β
% change
−10%
+10%
−10%
+10%
Algorithm
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
h
I
T
A
E
θ
I
T
A
E
h
I
T
A
E
θ
GA
3.3648
113.0802
3.4791
121.0771
3.5214
118.4546
8.0568
242.6463
PSO
2.2940
25.4754
2.4915
28.4621
3.5588
41.9787
6.3555
89.2480
NPSO
0.4015
20.1432
0.5364
21.4972
0.4227
19.6518
0.9005
24.0710
GPSO
0.9014
19.0715
2.3205
37.1065
2.0433
30.8812
2.9261
48.9147
DE–NLJ
0.3434
17.7616
0.5704
19.8848
0.4441
18.3618
1.0579
20.2744
GADE
3.7324
123.2064
4.1695
142.4313
4.4529
144.3820
15.5231
463.8153
GNPSO
0.2510
15.2886
0.2846
17.5032
0.2519
16.2986
0.3190
18.6249


## 4.5. Computational Efficiency

In order to further analyze the computational efficiency of GNPSO, the statistical results of the computation time obtained by GA, PSO, NPSO, GPSO, DE–NLJ, GADE and GNPSO on 25 independent runs are summarized in Table18, Table 19 and Table 20 for Case 1, Case 2, and Case 3 respectively. The statistical analysis of computation time for each case is implemented based on performance indices, including the minimum (t
m
i
n), the maximum (t
m
a
x), average (t
a
v
e
r
.), median (t
m
e
d
i
.) and standard deviation (t
sd). In practice, the fast tuning of the IMC-PID controller is favored to avoid the large time delay of the control signal and the uncontrolled system in possible damage. For the previous study of GPSO based PID control issue [24], the computation complexity is nearly O(PGn) with the population size of P, max. iteration number of G and the number of control parameters n, while those settings could cost large computation time of more than 200 seconds for complex wing vibration control systems. In this paper, to ensure the fast tuning of the IMC-PID controller, the computational cost is approximated as O(0.3
P
0.4
G
0.3
n) with reduced population size/iteration number and less control parameter. Tables 18–20 reveal that computation time obtained by the algorithms for complex wing vibration control systems are relatively small, due to reduced computation complexity. Clearly, among all algorithms, GNPSO performs best in terms of all performance indices for Cases 2 and 3, and three out of five indices for Case 1. Thus, GNPSO has better computational efficiency than other six algorithms on diverse wing vibration systems.Table 18
Statistical results of computation time (sec) using different algorithms for Case 1: linear wing vibration system with single control surface.
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
t
m
i
n
38.0083
36.4092
36.9464
41.9241
46.0005
35.4421
t
m
a
x
45.8949
38.5420
46.0858
49.6778
51.8526
46.4936
t
a
v
e
r
.
43.0086
37.4733
37.9925
44.1037
47.8175
36.6189
t
m
e
d
i
.
43.3580
37.4422
37.6460
43.6939
47.2339
36.3542
t
s
d
1.8242
0.4833
1.7357
1.7207
1.8155
2.1087
The bold values mean the best values among all the algorithms.Table 19
Statistical results of computation time (sec) using different algorithms for Case 2: nonlinear wing vibration system with single control surface.
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
t
m
i
n
41.5313
36.8535
36.8671
44.4655
53.3908
35.4934
t
m
a
x
44.2525
39.8687
40.1606
48.1763
60.1680
38.3784
t
a
v
e
r
.
42.4075
37.9693
38.4001
46.1848
55.0292
37.2044
t
m
e
d
i
.
42.3507
37.7893
38.4244
46.1072
54.5242
37.2493
t
s
d
0.7156
0.7506
0.7647
1.1402
1.5666
0.6607
The bold values mean the best values among all the algorithms.Table 20
Statistical results of computation time (sec) using different algorithms for Case 3: MIMO nonlinear wing vibration system with two control surfaces.
PSO
NPSO
GPSO
DE–NLJ
GADE
GNPSO
t
m
i
n
51.9130
46.5234
48.7196
55.4146
59.2800
45.8710
t
m
a
x
73.7114
61.4269
62.4517
64.0192
83.6449
51.5678
t
a
v
e
r
.
62.6266
52.2032
51.5982
58.5969
64.6545
47.6273
t
m
e
d
i
.
61.2508
49.2194
50.9832
58.2329
62.6389
47.3272
t
s
d
6.0806
5.0044
2.7682
1.9535
5.8313
1.0530
The bold values mean the best values among all the algorithms.


## 5. Conclusions

In this paper, a hybrid algorithm of modified Grey-PSO and NLJ, called GNPSO is proposed to search for optimal parameters of IMC-PID controller for different types of wing vibration systems. The key features of proposed GNPSO involve efficient global search by designing new grey analysis rule, improved balance between exploration and exploitation by employing new updating formula of acceleration coefficients, enhanced local search by involving NLJ, and reduced computation cost by using small population/iteration number. The superiority of proposed GNPSO to classic algorithms (GA, PSO), other hybrid algorithms (DE–NLJ, NPSO) and recently published grey-based algorithms (GPSO, GADE) is demonstrated by simulation tests on a wide range of wing vibration control systems, including linear system, nonlinear system and MIMO system.Simulation and statistical results reveal that GNPSO optimized method outperforms other six algorithms based methods on better convergence behavior, improved vibration control performance, and good applicability to a wide range of wing vibration systems. Especially, for nonlinear wing vibration system, GNPSO could meet most of control aims while existing algorithms, such as GA, PSO, NPSO and DE–NLJ tend to meet one control objective by sacrificing others. For MIMO nonlinear wing vibration system, proposed technique ensures good vibration control performance while the unsatisfied control performance is obtained by GA, PSO, GADE, GPSO and DE–NLJ. Besides, reduced computation time of parameters optimization is obtained by GNPSO on all system cases, compared to other existing optimization algorithms. In future work, it is worth to further develop the GNPSO based tuning method for under-actuated wing vibration control systems, regarding cost function, algorithm design and operation selections.



---

*Source: 8283178-2019-12-15.xml*