# Active Disturbance Rejection Attitude Control for Hypersonic Vehicle Based on Intelligent Stochastic Robust Optimization Method

Authors: Bo Wang; Wei Liu; Zhongtao Cheng; Lei Liu; Yongji Wang
Journal: Complexity
 (2020)

---

## Abstract

This paper designs a double-loop cascade active disturbance rejection control (ADRC) to overcome the external disturbances and parameter uncertainty during hypersonic vehicle flight. The vehicle attitude angle and attitude angular velocity are regulated in outer loop and inner loop, respectively. A stochastic robust approach is employed to further tune the ADRC parameters for better control performances. The Monte Carlo sampling of uncertain parameter is adopted to evaluate the stochastic robust performance. An improved differential evolution algorithm that combines neighborhood field optimization and triangular mutation is employed as the numerical solver. Simulation results show that the ADRC controller with optimized parameters manifests improved robustness as well as good control performances.


---

## Body

## 1. Introduction

The development of hypersonic vehicle (HV) technology has received increasing attention over the past decades. One of the key issues in this area is the attitude control. The difficulty of HV attitude control is incurred from the high nonlinearity, strong coupling, uncertainty, and external disturbance. Due to the high speed and the significant interactions across the aerodynamics, propulsion system, and structural dynamics, HVs are sensitive to uncertainties [1], which are inevitable in the complicated and uncertain flight environment. As a result, developing a robust and adaptive control for HVs under uncertainties and external disturbance becomes one of the major focuses.Wang and Stengel [2] summarized the robust flight control system of HVs. There are several kinds of control approaches proposed in HV attitude control, including the robust control, adaptive control, sliding mode control, and active disturbance rejection control (ADRC). In order to control the nonlinear system, feedback linearization was employed to transform the HV model into an equivalent linear model. Parker et al. [3] proposed a feedback linearization-based robust minimax LQR control method for HV under uncertainties. The sliding mode control (SMC) is also an effective method to control the system with uncertainties. Shtessel et al. [4] proposed a double-loop SMC with satisfied robust performance, in which the angular velocity and attitude angular are tracked in inner loop and outer loop, respectively. In presence of both uncertainties and external disturbances, stabilized velocity and altitude tracking was achieved by the adaptive terminal sliding mode control (SMC) of HV, where disturbances were estimated and compensated with a nonlinear disturbance observer (DOB) [5, 6].Han [7] proposed the scheme of ADRC which estimates and compensates the uncertainties and disturbances by an extended state observer (ESO). The ADRC becomes more and more popular due to its inherent robustness against nonlinearity and disturbances without relying on an accurate mathematical model of the system [8, 9]. There have been a number of achievements of employing ADRC in the HV control area. Yao and Wang applied ADRC to flight vehicle attitude control [10], Chen studied a robust ADRC control of re-entry vehicle [11], Zhang designed an optimization approach for HVs’ trajectory linearization ADRC [12], Huo proposed a drag profile tracking law based on ADRC control for re-entry vehicle [13], and Chao studied the velocity control based on active disturbance rejection for air-breathing supersonic vehicles [14].In ADRC, the ESO estimates the unmeasured uncertainties and disturbances of the controlled system, such that they can be compensated in the control law design, so as to improve the performance of the control system. Designing an ESO in ADRC for HV control involves tuning a lot of parameters. Although ADRC is inherent robust, the ADRC parameters can influence the ADRC stability and ESO convergence, which is the key issue in its successful applications. Wan [15] analyzed the stability of the active disturbance controller based on the absolute stability theory, and gave the necessary and sufficient conditions for the stability of the first-order ADRC system. The ADRC stability in HV control is even more complicated, due to the multiinput multioutput (MIMO) and nonholonomic nature of the system. Guo et al. [16, 17] have proposed the proof of MIMO ADRC stability. Chen [18] studied robust stabilization of extended nonholonomic chained-form systems with dynamic nonlinear uncertain terms by using ADRC.For flight control, especially HV control, it remains an open problem to choose proper parameters to guarantee the system stability and further dynamic performances. It is worth studying introducing the optimization technologies to improve the ADRC performances for flight control [19]. Specifically, the stochastic robust optimization is a promising perspective to find optimal ADRC parameters that satisfy more control performance requirements. Based on the Monte Carlo method, the chance of a set of parameters to manifest satisfying performance against uncertainties and disturbances can be estimated and employed as the criterion of choosing parameters, such that the system performance in most situations can be guaranteed. However, the involvement of nonlinearity, MIMO, fast time variation, and strong coupling in the controlled system incurs a very complicated optimization problem, where the effectiveness of conventional mathematical programming approaches is limited.The bioinspired optimization, a.k.a. the intelligent optimization technique, has been proved to be effective against complicated problems in HV control area [20, 21] and many other engineering fields [22–29]. Among a variety of methods, the differential evolution algorithm is one of the most widely adopted algorithms. It is easy to implement by computers, manifests excellent convergence speed in complicated optimization problems [30]. Some recent results of the DE algorithm further enhanced its superiority. Zhang and Wu [31, 32] proposed a neighborhood field cooperation mechanism to improve the convergence speed and global search ability. Mohamed [33] proposed a triangular mutation in the DE algorithm to accelerate the convergence while avoiding premature convergence. The combination of stochastic robust approach and intelligent optimization can facilitate the stochastic robust ADRC parameter optimization problem.This paper mainly investigates the design of a double-loop cascade active disturbance rejection control for hypersonic vehicles, where the angular velocity ADRC and attitude angular ADRC are designed, respectively. During the controller design process, the coefficients of inner and outer ESO and the parameters of related nonlinear controller are the tuned via a stochastic robust optimization approach, which guarantees the good performance of ADRC system against unmeasured uncertainties and disturbances. In order to solve the optimization problem, an improved differential evolution (DE) algorithm that incorporates neighborhood field optimization and triangular mutation is employed as the numerical solver. Simulation results illustrate the effectiveness of the proposed approach.The structure of this paper is as follows. The dynamic model of hypersonic vehicle is introduced in Section2, the design of cascade ADRC is discussed in Section 3, the control optimization based on stochastic robust is carried in Section 4. The simulation results are illustrated in Section 5. Finally, the conclusion of the work is proposed in Section 6.


## 2. Hypersonic Vehicle and Uncertainty Modelling

Before the implement of the controller design, the hypersonic vehicle should be established. A generic hypersonic vehicle (HV) model proposed by the American Langley Institute is hereby employed [34]. It should be noted that the inertia product between different directions is ignored, and the model established in the speed coordinate can be expressed as follows:(1)α˙=−pcosαtanβ+q+rsinαtanβ−LMVcosβ+GvzMVcosβ,β˙=psinαtanβ−rcosα+YMV+GvyMV,ϕ˙=p+rcosϕtanθ+qsinϕtanθ,q˙=Iz−IxIyrp+mIy,r˙=Ix−IyIzpq+nIz,p˙=Iy−IzIxqr+lIx,with(2)D=CDq¯S,Y=CYq¯S,L=CLq¯S,m=Cmq¯cS+xcgDsinα+Lcosα,n=Cnq¯bS−xcgY,l=Clq¯bS,where α and β denote the angle of attack and angle of side slip, respectively. θ and ϕ denote the pitch angle and roll angle, respectively. q, r, and p denote the pitch rate, yaw rate, and roll rate, respectively. D, Y, and L denote the drag force, yawing force, and lift force, respectively. m, n, and l denote the pitch moment, yawing moment, and rolling moment, respectively. q¯=1/2ρV2 is the dynamic pressure, with ρ, the density of air. S is the reference area of the HV, and M is the mass. V denotes the HV speed. Gvy and Gvz are the gravity vector components. b is the wing length, and c is the mean aerodynamic chord. CD,CY,CL,Cm,Cn,andCl are a set of coefficients for D,Y,L,m,n,andl. The uncertainties come with the mass M, air density ρ, and the coefficients CD,CY,CL,Cm,Cn,andCl. The following equation is employed to represent the calculation of the eight considered influence factors:(3)M=1+ΔMM0,ρ=1+Δρρ0,CD=1+ΔCDCD0,CY=1+ΔCYCY0,CL=1+ΔCLCL0,Cm=1+ΔCmCm0,Cn=1+ΔCnCn0,Cl=1+ΔClCl0,i.e., let ⋅ denote an arbitrary influenced parameter, then the parameter is computed by ⋅=1+Δ⋅⋅0, with the parametric uncertainty Δ⋅ and nominal value ⋅0.


## 3. Cascade ADRC Design

A general ADRC is composed of four parts, which are the scheduling transition process, the ESO, the nonlinear combination, and the disturbance compensation. The main idea of the ADRC is to treat all the parametric uncertainty estimations as unknown disturbances, which will be compensated upon the design of the control input [7].


### 3.1. Design of Cascade ADRC

When applied to the attitude angle control of HV, the classical second-order single-loop ADRC ignores a large number of uncertainties, and the disturbance observed by the ESO is too complex, resulting in undesirable large chattering in the control output. Hence, a cascade ADRC design is employed, where the inner loop and outer loop are both cascade first-order systems, such that the influence of uncertainty on the control output is reduced. The basic structure of the double-loop ADRC is shown in Figure1.Figure 1
Basic structure of the double-loop cascade ADRC.This structure is based on the idea of the backstepping method. Firstly, the angular velocity of HV in each channel is determined by the ADRC in the outer loop, which will be taken as the set point for the inner loop. Thereafter, according to the ADRC in the inner loop, the control quantity of rudder is determined. Finally, the control input is acquired from reverse recursion.To fulfill the design of the double-loop cascade controller, 1 is divided into inner and outer loop control systems. The inner loop system controls the angular velocity of HV, and the outer loop control system controls the attitude angle of HV. Control systems can be simplified into(4)α˙=f1α,β,δe,δa,θ,ϕ,M,V+b11q,β˙=f2α,β,δe,δa,δr,θ,ϕ,M,V+b21r,ϕ˙=f3α,β,δe,δa,δr,θ,ϕ,M,V,r,q+b31p,q˙=g1α,β,δe,δa,δr,θ,ϕ,M,V,r,p+b12δα,r˙=g2α,β,δe,δa,δr,θ,ϕ,M,V,p,q+b22δβ,p˙=g3α,β,δe,δa,δr,θ,ϕ,M,V,r,q+b32δϕ,where f1,f2,f3,g1,g2,andg3 are unknown uncertainty functions. b11,b21,b31,b12,b22,andb32 are the compensation factors of control. And δα,δβ,andδϕ are the virtual rudder angle of the pitch channel, yaw channel, and rolling channel, respectively.The pitch channel is adopted as an example to demonstrate the controller design process. Firstly, an outer loop ADRC is designed in a continuous manner. Since the system is a first-order system, the outer loop ADRC controller only needs position information for feedback. Therefore, the transition process only needs the transition signal of the given signal and does not need the differential signal of the given signal, such that the FAL function can be adopted [7]. The specific form of the designed transition process is as follows:(5)e01=x1−α0,fh=fale01,r01,h,x˙1=−fh,where the FAL function is defined as follows:(6)fale,r,h=e/hr−1,e≤h,ersigne,e>h.The extended state observer along with the nonlinear feedback and disturbance compensation are thereby designed:(7)e11=z11−α,fe1=fale11,0.5,h,z˙11=z12−β11e11+b11u11,z˙12=−β12fe1,es11=x1−z11,u11=k11fales11,r11,h−z12b11.In this way, the outer loop ADRC controller for pitch channel is obtained from the transition process of the outer loop, the extended state observer, nonlinear feedback, and disturbance compensation, and the dynamic analysis is carried out for the above parts:(8)α˙=f1+b11⋅k11fales11,r11,h−z12≈b11⋅k11fales11,r11,h.Equations (5)–(7) reveal that the original nonlinear system is transformed into a first-order pure integral system, which is controlled by nonlinear feedback antidisturbance controller. Therefore, the outer loop’s ADRC controller can better design the virtual control. The nonlinear system for controlling the inner loop is as follows:(9)e12=s11−q,fe2=fale12,0.5,h,s˙11=s12−β13e12+b12u1,s˙12=−β14⋅fe2,es12=u11−s11,u1=k12fales12,r12,h−s12b12.It should be noted that the transition process is eliminated in the ADRC of the inner loop to reduce unnecessary chattering. Further explanations of the ADRC design can be referred to [7, 9].Combining the abovementioned ADRC for the outer and inner loops, the cascade ADRC structure is determined, as shown in Figure2. The double-loop ADRC converts the original system into a pure integral link system. By controlling the pure integral system through nonlinear feedback method, the inner and outer loops of the system remain stable. In this way, the entire flight system can be kept stable.Figure 2
Cascade ADRC structure.


## 3.1. Design of Cascade ADRC

When applied to the attitude angle control of HV, the classical second-order single-loop ADRC ignores a large number of uncertainties, and the disturbance observed by the ESO is too complex, resulting in undesirable large chattering in the control output. Hence, a cascade ADRC design is employed, where the inner loop and outer loop are both cascade first-order systems, such that the influence of uncertainty on the control output is reduced. The basic structure of the double-loop ADRC is shown in Figure1.Figure 1
Basic structure of the double-loop cascade ADRC.This structure is based on the idea of the backstepping method. Firstly, the angular velocity of HV in each channel is determined by the ADRC in the outer loop, which will be taken as the set point for the inner loop. Thereafter, according to the ADRC in the inner loop, the control quantity of rudder is determined. Finally, the control input is acquired from reverse recursion.To fulfill the design of the double-loop cascade controller, 1 is divided into inner and outer loop control systems. The inner loop system controls the angular velocity of HV, and the outer loop control system controls the attitude angle of HV. Control systems can be simplified into(4)α˙=f1α,β,δe,δa,θ,ϕ,M,V+b11q,β˙=f2α,β,δe,δa,δr,θ,ϕ,M,V+b21r,ϕ˙=f3α,β,δe,δa,δr,θ,ϕ,M,V,r,q+b31p,q˙=g1α,β,δe,δa,δr,θ,ϕ,M,V,r,p+b12δα,r˙=g2α,β,δe,δa,δr,θ,ϕ,M,V,p,q+b22δβ,p˙=g3α,β,δe,δa,δr,θ,ϕ,M,V,r,q+b32δϕ,where f1,f2,f3,g1,g2,andg3 are unknown uncertainty functions. b11,b21,b31,b12,b22,andb32 are the compensation factors of control. And δα,δβ,andδϕ are the virtual rudder angle of the pitch channel, yaw channel, and rolling channel, respectively.The pitch channel is adopted as an example to demonstrate the controller design process. Firstly, an outer loop ADRC is designed in a continuous manner. Since the system is a first-order system, the outer loop ADRC controller only needs position information for feedback. Therefore, the transition process only needs the transition signal of the given signal and does not need the differential signal of the given signal, such that the FAL function can be adopted [7]. The specific form of the designed transition process is as follows:(5)e01=x1−α0,fh=fale01,r01,h,x˙1=−fh,where the FAL function is defined as follows:(6)fale,r,h=e/hr−1,e≤h,ersigne,e>h.The extended state observer along with the nonlinear feedback and disturbance compensation are thereby designed:(7)e11=z11−α,fe1=fale11,0.5,h,z˙11=z12−β11e11+b11u11,z˙12=−β12fe1,es11=x1−z11,u11=k11fales11,r11,h−z12b11.In this way, the outer loop ADRC controller for pitch channel is obtained from the transition process of the outer loop, the extended state observer, nonlinear feedback, and disturbance compensation, and the dynamic analysis is carried out for the above parts:(8)α˙=f1+b11⋅k11fales11,r11,h−z12≈b11⋅k11fales11,r11,h.Equations (5)–(7) reveal that the original nonlinear system is transformed into a first-order pure integral system, which is controlled by nonlinear feedback antidisturbance controller. Therefore, the outer loop’s ADRC controller can better design the virtual control. The nonlinear system for controlling the inner loop is as follows:(9)e12=s11−q,fe2=fale12,0.5,h,s˙11=s12−β13e12+b12u1,s˙12=−β14⋅fe2,es12=u11−s11,u1=k12fales12,r12,h−s12b12.It should be noted that the transition process is eliminated in the ADRC of the inner loop to reduce unnecessary chattering. Further explanations of the ADRC design can be referred to [7, 9].Combining the abovementioned ADRC for the outer and inner loops, the cascade ADRC structure is determined, as shown in Figure2. The double-loop ADRC converts the original system into a pure integral link system. By controlling the pure integral system through nonlinear feedback method, the inner and outer loops of the system remain stable. In this way, the entire flight system can be kept stable.Figure 2
Cascade ADRC structure.


## 4. Stochastic Robust Method Design

Maintaining system stability is a basic requirement for automatic control. The cascade ADRC can provide good antidisturbance performance. However, the flight control is such a complicated problem that a series of dynamic performances are also taken into account. As aforementioned, in order to further improve the control performance of ADRC, a stochastic robust method is employed to find the optimal ADRC parameters that satisfy as many dynamic performance requirements as possible.


### 4.1. Stochastic Robust Method

The stochastic robust method is composed of stochastic robust analysis (SRA) and stochastic robust design (SRD), as shown in Figure3. The Monte Carlo method is applied to randomly sample uncertain parameters, such that the system dissatisfaction probability under uncertainties can be identified. Such a probability is adopted to characterize the robustness of the system, and the confidence interval of the system dissatisfaction is adopted to describe the estimation accuracy.Figure 3
The stochastic robust optimization method.Based on uncertainty analysis of the HV in flight, the uncertainty parameter vector of the system is tabulated as Table1. The uncertainty parameter obeys the normal distribution of Nμ,σ2, where μ=0 and 3σ are uncertainty boundaries. Then, the uncertainty factor parameter vector is defined as(10)ν=m,ρ,CD,CY,CL,Cm,Cn,Cl,and parameters to be optimized:(11)d=r01,r02,r03,k11,k12,k21,k22,k31,k32,b11,b12,b21,b22,b31,b32.Table 1
Uncertain parameters and their distribution range.
Uncertainty parameterUncertainty range3σ (%)m, ρ10CD, CY, CL, Cm, Cn, Cl,30In order to make the optimization result more general, the step responses of the attack angle and the roll angle are employed. The parameters of the ESO are excluded at the current stage, as our major concern is to improve the transient performances of the controller. Including the ESO parameters can result in a much larger optimization problem, with less benefits to the investigated problem. The objective function is designed according to the step responses. For the investigated time-varying nonlinear flight system, the stability performance is judged from the integration of time and absolute error (ITAE) criterion:(12)SITAE=∫0ttetdt,for the sake of better transient oscillation performance. et denotes the error of the system output upon time t. Generally, a minimized SITAE is encouraged, where system stability is required.There are a number of performance requirements to be satisfied in HV control. LetSn denote the n-th performance obtained from a set of investigated parameters, and In is the performance index to characterize the system dissatisfaction. Conventionally, In should be binary, i.e., if Sn satisfies the requirement, In should be 0; otherwise, In=1. In order to improve the optimization effectiveness, a rising type judgment is here employed to decide In for the sake of reducing conservatism, as follows:(13)In=0,Sn≤an,12+12sinπbn−anSn−bn+an2,an<Sn≤bn,1,Sn>bn,n=1,2,3,….For convenience of derivation, assuming that all involved performance indicators are to be minimized.an and bn are boundaries of satisfaction and dissatisfaction. Sn≤an represents that performance Sn completely satisfies the requirement. Sn>bn represents the complete dissatisfaction. The sine type function is adopted such that a smooth, monotone increasing curve can be achieved between an and bn. The value of the function ranges between 0 and 1. As an example, Figure 4 depicts such a rising type judgment, with an=1 and bn=3.Figure 4
An example of the rising type judgment, wherean=1 and bn=3.A weighted sum approach is employed to obtain the following objective function:(14)J=∑i=1kwi∑n=1NIi,nN,where N denotes the number of samples, k is the number of performance requirements, and wi is the weight for the i-th performance requirements. For the investigated HV control problem, 16 performance indexes, corresponding to the step response of angle of attack and roll angle, are involved as follows:Step response of angle of attack:I1: Index of Angle of Attack ITAEI2: Index of Side-Slip Angle ITAEI3: Index of Roll Angle ITAEI4: Index of Pitch Angular Velocity Change ITAEI5: Index of Yaw Angular Velocity Change ITAEI6: Index of Roll Angular Velocity Change ITAEI7: Index of Adjustment time of angle of attackI8: Index of Overshoot of angle of attackStep response of roll angle:I9: Index of Angle of Attack ITAEI10: Index of Side-Slip Angle ITAEI11: Index of Roll Angle ITAEI12: Index of Pitch Angular Velocity Change ITAEI13: Index of Yaw Angular Velocity Change ITAEI14: Index of Roll Angular Velocity Change ITAEI15: Index of Adjustment time of roll angleI16: Index of Overshoot of roll angleThe weightsw−i imply the relative importance of the performance indexes. For the investigated HV control problem, the importance of the 16 performance indexes is considered to be the same, resulting in the same weights. For simplicity, all the weights are set to be 1. For each performance indicator, the ideal value is 0; therefore, an=0 with n=1,2,…,16. bn is set to be Sst,n, where Sst,n represents the standard value. Such standard values are obtained from a set of predecided parameters, which are usually manually tuned.If the parameter dimension to be optimized increases, the optimization of parameters will become larger. Hence, minimizing the parameters to be optimized is beneficial to obtain better control effects. According to the cascade ADRC structure, the controller parameters to be optimized are illustrated in Table2. The dimension for the parameters to be optimized is 15. Set the accuracy coefficient ε=0.01 and the confidence coefficient α=0.01. The following equation can be derived from the minimum sampling boundary theorem:(15)N≥ln1/αln1/1−ε=ln1/0.01ln1/1−0.01≈458.2,i.e., the uncertainty sample is at least 459 for the Monte Carlo sampling algorithm. If the number of samples reaches 459 and all sampling points meet the requirements, the probability that the HV control system meets the requirements is 99%. Hence, the Monte Carlo sampling size is set to be 460. Due to the computational complexity involved in the stochastic robust optimization, an improved DE algorithm is hereby employed as the numerical solver.Table 2
ADRC parameters to be optimized.
StructureParameters to be optimizedTransient Processr01r02r03Nonlineark11k21k31Combinationk12k22k32Disturbanceb11b21b31Compensationb12b22b32ChannelsPitchYawRolling


### 4.2. Improved DE Algorithm Incorporating NFO and Triangular Mutation

The complexity of the stochastic robust optimization results in demands in both efficient convergence and effective global search ability from the numerical solver. Furthermore, the objective function in (14) incurs high computational burden. Therefore, the neighborhood field optimization (NFO) and triangular mutation are incorporated into standard DE algorithm to improve the numerical solver efficiency.


#### 4.2.1. Neighborhood Filed Optimization

The NFO is inspired by the real world ecosystems, where a group of animals such as bees and birds are able to communicate and learn from their neighbors within limited perceptual range. Such a communication mechanism is introduced in the NFO method. During the search process, individuals learn from the local environment, instead of the globally “best” ones as in standard evolutionary algorithm. Specifically, an individual keeps following its superior neighbors and diverging from inferior neighbors [35]. In this way, an attractive field towards the superior neighbors and a repulsive field from the inferior neighbors collectively drive the individual to the global optima. According to Wu, the NFO is able to deliver promising results efficiently within acceptable computational time, thereby reduces the computational cost [31].


#### 4.2.2. Original Triangular Mutation

The original triangular mutation adopts the following mutation rule:(16)vg,j=x¯g,c+αL⋅xg,best−xg,better+αL⋅xg,better−xg,worst+αL⋅xg,best−xg,worst,with the triangular convex combination vector x¯g,c:(17)x¯g,c=w1⋅xg,best+w2⋅xg,better+w3⋅xg,worst,where αL represents the learning rate, which is a real number within 0,1. w1, w2, and w3 are nonnegative weights subject to ∑i=13wi=1. According to Mohamed [33], xg,best, xg,better, and xg,worst are the tournament best, better, and worst three randomly selected vectors. 16 can be rewritten into(18)vg,j=x¯g,c+2αL⋅xg,best−xg,worst,for the sake of amplifying the difference vector xg,best−xg,worst. When αL is larger, there are more mutant individuals, which is beneficial to population diversity, but it will reduce the convergence rate; when αL is smaller, the excellent attributes of the parent individuals are retained, but it is easy to fall into the local extreme point. The weights are computed as follows:(19)wi=pi∑i=13pi,i=1,2,3,where p1=1, p2=rand0.75,1, and p3=rand0.5,p2. randξl,ξu denotes that a real number is randomly generated between and without including a set of real numbers ξl and ξu. In this way, the weights of xg,best, xg,better, and xg,worst are in descending order to ensure greater influence from the superior individuals.


#### 4.2.3. Combined Mutation Strategy

The triangular mutation rule can be further improved via incorporating the NFO method. Letx^g,best, x^g,better, and x^g,worst denote the tournament best, better, and worst individuals that are randomly selected from the neighborhood field of the current individual x. The neighborhood field is defined as the nearest n superior individuals and m inferior individuals. The distance between individuals can vary upon different problems. In this issue, the Euclid distance is adopted.The mutation rule is thereby modified as follows:(20)vg,j=x¯gc+αL1⋅x^g,best−x^g,worst+αL2⋅x^g,best−x^g,better+αL3⋅x^g,better−x^g,worst,where αL1, αL2, and αL3 are learning rates that are computed as follows:(21)αLi=qi∑i=13qi,i=1,2,3,with q1=1, q2=rand0.75,1, and q3=rand0.5,q2, i.e., the learning rates are computed subject to the same strategy as the convex combination vector weights w1, w2, and w3. The convex combination vector x¯g,c is computed as follows:(22)x¯g,c=w1⋅x^g,best+w2⋅x^g,better+w3⋅x^g,worst,where the weights are computed as 19 indicates. To be noticed, pi and qi are independent and generated upon each iteration.


#### 4.2.4. Strategy Selection

In order to keep the balance between the global search and local search, a strategy selection mechanism is added, where the improved DE algorithm and standard DE algorithm are combined. The selection is made based on a nonlinear decreasing probability rule:(23)rand0,1≥1−GGEN2,where G denotes the current number of iterations, and GEN is the maximum iteration number. The probability that (23) is true keeps increasing upon the progress of iterations. Combined mutation strategy (20) applies if (23) is true. Otherwise, the mutation strategy from the original NFO-DE algorithm from [31] is adopted as follows:(24)vg,j=xg,j+αL⋅rand0,1⋅xcg,j−xg,j+αL⋅rand0,1⋅xg,j−xwg,j,where xcg,j denotes the nearest superior neighbor of xg,j, and xwg,j is the nearest inferior neighbor. The standard NFO-DE mutation strategy is more likely to take place at the beginning of the iteration, for the sake of global search ability. In the middle of the iteration, in order to improve the search ability and ensure that the algorithm can converge to optima, the probability of adopting the combined strategy is gradually increasing. In the middle and final stages, the possibility of choosing combined strategy is greater than that of standard NFO mutation strategy. In this way, the global search and fast convergence are considered in parallel.


#### 4.2.5. Pseudocode

The pseudocode of the DE algorithm that incorporates the combined mutation strategy is proposed in Algorithm1. For convenience, d¯x¯1,x¯2 is adopted to denote the Euclid distance between vectors x¯1 and x¯2.Algorithm 1: Pseudocode of the DE algorithm with combined mutation strategy.
Definition:np: the population size;GEN: the maximum number of generations for stopping criterion;d: dimension of the problem;x: the decision matrix with the size of GEN ∗ np ∗ d;J: the function value vector with the size of np ∗ 1;αL: the learning rate for conventional mutation strategy.(1)
BEGIN(2)
Set mutation probabilityCr and learning rate αL;(3)
Create a randomly initialized population {x1,jj=1,2,…,np};(4)
Letxbest=x1,1;(5)
whileg=1toGEN do(6)
whilej=1tonp do(7)
Locatexg,j in X and obtain its n nearest superior neighbors xcg,ji and inferior neighbors xwg,ji, with i=1,…,n;(8)
Selectx^g,best as the tournament best from (xcg,ji) and x^g,worst as the tournament worst from (xwg,ji), with i=1,…,n;(9)
Selectx^g,better as the argminid¯xg,j,xcg,ji;(10)
Generate three random valuesq1, q2, and q3, where q1=1, q2=rand0.75,1, and q3=rand0.5,q2;(11)
Compute the convex combination vector weightsw1, w2, and w3, according to wi=qi/∑i=13qi,i=1,2,3;(12)
Obtain the combination vector:x¯g,c=w1⋅x^g,best+w2⋅x^g,better+w3⋅x^g,worst;(13)
Generate the random valuesq1, q2, and q3, where q1=1, q2=rand0.75,1, and q3=rand0.5,q2 again;(14)
Compute the learning ratesαL1, αL2, and αL3, according to αLi=qi/∑i=13qi,i=1,2,3;(15)
ifrand0,1≥1−g/GEN2then(16)
The triangular mutation vectorvg,j=x¯g,c+αL1⋅x^g,best−x^g,worst+αL2⋅x^g,best−x^g,better+αL3⋅x^g,better−x^g,worst;(17)
else(18)
The conventional mutation vectorvg,j=xg,j+αL⋅rand0,1⋅xcg,j−xg,j+αL⋅rand0,1⋅xg,j−xwg,j;(19)
end if(20)
Repairvg,j if it violates the upperbound or lowerbound;(21)
Generateirand=randint1,d;(22)
whilei=1toddo(23)
ifi=irand or rand0,1<CR then(24)
vg,j,i=vg,j,i;(25)
else(26)
vg,j,i=xg,j,i;(27)
end if(28)
end while(29)
ifJvg,j≤Jxg,j then(30)
xg+1,j=vg,j;(31)
ifJxg+1,j<Jxbest then(32)
xbest=xg+1,j;(33)
end if(34)
else(35)
xg+1,j=xg,j;(36)
end if(37)
end while(38)
end while(39)
Returnxbest;(40)
ENDRemark 1.
The proposed algorithm combines a newly proposed mutation rule based on the convex combination vector of the triplet as defined by the three vectors and the difference vector between the best and the worst individuals among the three randomly selected vectors from the neighborhood field. The combined mutation strategy manifests advantages in accuracy and convergence speed on the investigated stochastic robust optimization problem. It has yet to verify its superiority on generic issues, e.g., the benchmark problems.


## 4.1. Stochastic Robust Method

The stochastic robust method is composed of stochastic robust analysis (SRA) and stochastic robust design (SRD), as shown in Figure3. The Monte Carlo method is applied to randomly sample uncertain parameters, such that the system dissatisfaction probability under uncertainties can be identified. Such a probability is adopted to characterize the robustness of the system, and the confidence interval of the system dissatisfaction is adopted to describe the estimation accuracy.Figure 3
The stochastic robust optimization method.Based on uncertainty analysis of the HV in flight, the uncertainty parameter vector of the system is tabulated as Table1. The uncertainty parameter obeys the normal distribution of Nμ,σ2, where μ=0 and 3σ are uncertainty boundaries. Then, the uncertainty factor parameter vector is defined as(10)ν=m,ρ,CD,CY,CL,Cm,Cn,Cl,and parameters to be optimized:(11)d=r01,r02,r03,k11,k12,k21,k22,k31,k32,b11,b12,b21,b22,b31,b32.Table 1
Uncertain parameters and their distribution range.
Uncertainty parameterUncertainty range3σ (%)m, ρ10CD, CY, CL, Cm, Cn, Cl,30In order to make the optimization result more general, the step responses of the attack angle and the roll angle are employed. The parameters of the ESO are excluded at the current stage, as our major concern is to improve the transient performances of the controller. Including the ESO parameters can result in a much larger optimization problem, with less benefits to the investigated problem. The objective function is designed according to the step responses. For the investigated time-varying nonlinear flight system, the stability performance is judged from the integration of time and absolute error (ITAE) criterion:(12)SITAE=∫0ttetdt,for the sake of better transient oscillation performance. et denotes the error of the system output upon time t. Generally, a minimized SITAE is encouraged, where system stability is required.There are a number of performance requirements to be satisfied in HV control. LetSn denote the n-th performance obtained from a set of investigated parameters, and In is the performance index to characterize the system dissatisfaction. Conventionally, In should be binary, i.e., if Sn satisfies the requirement, In should be 0; otherwise, In=1. In order to improve the optimization effectiveness, a rising type judgment is here employed to decide In for the sake of reducing conservatism, as follows:(13)In=0,Sn≤an,12+12sinπbn−anSn−bn+an2,an<Sn≤bn,1,Sn>bn,n=1,2,3,….For convenience of derivation, assuming that all involved performance indicators are to be minimized.an and bn are boundaries of satisfaction and dissatisfaction. Sn≤an represents that performance Sn completely satisfies the requirement. Sn>bn represents the complete dissatisfaction. The sine type function is adopted such that a smooth, monotone increasing curve can be achieved between an and bn. The value of the function ranges between 0 and 1. As an example, Figure 4 depicts such a rising type judgment, with an=1 and bn=3.Figure 4
An example of the rising type judgment, wherean=1 and bn=3.A weighted sum approach is employed to obtain the following objective function:(14)J=∑i=1kwi∑n=1NIi,nN,where N denotes the number of samples, k is the number of performance requirements, and wi is the weight for the i-th performance requirements. For the investigated HV control problem, 16 performance indexes, corresponding to the step response of angle of attack and roll angle, are involved as follows:Step response of angle of attack:I1: Index of Angle of Attack ITAEI2: Index of Side-Slip Angle ITAEI3: Index of Roll Angle ITAEI4: Index of Pitch Angular Velocity Change ITAEI5: Index of Yaw Angular Velocity Change ITAEI6: Index of Roll Angular Velocity Change ITAEI7: Index of Adjustment time of angle of attackI8: Index of Overshoot of angle of attackStep response of roll angle:I9: Index of Angle of Attack ITAEI10: Index of Side-Slip Angle ITAEI11: Index of Roll Angle ITAEI12: Index of Pitch Angular Velocity Change ITAEI13: Index of Yaw Angular Velocity Change ITAEI14: Index of Roll Angular Velocity Change ITAEI15: Index of Adjustment time of roll angleI16: Index of Overshoot of roll angleThe weightsw−i imply the relative importance of the performance indexes. For the investigated HV control problem, the importance of the 16 performance indexes is considered to be the same, resulting in the same weights. For simplicity, all the weights are set to be 1. For each performance indicator, the ideal value is 0; therefore, an=0 with n=1,2,…,16. bn is set to be Sst,n, where Sst,n represents the standard value. Such standard values are obtained from a set of predecided parameters, which are usually manually tuned.If the parameter dimension to be optimized increases, the optimization of parameters will become larger. Hence, minimizing the parameters to be optimized is beneficial to obtain better control effects. According to the cascade ADRC structure, the controller parameters to be optimized are illustrated in Table2. The dimension for the parameters to be optimized is 15. Set the accuracy coefficient ε=0.01 and the confidence coefficient α=0.01. The following equation can be derived from the minimum sampling boundary theorem:(15)N≥ln1/αln1/1−ε=ln1/0.01ln1/1−0.01≈458.2,i.e., the uncertainty sample is at least 459 for the Monte Carlo sampling algorithm. If the number of samples reaches 459 and all sampling points meet the requirements, the probability that the HV control system meets the requirements is 99%. Hence, the Monte Carlo sampling size is set to be 460. Due to the computational complexity involved in the stochastic robust optimization, an improved DE algorithm is hereby employed as the numerical solver.Table 2
ADRC parameters to be optimized.
StructureParameters to be optimizedTransient Processr01r02r03Nonlineark11k21k31Combinationk12k22k32Disturbanceb11b21b31Compensationb12b22b32ChannelsPitchYawRolling


## 4.2. Improved DE Algorithm Incorporating NFO and Triangular Mutation

The complexity of the stochastic robust optimization results in demands in both efficient convergence and effective global search ability from the numerical solver. Furthermore, the objective function in (14) incurs high computational burden. Therefore, the neighborhood field optimization (NFO) and triangular mutation are incorporated into standard DE algorithm to improve the numerical solver efficiency.


### 4.2.1. Neighborhood Filed Optimization

The NFO is inspired by the real world ecosystems, where a group of animals such as bees and birds are able to communicate and learn from their neighbors within limited perceptual range. Such a communication mechanism is introduced in the NFO method. During the search process, individuals learn from the local environment, instead of the globally “best” ones as in standard evolutionary algorithm. Specifically, an individual keeps following its superior neighbors and diverging from inferior neighbors [35]. In this way, an attractive field towards the superior neighbors and a repulsive field from the inferior neighbors collectively drive the individual to the global optima. According to Wu, the NFO is able to deliver promising results efficiently within acceptable computational time, thereby reduces the computational cost [31].


### 4.2.2. Original Triangular Mutation

The original triangular mutation adopts the following mutation rule:(16)vg,j=x¯g,c+αL⋅xg,best−xg,better+αL⋅xg,better−xg,worst+αL⋅xg,best−xg,worst,with the triangular convex combination vector x¯g,c:(17)x¯g,c=w1⋅xg,best+w2⋅xg,better+w3⋅xg,worst,where αL represents the learning rate, which is a real number within 0,1. w1, w2, and w3 are nonnegative weights subject to ∑i=13wi=1. According to Mohamed [33], xg,best, xg,better, and xg,worst are the tournament best, better, and worst three randomly selected vectors. 16 can be rewritten into(18)vg,j=x¯g,c+2αL⋅xg,best−xg,worst,for the sake of amplifying the difference vector xg,best−xg,worst. When αL is larger, there are more mutant individuals, which is beneficial to population diversity, but it will reduce the convergence rate; when αL is smaller, the excellent attributes of the parent individuals are retained, but it is easy to fall into the local extreme point. The weights are computed as follows:(19)wi=pi∑i=13pi,i=1,2,3,where p1=1, p2=rand0.75,1, and p3=rand0.5,p2. randξl,ξu denotes that a real number is randomly generated between and without including a set of real numbers ξl and ξu. In this way, the weights of xg,best, xg,better, and xg,worst are in descending order to ensure greater influence from the superior individuals.


### 4.2.3. Combined Mutation Strategy

The triangular mutation rule can be further improved via incorporating the NFO method. Letx^g,best, x^g,better, and x^g,worst denote the tournament best, better, and worst individuals that are randomly selected from the neighborhood field of the current individual x. The neighborhood field is defined as the nearest n superior individuals and m inferior individuals. The distance between individuals can vary upon different problems. In this issue, the Euclid distance is adopted.The mutation rule is thereby modified as follows:(20)vg,j=x¯gc+αL1⋅x^g,best−x^g,worst+αL2⋅x^g,best−x^g,better+αL3⋅x^g,better−x^g,worst,where αL1, αL2, and αL3 are learning rates that are computed as follows:(21)αLi=qi∑i=13qi,i=1,2,3,with q1=1, q2=rand0.75,1, and q3=rand0.5,q2, i.e., the learning rates are computed subject to the same strategy as the convex combination vector weights w1, w2, and w3. The convex combination vector x¯g,c is computed as follows:(22)x¯g,c=w1⋅x^g,best+w2⋅x^g,better+w3⋅x^g,worst,where the weights are computed as 19 indicates. To be noticed, pi and qi are independent and generated upon each iteration.


### 4.2.4. Strategy Selection

In order to keep the balance between the global search and local search, a strategy selection mechanism is added, where the improved DE algorithm and standard DE algorithm are combined. The selection is made based on a nonlinear decreasing probability rule:(23)rand0,1≥1−GGEN2,where G denotes the current number of iterations, and GEN is the maximum iteration number. The probability that (23) is true keeps increasing upon the progress of iterations. Combined mutation strategy (20) applies if (23) is true. Otherwise, the mutation strategy from the original NFO-DE algorithm from [31] is adopted as follows:(24)vg,j=xg,j+αL⋅rand0,1⋅xcg,j−xg,j+αL⋅rand0,1⋅xg,j−xwg,j,where xcg,j denotes the nearest superior neighbor of xg,j, and xwg,j is the nearest inferior neighbor. The standard NFO-DE mutation strategy is more likely to take place at the beginning of the iteration, for the sake of global search ability. In the middle of the iteration, in order to improve the search ability and ensure that the algorithm can converge to optima, the probability of adopting the combined strategy is gradually increasing. In the middle and final stages, the possibility of choosing combined strategy is greater than that of standard NFO mutation strategy. In this way, the global search and fast convergence are considered in parallel.


### 4.2.5. Pseudocode

The pseudocode of the DE algorithm that incorporates the combined mutation strategy is proposed in Algorithm1. For convenience, d¯x¯1,x¯2 is adopted to denote the Euclid distance between vectors x¯1 and x¯2.Algorithm 1: Pseudocode of the DE algorithm with combined mutation strategy.
Definition:np: the population size;GEN: the maximum number of generations for stopping criterion;d: dimension of the problem;x: the decision matrix with the size of GEN ∗ np ∗ d;J: the function value vector with the size of np ∗ 1;αL: the learning rate for conventional mutation strategy.(1)
BEGIN(2)
Set mutation probabilityCr and learning rate αL;(3)
Create a randomly initialized population {x1,jj=1,2,…,np};(4)
Letxbest=x1,1;(5)
whileg=1toGEN do(6)
whilej=1tonp do(7)
Locatexg,j in X and obtain its n nearest superior neighbors xcg,ji and inferior neighbors xwg,ji, with i=1,…,n;(8)
Selectx^g,best as the tournament best from (xcg,ji) and x^g,worst as the tournament worst from (xwg,ji), with i=1,…,n;(9)
Selectx^g,better as the argminid¯xg,j,xcg,ji;(10)
Generate three random valuesq1, q2, and q3, where q1=1, q2=rand0.75,1, and q3=rand0.5,q2;(11)
Compute the convex combination vector weightsw1, w2, and w3, according to wi=qi/∑i=13qi,i=1,2,3;(12)
Obtain the combination vector:x¯g,c=w1⋅x^g,best+w2⋅x^g,better+w3⋅x^g,worst;(13)
Generate the random valuesq1, q2, and q3, where q1=1, q2=rand0.75,1, and q3=rand0.5,q2 again;(14)
Compute the learning ratesαL1, αL2, and αL3, according to αLi=qi/∑i=13qi,i=1,2,3;(15)
ifrand0,1≥1−g/GEN2then(16)
The triangular mutation vectorvg,j=x¯g,c+αL1⋅x^g,best−x^g,worst+αL2⋅x^g,best−x^g,better+αL3⋅x^g,better−x^g,worst;(17)
else(18)
The conventional mutation vectorvg,j=xg,j+αL⋅rand0,1⋅xcg,j−xg,j+αL⋅rand0,1⋅xg,j−xwg,j;(19)
end if(20)
Repairvg,j if it violates the upperbound or lowerbound;(21)
Generateirand=randint1,d;(22)
whilei=1toddo(23)
ifi=irand or rand0,1<CR then(24)
vg,j,i=vg,j,i;(25)
else(26)
vg,j,i=xg,j,i;(27)
end if(28)
end while(29)
ifJvg,j≤Jxg,j then(30)
xg+1,j=vg,j;(31)
ifJxg+1,j<Jxbest then(32)
xbest=xg+1,j;(33)
end if(34)
else(35)
xg+1,j=xg,j;(36)
end if(37)
end while(38)
end while(39)
Returnxbest;(40)
ENDRemark 1.
The proposed algorithm combines a newly proposed mutation rule based on the convex combination vector of the triplet as defined by the three vectors and the difference vector between the best and the worst individuals among the three randomly selected vectors from the neighborhood field. The combined mutation strategy manifests advantages in accuracy and convergence speed on the investigated stochastic robust optimization problem. It has yet to verify its superiority on generic issues, e.g., the benchmark problems.


## 4.2.1. Neighborhood Filed Optimization

The NFO is inspired by the real world ecosystems, where a group of animals such as bees and birds are able to communicate and learn from their neighbors within limited perceptual range. Such a communication mechanism is introduced in the NFO method. During the search process, individuals learn from the local environment, instead of the globally “best” ones as in standard evolutionary algorithm. Specifically, an individual keeps following its superior neighbors and diverging from inferior neighbors [35]. In this way, an attractive field towards the superior neighbors and a repulsive field from the inferior neighbors collectively drive the individual to the global optima. According to Wu, the NFO is able to deliver promising results efficiently within acceptable computational time, thereby reduces the computational cost [31].


## 4.2.2. Original Triangular Mutation

The original triangular mutation adopts the following mutation rule:(16)vg,j=x¯g,c+αL⋅xg,best−xg,better+αL⋅xg,better−xg,worst+αL⋅xg,best−xg,worst,with the triangular convex combination vector x¯g,c:(17)x¯g,c=w1⋅xg,best+w2⋅xg,better+w3⋅xg,worst,where αL represents the learning rate, which is a real number within 0,1. w1, w2, and w3 are nonnegative weights subject to ∑i=13wi=1. According to Mohamed [33], xg,best, xg,better, and xg,worst are the tournament best, better, and worst three randomly selected vectors. 16 can be rewritten into(18)vg,j=x¯g,c+2αL⋅xg,best−xg,worst,for the sake of amplifying the difference vector xg,best−xg,worst. When αL is larger, there are more mutant individuals, which is beneficial to population diversity, but it will reduce the convergence rate; when αL is smaller, the excellent attributes of the parent individuals are retained, but it is easy to fall into the local extreme point. The weights are computed as follows:(19)wi=pi∑i=13pi,i=1,2,3,where p1=1, p2=rand0.75,1, and p3=rand0.5,p2. randξl,ξu denotes that a real number is randomly generated between and without including a set of real numbers ξl and ξu. In this way, the weights of xg,best, xg,better, and xg,worst are in descending order to ensure greater influence from the superior individuals.


## 4.2.3. Combined Mutation Strategy

The triangular mutation rule can be further improved via incorporating the NFO method. Letx^g,best, x^g,better, and x^g,worst denote the tournament best, better, and worst individuals that are randomly selected from the neighborhood field of the current individual x. The neighborhood field is defined as the nearest n superior individuals and m inferior individuals. The distance between individuals can vary upon different problems. In this issue, the Euclid distance is adopted.The mutation rule is thereby modified as follows:(20)vg,j=x¯gc+αL1⋅x^g,best−x^g,worst+αL2⋅x^g,best−x^g,better+αL3⋅x^g,better−x^g,worst,where αL1, αL2, and αL3 are learning rates that are computed as follows:(21)αLi=qi∑i=13qi,i=1,2,3,with q1=1, q2=rand0.75,1, and q3=rand0.5,q2, i.e., the learning rates are computed subject to the same strategy as the convex combination vector weights w1, w2, and w3. The convex combination vector x¯g,c is computed as follows:(22)x¯g,c=w1⋅x^g,best+w2⋅x^g,better+w3⋅x^g,worst,where the weights are computed as 19 indicates. To be noticed, pi and qi are independent and generated upon each iteration.


## 4.2.4. Strategy Selection

In order to keep the balance between the global search and local search, a strategy selection mechanism is added, where the improved DE algorithm and standard DE algorithm are combined. The selection is made based on a nonlinear decreasing probability rule:(23)rand0,1≥1−GGEN2,where G denotes the current number of iterations, and GEN is the maximum iteration number. The probability that (23) is true keeps increasing upon the progress of iterations. Combined mutation strategy (20) applies if (23) is true. Otherwise, the mutation strategy from the original NFO-DE algorithm from [31] is adopted as follows:(24)vg,j=xg,j+αL⋅rand0,1⋅xcg,j−xg,j+αL⋅rand0,1⋅xg,j−xwg,j,where xcg,j denotes the nearest superior neighbor of xg,j, and xwg,j is the nearest inferior neighbor. The standard NFO-DE mutation strategy is more likely to take place at the beginning of the iteration, for the sake of global search ability. In the middle of the iteration, in order to improve the search ability and ensure that the algorithm can converge to optima, the probability of adopting the combined strategy is gradually increasing. In the middle and final stages, the possibility of choosing combined strategy is greater than that of standard NFO mutation strategy. In this way, the global search and fast convergence are considered in parallel.


## 4.2.5. Pseudocode

The pseudocode of the DE algorithm that incorporates the combined mutation strategy is proposed in Algorithm1. For convenience, d¯x¯1,x¯2 is adopted to denote the Euclid distance between vectors x¯1 and x¯2.Algorithm 1: Pseudocode of the DE algorithm with combined mutation strategy.
Definition:np: the population size;GEN: the maximum number of generations for stopping criterion;d: dimension of the problem;x: the decision matrix with the size of GEN ∗ np ∗ d;J: the function value vector with the size of np ∗ 1;αL: the learning rate for conventional mutation strategy.(1)
BEGIN(2)
Set mutation probabilityCr and learning rate αL;(3)
Create a randomly initialized population {x1,jj=1,2,…,np};(4)
Letxbest=x1,1;(5)
whileg=1toGEN do(6)
whilej=1tonp do(7)
Locatexg,j in X and obtain its n nearest superior neighbors xcg,ji and inferior neighbors xwg,ji, with i=1,…,n;(8)
Selectx^g,best as the tournament best from (xcg,ji) and x^g,worst as the tournament worst from (xwg,ji), with i=1,…,n;(9)
Selectx^g,better as the argminid¯xg,j,xcg,ji;(10)
Generate three random valuesq1, q2, and q3, where q1=1, q2=rand0.75,1, and q3=rand0.5,q2;(11)
Compute the convex combination vector weightsw1, w2, and w3, according to wi=qi/∑i=13qi,i=1,2,3;(12)
Obtain the combination vector:x¯g,c=w1⋅x^g,best+w2⋅x^g,better+w3⋅x^g,worst;(13)
Generate the random valuesq1, q2, and q3, where q1=1, q2=rand0.75,1, and q3=rand0.5,q2 again;(14)
Compute the learning ratesαL1, αL2, and αL3, according to αLi=qi/∑i=13qi,i=1,2,3;(15)
ifrand0,1≥1−g/GEN2then(16)
The triangular mutation vectorvg,j=x¯g,c+αL1⋅x^g,best−x^g,worst+αL2⋅x^g,best−x^g,better+αL3⋅x^g,better−x^g,worst;(17)
else(18)
The conventional mutation vectorvg,j=xg,j+αL⋅rand0,1⋅xcg,j−xg,j+αL⋅rand0,1⋅xg,j−xwg,j;(19)
end if(20)
Repairvg,j if it violates the upperbound or lowerbound;(21)
Generateirand=randint1,d;(22)
whilei=1toddo(23)
ifi=irand or rand0,1<CR then(24)
vg,j,i=vg,j,i;(25)
else(26)
vg,j,i=xg,j,i;(27)
end if(28)
end while(29)
ifJvg,j≤Jxg,j then(30)
xg+1,j=vg,j;(31)
ifJxg+1,j<Jxbest then(32)
xbest=xg+1,j;(33)
end if(34)
else(35)
xg+1,j=xg,j;(36)
end if(37)
end while(38)
end while(39)
Returnxbest;(40)
ENDRemark 1.
The proposed algorithm combines a newly proposed mutation rule based on the convex combination vector of the triplet as defined by the three vectors and the difference vector between the best and the worst individuals among the three randomly selected vectors from the neighborhood field. The combined mutation strategy manifests advantages in accuracy and convergence speed on the investigated stochastic robust optimization problem. It has yet to verify its superiority on generic issues, e.g., the benchmark problems.


## 5. Simulation Verification

### 5.1. Cascade ADRC Simulation

Firstly, the transition process parameters are set as follows:r01=0.1, r02=0.1, and r03=0.01. ESO parameters are set as follows: βi1=100, βi2=300, βi3=100, and βi4=300, with i=1,2,3. Nonlinear feedback parameters are set as follows: k11=0.2, k21=0.01, k31=0.15, k12=20, k22=0.1, k32=10. ri1=0.5, and ri2=0.01, with i=1,2,3. The compensation factors are set as follows: b11=1, b12=0.007, b21=1, b22=−0.015, b31=1, and b32=0.0248. These parameters are manually designed and taken as the baseline parameters.


### 5.2. Stochastic Robust Optimization

The uncertainty parameter vector is shown in (10), and parameters to be optimized are shown in Table 2. The optimization upper bounds of the parameters to be optimized are the double of the values in the cascaded ADRC simulation, and the lower bounds are 50% of the values. Based on the proposed improved DE algorithm, the final objective function is obtained through performance indicators and weights. Taking the final objective function as the fitness value, the optimization parameters are optimized. The adopted DE algorithm is configured as follows. The population size and iteration numbers are decided according to the dimension of the decision variable vector. For the learning rates, multiple runs have been implemented to figure out the ones with better convergence. The problem dimension is 15; therefore, the population size is 75, and maximum iteration is 250. The learning rate αL is 0.5.Givend in (11), the final optimization result d∗ is as follows:(25)d∗=0.4581,0.0872,0.0056,0.2936,23.9940,0.0180,0.1026,0.4674,14.1997,1.0036,0.0128,1.0451,−0.0113,0.9556,0.0402,namely, the optimal parameters. Other parameters apart from d∗ adopt the same values in the previous subsection. The selected optimal parameters are the second-best result from 20 runs.In order to verify the robustness of the cascaded ADRC controller optimized by the stochastic robust method, the aerodynamic parameters are deviated in real time during the simulation. The adopted biased forms for the defection are as follows:(26)ΔM=0.1sint,Δρ=0.1sint,ΔCD=0.3sint+π3,ΔCY=0.3sint+2π3,ΔCL=0.3sint+π,ΔCm=0.3sint+4π3,ΔCn=0.3sint+5π3,ΔCl=0.3sint+2π.


### 5.3. Results and Analysis

According to above, there are four simulation cases. Firstly, the cascade ADRC simulation with baseline parameters is performed, namely, case (a), where uncertainties are ignored. This is to verify the effectiveness of the cascade ADRC design. For comparison, the uncertainties as characterized by 25 are introduced to perform the ADRC simulation with the baseline parameters, namely, case (b). In this way, the influence of the uncertainties is illustrated. Then, the stochastic robust optimization result, i.e., the optimal parameterd∗, is adopted, and the ADRC simulation is performed without uncertainties, namely, case (c). Finally, to investigate the robustness of d∗, the simulation with d∗ under uncertainties is performed, namely, case (d).There are three categories of control curves, as the output of the simulations. The first category is the attitude angle-tracking curves, including the attack angle-tracking curve, side-slip angle-tracking curve, and rolling angle-tracking curve. The results from the four cases are depicted in Figure5. The solid red lines represent the angle instruction signals, and the blue dashed lines represent the tracking curves.Figure 5
Attitude angle-tracking curves in the four cases: (a) attitude angle-tracking curve with baseline without uncertainties, (b) attitude angle-tracking curve with baseline under uncertainties, (c) attitude angle-tracking curve with optimal parameters without uncertainties, and (d) attitude angle-tracking curve with optimal parameters under uncertainties.
(a)(b)(c)(d)From Figure5, the following performances can be observed. In cases (a) and (b), the tracking tasks can be fulfilled, although overshoots and oscillations can take place. Furthermore, the parametric uncertainties deliver further oscillations during the tracking process. In cases (c) and (d), with optimal parameters, the control performances are further improved. The convergence speed is improved without overshoots, and oscillations are reduced. Even in case (d) where uncertainties take place, the control performances remain smooth.The angular velocity curves in the four cases are also depicted, as shown in Figure6. In case (a) and case (b), strong oscillations can be observed. Similar to the attitude angle-tracking curves in Figure 5, the convergence is guaranteed, whereas the transient processes are less smooth. In cases (c) and (d), the oscillations are reduced, curves become smoother, and introducing uncertainties only results in trivial changes. Furthermore, the steering curves of the servo in the four cases are depicted in Figure 7. The rudder deflections are significant with baseline parameters. When applying optimal parameters, the rudder actions are smoothen, where parametric uncertainties incur slight swing of the left and right rudders. The performances of the rudders are in accordance with the angular velocities against uncertainties.Figure 6
Angular velocity curves in the four cases: (a) angular velocity curve with baseline without uncertainties, (b) angular velocity curve with baseline under uncertainties, (c) angular velocity curve with optimal parameters without uncertainties, and (d) angular velocity curve with optimal parameters under uncertainties.
(a)(b)(c)(d)Figure 7
Steering curves of the servo in the four cases: (a) rudder deflection with baseline without uncertainties, (b) rudder deflection with baseline under uncertainties, (c) rudder deflection with optimal parameters without uncertainties, and (d) rudder deflection with optimal parameters under uncertainties.
(a)(b)(c)(d)According to the above figures, some advantages of the proposed cascade ADRC strategy and the stochastic robust optimization based controller parameter tuning in HVs control can be revealed. Firstly, the stability of the cascade ADRC strategy for HVs in complicated aerodynamic environments is verified. Secondly, the transient performances, even under the influence of aerodynamic parametric uncertainties, can be improved via the parameter optimization. While guaranteeing the convergence, overshoots are removed, oscillations are reduced, and control curves are smoothen with the optimal parameters. In this way, the design objective is achieved.


## 5.1. Cascade ADRC Simulation

Firstly, the transition process parameters are set as follows:r01=0.1, r02=0.1, and r03=0.01. ESO parameters are set as follows: βi1=100, βi2=300, βi3=100, and βi4=300, with i=1,2,3. Nonlinear feedback parameters are set as follows: k11=0.2, k21=0.01, k31=0.15, k12=20, k22=0.1, k32=10. ri1=0.5, and ri2=0.01, with i=1,2,3. The compensation factors are set as follows: b11=1, b12=0.007, b21=1, b22=−0.015, b31=1, and b32=0.0248. These parameters are manually designed and taken as the baseline parameters.


## 5.2. Stochastic Robust Optimization

The uncertainty parameter vector is shown in (10), and parameters to be optimized are shown in Table 2. The optimization upper bounds of the parameters to be optimized are the double of the values in the cascaded ADRC simulation, and the lower bounds are 50% of the values. Based on the proposed improved DE algorithm, the final objective function is obtained through performance indicators and weights. Taking the final objective function as the fitness value, the optimization parameters are optimized. The adopted DE algorithm is configured as follows. The population size and iteration numbers are decided according to the dimension of the decision variable vector. For the learning rates, multiple runs have been implemented to figure out the ones with better convergence. The problem dimension is 15; therefore, the population size is 75, and maximum iteration is 250. The learning rate αL is 0.5.Givend in (11), the final optimization result d∗ is as follows:(25)d∗=0.4581,0.0872,0.0056,0.2936,23.9940,0.0180,0.1026,0.4674,14.1997,1.0036,0.0128,1.0451,−0.0113,0.9556,0.0402,namely, the optimal parameters. Other parameters apart from d∗ adopt the same values in the previous subsection. The selected optimal parameters are the second-best result from 20 runs.In order to verify the robustness of the cascaded ADRC controller optimized by the stochastic robust method, the aerodynamic parameters are deviated in real time during the simulation. The adopted biased forms for the defection are as follows:(26)ΔM=0.1sint,Δρ=0.1sint,ΔCD=0.3sint+π3,ΔCY=0.3sint+2π3,ΔCL=0.3sint+π,ΔCm=0.3sint+4π3,ΔCn=0.3sint+5π3,ΔCl=0.3sint+2π.


## 5.3. Results and Analysis

According to above, there are four simulation cases. Firstly, the cascade ADRC simulation with baseline parameters is performed, namely, case (a), where uncertainties are ignored. This is to verify the effectiveness of the cascade ADRC design. For comparison, the uncertainties as characterized by 25 are introduced to perform the ADRC simulation with the baseline parameters, namely, case (b). In this way, the influence of the uncertainties is illustrated. Then, the stochastic robust optimization result, i.e., the optimal parameterd∗, is adopted, and the ADRC simulation is performed without uncertainties, namely, case (c). Finally, to investigate the robustness of d∗, the simulation with d∗ under uncertainties is performed, namely, case (d).There are three categories of control curves, as the output of the simulations. The first category is the attitude angle-tracking curves, including the attack angle-tracking curve, side-slip angle-tracking curve, and rolling angle-tracking curve. The results from the four cases are depicted in Figure5. The solid red lines represent the angle instruction signals, and the blue dashed lines represent the tracking curves.Figure 5
Attitude angle-tracking curves in the four cases: (a) attitude angle-tracking curve with baseline without uncertainties, (b) attitude angle-tracking curve with baseline under uncertainties, (c) attitude angle-tracking curve with optimal parameters without uncertainties, and (d) attitude angle-tracking curve with optimal parameters under uncertainties.
(a)(b)(c)(d)From Figure5, the following performances can be observed. In cases (a) and (b), the tracking tasks can be fulfilled, although overshoots and oscillations can take place. Furthermore, the parametric uncertainties deliver further oscillations during the tracking process. In cases (c) and (d), with optimal parameters, the control performances are further improved. The convergence speed is improved without overshoots, and oscillations are reduced. Even in case (d) where uncertainties take place, the control performances remain smooth.The angular velocity curves in the four cases are also depicted, as shown in Figure6. In case (a) and case (b), strong oscillations can be observed. Similar to the attitude angle-tracking curves in Figure 5, the convergence is guaranteed, whereas the transient processes are less smooth. In cases (c) and (d), the oscillations are reduced, curves become smoother, and introducing uncertainties only results in trivial changes. Furthermore, the steering curves of the servo in the four cases are depicted in Figure 7. The rudder deflections are significant with baseline parameters. When applying optimal parameters, the rudder actions are smoothen, where parametric uncertainties incur slight swing of the left and right rudders. The performances of the rudders are in accordance with the angular velocities against uncertainties.Figure 6
Angular velocity curves in the four cases: (a) angular velocity curve with baseline without uncertainties, (b) angular velocity curve with baseline under uncertainties, (c) angular velocity curve with optimal parameters without uncertainties, and (d) angular velocity curve with optimal parameters under uncertainties.
(a)(b)(c)(d)Figure 7
Steering curves of the servo in the four cases: (a) rudder deflection with baseline without uncertainties, (b) rudder deflection with baseline under uncertainties, (c) rudder deflection with optimal parameters without uncertainties, and (d) rudder deflection with optimal parameters under uncertainties.
(a)(b)(c)(d)According to the above figures, some advantages of the proposed cascade ADRC strategy and the stochastic robust optimization based controller parameter tuning in HVs control can be revealed. Firstly, the stability of the cascade ADRC strategy for HVs in complicated aerodynamic environments is verified. Secondly, the transient performances, even under the influence of aerodynamic parametric uncertainties, can be improved via the parameter optimization. While guaranteeing the convergence, overshoots are removed, oscillations are reduced, and control curves are smoothen with the optimal parameters. In this way, the design objective is achieved.


## 6. Conclusion

In this paper, a cascade ADRC is designed for high-speed, strong-coupling, fast time-varying, and strong-nonlinear high-speed HV. In order to obtain better dynamic performances and further enhance ADRC robustness, a stochastic robust optimization method is employed to further tune the ADRC parameters. The robust optimization objective function is designed to characterize the ability of the ADRC parameters to satisfy as many performance requirements as possible. An adaptive DE algorithm based on triangular mutation is employed to solve the complicated robust optimization problem. In simulations, the control performances from the baseline ADRC parameters and the robust optimization results are thoroughly compared. The control curves reveal that significant robustness improvement can be achieved via further parameter tuning in ADRC, such that the effectiveness and necessity of incorporating robust optimization into the ADRC controller design is identified and verified.



---

*Source: 3240405-2020-10-05.xml*