# Improved Hybrid Fireworks Algorithm-Based Parameter Optimization in High-Order Sliding Mode Control of Hypersonic Vehicles

Authors: Xiaomeng Yin; Xing Wei; Lei Liu; Yongji Wang
Journal: Complexity
 (2018)

---

## Abstract

With respect to the nonlinear hypersonic vehicle (HV) dynamics, achieving a satisfactory tracking control performance under uncertainties is always a challenge. The high-order sliding mode control (HOSMC) method with strong robustness has been applied to HVs. However, there are few methods for determining suitable HOSMC parameters for an efficacious control of HV, given that the uncertainties are randomly distributed. In this study, we introduce a hybrid fireworks algorithm- (FWA-) based parameter optimization into HV control design to satisfy the design requirements with high probability. First, the complex relation between design parameters and the cost function that evaluates the likelihood of system instability and violation of design requirements is modeled via stochastic robustness analysis. Subsequently, we propose an efficient hybrid FWA to solve the complex optimization problem concerning the uncertainties. The efficiency of the proposed hybrid FWA-based optimization method is demonstrated in the search of the optimal HV controller, in which the proposed method exhibits a better performance when compared with other algorithms.


---

## Body

## 1. Introduction

Hypersonic vehicles (HVs) have attracted increasing interest given their characteristics of high speed and excellent cost effectiveness to access the space. HVs usually fight in near space at a high speed, in which the aerodynamic properties are difficult to predict [1]. Additionally, owing to the peculiar structure of HVs, the couplings related to aerodynamics, propulsion, and structural dynamics are strong, and this makes HV sensitive to uncertainties [2]. In this study, we focus on the efficacious control design of nonlinear HV dynamics given that uncertainties are randomly distributed.As members of sliding mode control methods [3–5], high-order sliding mode control (HOSMC) methods [6–8] exhibit strong robustness and a reduced chattering effect while dealing with uncertainties. For example, Zhang et al. [8] proposed a quasi-continuous HOSMC for HV to effectively alleviate the chattering phenomena. In addition to the chattering effect, several design requirements also should be considered for practical HV control under the effects of uncertainties. The priority is guaranteeing the stability. Furthermore, in order to ensure a satisfactory control performance, high-accuracy tracking of trajectory commands and lower fuel consumption are desired. However, when uncertainties are involved in the nonlinear control structure of HV, it is a challenge to adjust design parameters to reach a satisfied level of tracking performance. Two problems have appeared because of introducing uncertainties into the HOSM control of HV.The first problem is that the modeling of the relation between the design parameters and the HV tracking performance under the effect of uncertain parameters is complex. Dealing with uncertainty in a probabilistic way, stochastic robustness analysis (SRA) was first proposed by Stengel and Ray [9], and it is an effective method to evaluate the extent to which the specified design requirements are satisfied. A cost function for SRA is formulated to estimate the likelihood that the design requirements are not satisfied. Subsequently, the design parameter space is searched to minimize the cost function to obtain the optimal performance in the presence of uncertainties [10]. Cao et al. [11] optimized the HV controller parameters by using SRA and hybrid PSO algorithm. However, only the dynamic response indices of step command were concerned in the cost function for SRA [11–13]. In order to achieve a desired tracking performance despite uncertainties, it is necessary to introduce appropriate indices that characterize the command tracking process and corresponding indicator functions into the optimization problem modeling of HV.The second important problem in the HOSM control of HV involves solving the optimization problem. Conventional optimization methods, such as the gradient search method, are no longer suitable given that the partial derivative of the cost function in SRA is difficult to obtain. For complex optimization problem involving uncertainties, a high efficiency computational intelligence optimization algorithm is required to determine the optimal controller parameters of HV to achieve a satisfied level of tracking performance under the influence of uncertainties. Nowadays, various computational intelligence techniques [14, 15], such as genetic algorithm (GA) [16], particle swarm optimization (PSO) [17], and differential evolutionary (DE), have been proposed for complex optimization problems with the development of computation technology.Among computational algorithms, the fireworks algorithm (FWA) is a relatively new swarm intelligence-based algorithm proposed by Tan and Zhu [18]. It simulates the process of fireworks explosion, in which the “good” fireworks generate more sparks in smaller explosion areas. Numerical experiments indicated that FWA converges to a global optimum with a smaller number of function evaluations than PSO and GA [19]. Li et al. [20] proposed an adaptive fireworks algorithm (AFWA) in which the explosion amplitude of fireworks that fails to produce a better spark increases. To improve interaction of solutions, hybrid algorithm of FWA-DE was developed by Zheng et al. [21]. Zhang et al. [22] proposed an improved FWA by enhancing fireworks interaction. With respect to improvements in the FWA [20–23], it is recognized that the diversification mechanism of FWA does not utilize more information on other qualified solutions in the swarm. Therefore, with respect to the HV control under uncertainties that are randomly distributed, it is necessary to develop an improved FWA with enhanced solutions interaction to effectively solve the complex optimization problem of searching for the optimal controller.In this study, an improved hybrid FWA-based parameter optimization method is proposed for HV control to achieve an excellent tracking performance in the presence of uncertainties. The main contributions are as follows:( 1 ) The uncertainties that are randomly distributed are considered in the modeling phase via SRA. The cost function evaluating the probability of design requirements violation is formulated to model the complex relation between design parameters and tracking performance of the uncertain HV system. Appropriate indices of the command tracking response are developed.( 2 ) A hybrid FWA to search for the optimal design parameters is proposed for the complex optimization problem involving uncertainties to satisfy design requirements with high probability. The introduction of the hybrid FWA into SRA effectively optimizes the tracking performance of the nonlinear HV system under uncertainties.This study is organized as follows: In Section2, the optimization problem in the HOSM control of HV is introduced. In Section 3, the complex relation between design parameters and HV performance under uncertainties is modeled. Section 4 proposes a new hybrid FWA to determine the optimal parameters of HV. Section 5 investigates the global convergence of the proposed hybrid FWA, and the simulation and comparison results are demonstrated. A few conclusions are made in Section 6.


## 2. HOSM Control Structure of HV with Uncertainties

The control-oriented model of a generic hypersonic vehicle (HV) is described by [24]. An inverse-square-law gravitational model and centripetal acceleration are considered, and the dynamic differential equations for velocity V, altitude h, flight-path angle γ, angle of attack α, and pitch rate q of HV are as follows:(1)V˙=Tcos⁡α-Dm-μsin⁡γr2,h˙=Vsin⁡γ,γ˙=L+Tsin⁡αmV-μ-V2rcos⁡γVr2,α˙=q-γ˙,q˙=MyyIyy,with(2)L=12ρV2SCLα,D=12ρV2SCDα,T=12ρV2SCTβ,Myy=12ρV2Sc-CMα+CMα,δe+CMα,q,r=h+RE,where L is the lift, D is the drag, T is the thrust, and Myy is the pitching moment. m, r, RE, μ, and ρ denote the mass, radial distance, radius of the Earth, gravitational constant, and density of air, respectively. Additionally, S, c-, and Iyy denote the reference area, mean aerodynamic chord, and the moment of inertia about y-body axes, respectively. δe denotes the elevator deflection, and β denotes the engine throttle setting.The thrustT in (2) is provided by the engine dynamics, and this is represented as follows [10]:(3)β¨=k1β˙+k2β+k3βc,where βc denotes the engine throttle setting command. It is adopted that k1=k2=0 and k3=1 for proper modeling of engine dynamics.In order to guarantee the robustness of the HV flight control system, the parametric uncertainties in (1)-(2) are considered as follows:(4)m=m01+Δm,ρ=ρ01+Δρ,CL=CL01+ΔCL,CD=CD01+ΔCD,CT=CT01+ΔCT,CM=CM01+ΔCM,where the uncertainties Δm, Δρ, ΔCL, ΔCD, ΔCT, and ΔCM are bounded.HV system (1) with engine dynamics is highly nonlinear. The relationship between input variables u=δe,βcT and the output variables y=V,hT is apparently expressed by the feedback linearization method [10]. We differentiate V three times and differentiate h four times, and we obtain the following expressions:(5)V˙=Tcos⁡α-Dm-μsin⁡γr2≜fx,V¨=Ω1x˙m,V⃛=Ω1x¨+x˙TΩ2x˙m,(6)h˙=Vsin⁡γ,h¨=V˙sin⁡γ+Vγ˙cos⁡γ,h⃛=V¨sin⁡γ+2V˙γ˙cos⁡γ-Vγ˙2sin⁡γ+Vγ¨cos⁡γ,(7)h4=V⃛sin⁡γ+3V¨γ˙cos⁡γ-3V˙γ˙2sin⁡γ+3V˙γ¨cos⁡γ-3Vγ˙γ¨sin⁡γ-Vγ˙3cos⁡γ+Vγ⃛cos⁡γ,(8)γ˙=L+Tsin⁡α/mV-μ-V2rcos⁡γ/Vr2≜gx,γ¨=π1x˙,γ⃛=π1x¨+x˙Tπ2x˙,where x=V,γ,α,β,hT, Ω1=∂fx/∂x, Ω2=∂Ω1/∂x, and π1=∂gx/∂x, π2=∂π1/∂x.In order to force the velocity and altitude to track the time-varying commanded outputyc=Vc,hcT, we define the velocity sliding tracking error and the altitude sliding tracking error as σv=V-Vc and σh=h-hc, respectively. Based on (5) and (7), we have(9)σ⃛vσh4=Fx+Gxu,where the formulations of F, G, Ω1, Ω2, π1, and π2 are the same as those in [10].As stated in [10], the matrix Gx in (9) is nonsingular over the entire flight envelope of HV, so (9) is decoupled with the auxiliary control input v as follows:(10)u=Gx-1-Fx+v,v=v1v2T=σ⃛vσh4T.A previous study [6] indicates that if appropriate control parameters are designed, then the finite time stabilization of system (9) is guaranteed by the quasi-continuous HOSMC v1 and v2, and this is given as follows:(11)v1=-αξvφ2,3N2,3,v2=-αξhφ3,4N3,4,with(12)φ2,3=σ¨v+βξv2σ˙v+βξv1σv2/3-1/2σ˙v+βξv1σv2/3sign⁡σv,N2,3=σ¨v+βξv2σ˙v+βξv1σv2/31/2,φ3,4=σ⃛h+βξh3σ¨h+βξh2σ˙h+βξh1σh3/42/3-1/2σ¨h+βξh2σ˙h+βξh1σh3/4-1/3σ˙h+βξh1σh3/4signσh,N3,4=σ⃛h+βξh3σ¨h+βξh2σ˙h+βξh1σh3/42/31/2.The HV control structure based on HOSM is shown in Figure1.Figure 1
HV control structure based on HOSM.For the quasi-continuous HOSM controller (11), the design parameters αξi and βξij(i=h,v,j=1,2,…) define the output trajectory of the HV system, which is shown in Figure 2. In the figure, altitude commands are in the dotted lines, and tracking trajectories are in the solid lines. In order to satisfy the design requirements, it is necessary to optimize the controller parameters.Figure 2
Altitude tracking trajectories using different design parameters.Furthermore, it is more appealing to satisfy the HV control design requirements under the effects of uncertainties. In Figure3, the tracking trajectories with uncertain parameters generated randomly are depicted by the solid lines, and the altitude commands are shown by the dotted lines. Within two dashed lines are the trajectories that meet the design requirements.Figure 3
Altitude tracking trajectories under uncertainties.The simulations indicate that the same set of design parameters will generate various trajectories in the presence of uncertainties. Therefore, it is necessary to employ a proper measure to quantify a set of data values. In this study, the probability that the design requirements are not satisfied is used for the HV performance evaluation with uncertainties.Therefore, the target of HV control design involves determining the optimal HOSM control parameters to satisfy design requirements with high probability. It is necessary to solve the following two problems in the HOSM parameter optimization:(1) to develop a cost function that evaluates the likelihood of system instability and the violation of the design requirements, so that the complex relation between HV design parameters and the performance under uncertainties is modeled; (2) to solve the complex optimization problem related to the uncertainties by a high efficient computational intelligence optimization algorithm.


## 3. Stochastic Robustness Analysis of HV

The concept of stochastic robustness was proposed by Stengel and Ray [9], and this is effective in evaluating the extent to which the specified design requirements are satisfied. We deal with uncertainties in a probabilistic way, and thus a cost function to evaluate the likelihood of system instability and the violation of design requirements is formulated via SRA.The flowchart of HOSM control design for HV based on SRA is shown in Figure4.Figure 4
Flowchart of HOSM control design based on SRA.In Figure4, a closed-loop HV system with uncertain parameters u~k is denoted by the dotted box. I· denotes the indicator function corresponding to the design requirement. The value of I· is within 0,1, and this is 0 if an acceptable performance appears and is 1 otherwise.With the indicator functionI·, the probability of satisfying a certain performance requirement P is defined by an integral of the corresponding indicator function over the expected variation space of parametric uncertainties. It is a practical method to estimate the probability P by Monte Carlo evaluation (MCE) as follows:(13)P^d=1N∑k=1NIHu~k,Cd,u~k∈Q,where H(u~k) represents the HV system with uncertain parameters u~k that are randomly selected within the parameter space Q. C(d) represents the HOSM controller with the design parameter vector d, and N denotes the sampling numbers.Thus, the cost functionJ^(d) for SRA is formed by combining the probability of various design requirements with weights as follows:(14)J^d=fP^1d,P^2d,…,where the estimated value of the cost function J^(d) approaches the true value J when the sampling number N→∞.As shown in Figure4, the optimal design parameters of HV are determined under the guidance of the cost function J^(d). Therefore, it is vital to define appropriate stochastic robustness measurements for the cost function to achieve the desired tracking performance despite uncertainties.


### 3.1. Stochastic Robustness Indices and Indicators

In this section, the stochastic robustness indices and indicators are introduced to evaluate the HV tracking performance in the presence of uncertainties.According to the requirements of HV control design, the first index is set to guarantee system stability in the presence of uncertainties. Additionally, it is necessary to develop performance indices to characterize the command tracking trajectories of HV. The tracking trajectory of a general reference signal is not standardized as that of step signal, and thus common indices, such as setting time, overshoot, and steady error, are no longer suitable. Thus, the following performance indices are introduced.(i)
Transient tracking performance:(15)TTPh=1tf∫0tfe-ahtht-hct·dt,TTPv=1tf∫0tfe-avtVt-Vct·dt,whereTTPh and TTPv represent the transient tracking performance indices for the altitude response and the velocity response, respectively. tf is the terminal time of the tracking command. ah and av are small positive constants that define the duration of the interested transient stage. A decrease in the value of TTP decreases the tracking error in the transient stage.(ii)
Steady tracking performance:(16)STPh=1tf∫0tf1-e-ahtht-hct·dt,STPv=1tf∫0tf1-e-avtVt-Vct·dt,whereSTPh and STPv represent the steady tracking performance indices for the altitude response and the velocity response, respectively. A decrease in the value of STP decreases the tracking error in the steady stage.(ii)
Fuel consumption performance:(17)FCPi=1tf∫0tfβt·dt,i=h,v,whereFCPh and FCPv represent the fuel consumption performance indices for the altitude response and the velocity response, respectively. βt denotes engine throttle setting during the flight. It is necessary to limit FCP within reasonable bounds.(iv)
Chattering effect:(18)CEi=max⁡δet,t∈ts,tf,i=h,v,wherets denotes the time when the sliding tracking errors σv and σh both tend to zero. CEh and CEv represent the maximum chatter amplitude of elevator for the altitude and velocity commands when t∈ts,tf, respectively. The chattering effect can severely deteriorate the flight control performance, and thus it is necessary to attenuate it.Through Monte Carlo sampling, the distribution of aforementioned index values is obtained from the tracking trajectories under uncertainties. In Figure5, after 200 times of random sampling, the distributions corresponding to the altitude tracking performance indices are shown. In order to evaluate the extent to which the design requirements are satisfied in the presence of uncertainties, the indicator function corresponding to the index should be carefully defined.Figure 5
Distribution of index values.The commonly used indicator is a binary function with two values of 0 and 1 to represent whether the design requirement is satisfied or not. However, for a practical engineering system, there exists an interval between the satisfied and unsatisfied performance. Thus, the following continuous function is employed as the indicatorI· as follows:(19)Ix=1x>x1x-x2x1-x22x2<x≤x10x≤x2,where x denotes the value of the performance index, such as TTPi,STPi,FCPi, and CEi, i=h,v. The positive constant x1 represents TTPi,1, STPi,1, FCPi,1, or CEi,1. The positive constant x2 represents TTPi,2, STPi,2, FCPi,2, or CEi,2, and x1 and x2 are set by the designer to define the interval between the satisfied and unsatisfied performance.


### 3.2. Optimization Problem

In order to evaluate the HV tracking performance under uncertainties, the aforementioned indices and indicator functions are employed to formulate the cost functionJ^d in (14), and they are listed in Table 1.Table 1
Stochastic robustness stability and performance indices of HV.
Metric number Weight inJd Indicator functionI⋅ Design requirements 1(2) 10.0 (10.0) I1(I2) System stability in altitude response (velocity response) 3(4) 1.0 (1.0) I 3   ( I 4 ) Transient tracking performance in altitude (velocity) response is less thanTTPh,1 (TTPv,1) 5(6) 1.0 (1.0) I 5   ( I 6 ) Steady tracking performance in altitude (velocity) response is less thanSTPh,1 (STPv,1) 7(8) 1.0 (1.0) I 7   ( I 8 ) Fuel consumption performance in altitude (velocity) response is less thanFCPh,1 (FCPv,1) 9(10) 0.5 (0.5) I 9   ( I 10 ) Chattering effect in altitude (velocity) response is less thanCEh,1 (CEv,1)By formulating the cost functionJ^d, the complex relation between the HOSM controller parameters and the HV tracking performance under uncertainties is modeled. The optimal controller parameters are obtained by solving the following optimization problem:(20)MinimizeJd=∑i=110ωiP^idsubjecttod=αξv,βξv1,βξv2,αξh,βξh1,βξh2,βξh3∈dmin,dmax,where d is the design parameter vector in the HV controller (11), and this is searched within dmin,dmax. ωi denotes the weights for the probabilities P^i(d) of various design requirements. The weight in cost function allows a trade-off between design requirements.Optimization problem (20) is a constrained nonlinear and nonconvex optimization problem, in which the cost function value is calculated with the Monte Carlo method. It is very difficult and time-consuming to determine the optimal solution.Therefore, for complex optimization problem (20) related to the uncertainties, it is necessary to develop a high efficient computational intelligence optimization algorithm to determine the optimal HV control parameters, so that an excellent tracking performance can be achieved despite uncertainties.


## 3.1. Stochastic Robustness Indices and Indicators

In this section, the stochastic robustness indices and indicators are introduced to evaluate the HV tracking performance in the presence of uncertainties.According to the requirements of HV control design, the first index is set to guarantee system stability in the presence of uncertainties. Additionally, it is necessary to develop performance indices to characterize the command tracking trajectories of HV. The tracking trajectory of a general reference signal is not standardized as that of step signal, and thus common indices, such as setting time, overshoot, and steady error, are no longer suitable. Thus, the following performance indices are introduced.(i)
Transient tracking performance:(15)TTPh=1tf∫0tfe-ahtht-hct·dt,TTPv=1tf∫0tfe-avtVt-Vct·dt,whereTTPh and TTPv represent the transient tracking performance indices for the altitude response and the velocity response, respectively. tf is the terminal time of the tracking command. ah and av are small positive constants that define the duration of the interested transient stage. A decrease in the value of TTP decreases the tracking error in the transient stage.(ii)
Steady tracking performance:(16)STPh=1tf∫0tf1-e-ahtht-hct·dt,STPv=1tf∫0tf1-e-avtVt-Vct·dt,whereSTPh and STPv represent the steady tracking performance indices for the altitude response and the velocity response, respectively. A decrease in the value of STP decreases the tracking error in the steady stage.(ii)
Fuel consumption performance:(17)FCPi=1tf∫0tfβt·dt,i=h,v,whereFCPh and FCPv represent the fuel consumption performance indices for the altitude response and the velocity response, respectively. βt denotes engine throttle setting during the flight. It is necessary to limit FCP within reasonable bounds.(iv)
Chattering effect:(18)CEi=max⁡δet,t∈ts,tf,i=h,v,wherets denotes the time when the sliding tracking errors σv and σh both tend to zero. CEh and CEv represent the maximum chatter amplitude of elevator for the altitude and velocity commands when t∈ts,tf, respectively. The chattering effect can severely deteriorate the flight control performance, and thus it is necessary to attenuate it.Through Monte Carlo sampling, the distribution of aforementioned index values is obtained from the tracking trajectories under uncertainties. In Figure5, after 200 times of random sampling, the distributions corresponding to the altitude tracking performance indices are shown. In order to evaluate the extent to which the design requirements are satisfied in the presence of uncertainties, the indicator function corresponding to the index should be carefully defined.Figure 5
Distribution of index values.The commonly used indicator is a binary function with two values of 0 and 1 to represent whether the design requirement is satisfied or not. However, for a practical engineering system, there exists an interval between the satisfied and unsatisfied performance. Thus, the following continuous function is employed as the indicatorI· as follows:(19)Ix=1x>x1x-x2x1-x22x2<x≤x10x≤x2,where x denotes the value of the performance index, such as TTPi,STPi,FCPi, and CEi, i=h,v. The positive constant x1 represents TTPi,1, STPi,1, FCPi,1, or CEi,1. The positive constant x2 represents TTPi,2, STPi,2, FCPi,2, or CEi,2, and x1 and x2 are set by the designer to define the interval between the satisfied and unsatisfied performance.


## 3.2. Optimization Problem

In order to evaluate the HV tracking performance under uncertainties, the aforementioned indices and indicator functions are employed to formulate the cost functionJ^d in (14), and they are listed in Table 1.Table 1
Stochastic robustness stability and performance indices of HV.
Metric number Weight inJd Indicator functionI⋅ Design requirements 1(2) 10.0 (10.0) I1(I2) System stability in altitude response (velocity response) 3(4) 1.0 (1.0) I 3   ( I 4 ) Transient tracking performance in altitude (velocity) response is less thanTTPh,1 (TTPv,1) 5(6) 1.0 (1.0) I 5   ( I 6 ) Steady tracking performance in altitude (velocity) response is less thanSTPh,1 (STPv,1) 7(8) 1.0 (1.0) I 7   ( I 8 ) Fuel consumption performance in altitude (velocity) response is less thanFCPh,1 (FCPv,1) 9(10) 0.5 (0.5) I 9   ( I 10 ) Chattering effect in altitude (velocity) response is less thanCEh,1 (CEv,1)By formulating the cost functionJ^d, the complex relation between the HOSM controller parameters and the HV tracking performance under uncertainties is modeled. The optimal controller parameters are obtained by solving the following optimization problem:(20)MinimizeJd=∑i=110ωiP^idsubjecttod=αξv,βξv1,βξv2,αξh,βξh1,βξh2,βξh3∈dmin,dmax,where d is the design parameter vector in the HV controller (11), and this is searched within dmin,dmax. ωi denotes the weights for the probabilities P^i(d) of various design requirements. The weight in cost function allows a trade-off between design requirements.Optimization problem (20) is a constrained nonlinear and nonconvex optimization problem, in which the cost function value is calculated with the Monte Carlo method. It is very difficult and time-consuming to determine the optimal solution.Therefore, for complex optimization problem (20) related to the uncertainties, it is necessary to develop a high efficient computational intelligence optimization algorithm to determine the optimal HV control parameters, so that an excellent tracking performance can be achieved despite uncertainties.


## 4. Optimization Technique with Improved Hybrid Fireworks Algorithm

In this section, we propose a hybrid FWA to solve the complex optimization problem of determining the optimal HV control parameters under uncertainties. First, by introducing the GA operators into the mutation process of AFWA, a hybrid FWA is developed with an improved diversification mechanism. Subsequently, the process of the hybrid FWA-based parameter optimization method is illustrated.


### 4.1. Adaptive Fireworks Algorithm

Inspired by the fireworks explosion, FWA is a relatively new swarm intelligence-based algorithm proposed by Tan and Zhu [18]. In FWA, the fireworks and sparks are considered as the potential solutions in the search space, and the explosion is viewed as a local search around the location of fireworks. The FWA converges to a global optimum with a lower number of function evaluations than those of the PSO and GA [19]. Subsequently, the AFWA [20] was developed to improve the local search capability of the best firework.The search process of AFWA is as follows:(1)
Initialization: randomly set the initial locations of fireworks.(2)
Explosion: each firework generates a set of sparks by executing the regular explosion operation.(3)
Gaussian mutation: select a few fireworks randomly, and execute the Gaussian explosion (mutation) operation on the selected fireworks to generate several sparks.(4)
Adaptive amplitude calculation: select the best individual as a firework in the next generation, and calculate its adaptive explosion amplitude.(5)
Selection: randomly select other fireworks from all individuals.(6)
Return to Step(2) until the stop criterion is fulfilled.In order to execute the regular explosion operation in Step(2), the number of sparks Ni of each firework is calculated as follows:(21)Ni=roundNc·fmax-fxi+ε∑i=1nfmax-fxi+ε.The explosion amplitudeAi is as follows:(22)Ai=Ac·fxi-fmin+ε∑i=1nfxi-fmin+ε,where n is the number of fireworks. Nc and Ac are two parameters that control the number of sparks and explosion amplitude, respectively. fxi represents the fitness value of xi, and fmax and fmin denote the maximum and minimum values of the cost function among the n fireworks, respectively. A small constant ε is to avoid zero-division error.In order to avoid the overwhelming effect of the best firework, the bound of the spark number is set as follows:(23)Ni=NminifNi<NminNmaxelseifNi>NmaxNielse,where Nmax and Nmin are the upper and lower bounds for Ni.For aD-dimension problem, after the calculation of spark number and explosion amplitude, the location of each spark xj is obtained by randomly setting approximately half of the dimensions (z dimensions), and for each dimension k, the value xjk (1≤k≤z, 1≤j≤Ni) is set based on xik (1≤i≤n). Therefore, the locations of the explosion sparks are set as follows:(24)xjk=xik+Ai·rand-1,1.In order to maintain the diversity, for a few randomly selected fireworks, approximately half of the dimensions are selected to change. The mutation sparks are generated by adding a Gaussian distribution coefficient toxik as follows:(25)xjk=xik+x∗k-xik·Gaussian0,1,where x∗k is the position ofkth dimension of the best firework x∗.If the new locations of the newly generated sparks are beyond the search space, they are mapped within the search space as follows:(26)xjk=xk,min+rand0,1·xk,max-xk,min,where xk,max and xk,min denote the upper and lower bounds of the kth dimension of the search space, respectively.In order to improve the local search capability, the best individual is selected as a firework in the next generation. It has adaptive explosion amplitude calculated by selecting an individual that satisfies the following conditions:(1) Its fitness is worse when compared with that of the best firework in the current generation. (2) Its distance to the best individual is minimal among all individuals that satisfy (1). This is expressed as follows:(27)s^=arg⁡minsi⁡dsi,s∗,fsi>fx∗,where si denotes all sparks; s∗ denotes the best individual among sparks and fireworks. x∗ is the best firework in the current generation, and d· represents the distance.The adaptive amplitude of best firework in next generation is calculated as follows:(28)A∗g+1=0.5·A∗g+λ·s^-x∗∞,where A∗g and A∗g+1 are the adaptive amplitude in current generation g and the next generation g+1, respectively. λ is a positive constant (usually higher than 1), and ·∞ represents the infinity norm.The search process indicates that the diversification mechanism of FWA does not utilize more information on all the qualified solutions, and thus it is necessary to enhance the interaction between fireworks and sparks. It is well known that GA is an efficient evolutionary algorithm that performs searches by combining possible solutions in different directions [16]. Additionally, GA exhibits potential parallelism, and thus individuals can be compared simultaneously. Therefore, we introduce GA into the mutation process of AFWA to generate more diverse and fitter solutions.


### 4.2. Hybrid Fireworks Algorithm with the Genetic Operator

In order to improve the search efficiency, the main idea in the proposed hybrid FWA involves utilizing all individuals (fireworks and sparks) to generate new individuals. In order to generate more diverse and fitter solutions, another idea involves selecting the father and mother from individuals with different features that correspond to “core individuals” and “noncore individuals.” Core individuals include the best firework and the sparks generated by the best firework. They exhibit better fitness values and closer locations. Noncore individuals include the other “bad” fireworks and sparks generated by them. They are more diverse.The process of the genetic operator is given as follows:(1)
Encoding: encode solutions to become chromosomes (individuals) with discrete units termed as genes.(2)
Recombination pool construction: construct recombination pool with qualified individuals.(3)
Parent selection: select parents from core individuals and noncore individuals, respectively.(4)
Crossover and mutation also exist.The process of the genetic operator is illustrated in Figure6.Figure 6
Process of the genetic operator.The steps in the genetic operator are stated in detail as follows.(i) Encoding. The D-dimension solutions are encoded to D-dimension chromosomes, in which each gene represents the value of corresponding dimension of a solution.(ii) Recombination Pool Construction. In order to improve the efficiency of crossover and mutation operations, two pools to select father and mother are constructed. The pool for the selection of the father is constructed by the core individuals from two sources, which include all the fathers (m fathers) in the last generation and several core individuals selected in the current generation (m core individuals). It aids in utilizing the information of the fitter individuals in a wider range. Similarly, the pool for the selection of the mother is constructed by all the mothers (m mothers) in the last generation and several noncore individuals selected in the current generation (m noncore individuals).With respect to the core individuals that are fitter and located closer, a random selection is applied among them to construct the pool for father selection. Conversely, the noncore individuals are diverse. Therefore, a roulette wheel is employed to select the fitter ones to construct the pool for the selection of the mother.The algorithm of constructing the recombination pool is shown in Algorithm1.Algorithm 1:Construct the recombination pool.
( 1 ) Construct the pool for the father selection as follows: (1.1) All the fathers of the last generation (m fathers) are reserved in the pool. (1.2) Randomly select m core individuals in the current generation to join the pool. ( 2 ) Construct the pool for the mother selection as follows: (2.1) All the mothers of the last generation (m mothers) are reserved in the pool. (2.2) Select m noncore individuals in the current generation by the roulette wheel to join the pool.(iii) Parent Selection. The individuals from the current generation are preferred to select the parents with a higher probability of generating diverse and fitter offspring. In order to select fathers from pool, the fathers of the last generation are replaced by other individuals that have better fitness, and the remaining fathers of the last generation may be replaced by other individuals again with a probability of Pr(0<Pr<1). Mothers are selected in the same way as the fathers.The algorithm of selecting parents is shown in Algorithm2.Algorithm 2:Select parents.
( 1 ) Select m fathers from the pool as follows: (1.1) Replace the last generation’s fathers by other individuals in pool that have better fitness. (1.2) Replace the remaining last generation’s fathers again with a probability Pr (0<Pr<1). ( 2 ) Select m mothers from pool as follows: (2.1) Replace the last generation’s mothers by other individuals in pool that have better fitness. (2.2) Replace the remaining last generation’s mothers again with a probability Pr.(iv) Crossover and Mutation. The selected parents are randomly paired to exchange information to generate new two individuals. In the crossover, the tails of a pair of chromosomes (individuals) are swapped at a random point along the gene sequence with a crossover probability Pc(0.7<Pc<1). After the crossover, the gene in sequence is mutated. This means the offspring are obtained by randomly setting approximately Pm of the dimensions of the individual within the search space, where Pm denotes mutation probability (0<Pm<0.2).Thus, a new hybrid FWA is proposed by introducing the GA into the mutation process of AFWA. The flowchart of the proposed optimization algorithm is shown in Figure7.Figure 7
The flowchart of hybrid FWA.Herexmax and xmin are the upper and the lower bounds of the search space, respectively.


### 4.3. Hybrid FWA-Based Parameter Optimization

The proposed hybrid FWA-based parameter optimization method combines the advantages of SRA and the hybrid FWA. By the SRA, the cost function is given to evaluate the HV tracking performance under uncertainties. Subsequently, the hybrid FWA is used to determine the optimal design parameters to satisfy the tracking performance requirements of HV with high probability. The flowchart of the proposed hybrid FWA-based parameter optimization method is given in Figure8.Figure 8
Flowchart of the hybrid FWA-based parameter optimization method.To illustrate the search process of the proposed hybrid FWA-based parameter optimization method in detail, the following steps are given:(1)
Generate several solutions by the hybrid FWA search process.
(a)
Randomly initialize a population ofn fireworks in the search space.(b)
For each firework, generateNi explosion sparks within the explosion amplitude Ai, and subsequently the positions xj of the explosion sparks are obtained.(c)
Encode all individuals as chromosomes.(d)
Select2m parents from all the chromosomes, and 2m diverse individuals are generated via the genetic operator.(2)
Evaluate the solution’s fitness with the cost function in SRA.
(a)
With the stochastic robustness indices listed in Table1 and the indicator function I· as defined in (19), calculate the indicator function value for the corresponding index.(b)
By the Monte Carlo simulation,N samples under uncertainties are generated to estimate the probability P^d in which the design requirements of HV control are not satisfied.(c)
For all the solutions generated in the search process, calculate the cost functionJ^d in (20).(3)
Prepare for the next step searching.
(a)
After the evaluation of all solution’s fitness, the optimal solution is selected as a firework in the next generation. Its adaptive amplitudeA∗ is calculated based on (28).(b)
Randomly selectn-1 fireworks among all the individuals.(4)
Check if the stop criterion is fulfilled.(5)
The optimal HOSM parameters are obtained, and an excellent HV tracking performance under uncertainties is achieved.


## 4.1. Adaptive Fireworks Algorithm

Inspired by the fireworks explosion, FWA is a relatively new swarm intelligence-based algorithm proposed by Tan and Zhu [18]. In FWA, the fireworks and sparks are considered as the potential solutions in the search space, and the explosion is viewed as a local search around the location of fireworks. The FWA converges to a global optimum with a lower number of function evaluations than those of the PSO and GA [19]. Subsequently, the AFWA [20] was developed to improve the local search capability of the best firework.The search process of AFWA is as follows:(1)
Initialization: randomly set the initial locations of fireworks.(2)
Explosion: each firework generates a set of sparks by executing the regular explosion operation.(3)
Gaussian mutation: select a few fireworks randomly, and execute the Gaussian explosion (mutation) operation on the selected fireworks to generate several sparks.(4)
Adaptive amplitude calculation: select the best individual as a firework in the next generation, and calculate its adaptive explosion amplitude.(5)
Selection: randomly select other fireworks from all individuals.(6)
Return to Step(2) until the stop criterion is fulfilled.In order to execute the regular explosion operation in Step(2), the number of sparks Ni of each firework is calculated as follows:(21)Ni=roundNc·fmax-fxi+ε∑i=1nfmax-fxi+ε.The explosion amplitudeAi is as follows:(22)Ai=Ac·fxi-fmin+ε∑i=1nfxi-fmin+ε,where n is the number of fireworks. Nc and Ac are two parameters that control the number of sparks and explosion amplitude, respectively. fxi represents the fitness value of xi, and fmax and fmin denote the maximum and minimum values of the cost function among the n fireworks, respectively. A small constant ε is to avoid zero-division error.In order to avoid the overwhelming effect of the best firework, the bound of the spark number is set as follows:(23)Ni=NminifNi<NminNmaxelseifNi>NmaxNielse,where Nmax and Nmin are the upper and lower bounds for Ni.For aD-dimension problem, after the calculation of spark number and explosion amplitude, the location of each spark xj is obtained by randomly setting approximately half of the dimensions (z dimensions), and for each dimension k, the value xjk (1≤k≤z, 1≤j≤Ni) is set based on xik (1≤i≤n). Therefore, the locations of the explosion sparks are set as follows:(24)xjk=xik+Ai·rand-1,1.In order to maintain the diversity, for a few randomly selected fireworks, approximately half of the dimensions are selected to change. The mutation sparks are generated by adding a Gaussian distribution coefficient toxik as follows:(25)xjk=xik+x∗k-xik·Gaussian0,1,where x∗k is the position ofkth dimension of the best firework x∗.If the new locations of the newly generated sparks are beyond the search space, they are mapped within the search space as follows:(26)xjk=xk,min+rand0,1·xk,max-xk,min,where xk,max and xk,min denote the upper and lower bounds of the kth dimension of the search space, respectively.In order to improve the local search capability, the best individual is selected as a firework in the next generation. It has adaptive explosion amplitude calculated by selecting an individual that satisfies the following conditions:(1) Its fitness is worse when compared with that of the best firework in the current generation. (2) Its distance to the best individual is minimal among all individuals that satisfy (1). This is expressed as follows:(27)s^=arg⁡minsi⁡dsi,s∗,fsi>fx∗,where si denotes all sparks; s∗ denotes the best individual among sparks and fireworks. x∗ is the best firework in the current generation, and d· represents the distance.The adaptive amplitude of best firework in next generation is calculated as follows:(28)A∗g+1=0.5·A∗g+λ·s^-x∗∞,where A∗g and A∗g+1 are the adaptive amplitude in current generation g and the next generation g+1, respectively. λ is a positive constant (usually higher than 1), and ·∞ represents the infinity norm.The search process indicates that the diversification mechanism of FWA does not utilize more information on all the qualified solutions, and thus it is necessary to enhance the interaction between fireworks and sparks. It is well known that GA is an efficient evolutionary algorithm that performs searches by combining possible solutions in different directions [16]. Additionally, GA exhibits potential parallelism, and thus individuals can be compared simultaneously. Therefore, we introduce GA into the mutation process of AFWA to generate more diverse and fitter solutions.


## 4.2. Hybrid Fireworks Algorithm with the Genetic Operator

In order to improve the search efficiency, the main idea in the proposed hybrid FWA involves utilizing all individuals (fireworks and sparks) to generate new individuals. In order to generate more diverse and fitter solutions, another idea involves selecting the father and mother from individuals with different features that correspond to “core individuals” and “noncore individuals.” Core individuals include the best firework and the sparks generated by the best firework. They exhibit better fitness values and closer locations. Noncore individuals include the other “bad” fireworks and sparks generated by them. They are more diverse.The process of the genetic operator is given as follows:(1)
Encoding: encode solutions to become chromosomes (individuals) with discrete units termed as genes.(2)
Recombination pool construction: construct recombination pool with qualified individuals.(3)
Parent selection: select parents from core individuals and noncore individuals, respectively.(4)
Crossover and mutation also exist.The process of the genetic operator is illustrated in Figure6.Figure 6
Process of the genetic operator.The steps in the genetic operator are stated in detail as follows.(i) Encoding. The D-dimension solutions are encoded to D-dimension chromosomes, in which each gene represents the value of corresponding dimension of a solution.(ii) Recombination Pool Construction. In order to improve the efficiency of crossover and mutation operations, two pools to select father and mother are constructed. The pool for the selection of the father is constructed by the core individuals from two sources, which include all the fathers (m fathers) in the last generation and several core individuals selected in the current generation (m core individuals). It aids in utilizing the information of the fitter individuals in a wider range. Similarly, the pool for the selection of the mother is constructed by all the mothers (m mothers) in the last generation and several noncore individuals selected in the current generation (m noncore individuals).With respect to the core individuals that are fitter and located closer, a random selection is applied among them to construct the pool for father selection. Conversely, the noncore individuals are diverse. Therefore, a roulette wheel is employed to select the fitter ones to construct the pool for the selection of the mother.The algorithm of constructing the recombination pool is shown in Algorithm1.Algorithm 1:Construct the recombination pool.
( 1 ) Construct the pool for the father selection as follows: (1.1) All the fathers of the last generation (m fathers) are reserved in the pool. (1.2) Randomly select m core individuals in the current generation to join the pool. ( 2 ) Construct the pool for the mother selection as follows: (2.1) All the mothers of the last generation (m mothers) are reserved in the pool. (2.2) Select m noncore individuals in the current generation by the roulette wheel to join the pool.(iii) Parent Selection. The individuals from the current generation are preferred to select the parents with a higher probability of generating diverse and fitter offspring. In order to select fathers from pool, the fathers of the last generation are replaced by other individuals that have better fitness, and the remaining fathers of the last generation may be replaced by other individuals again with a probability of Pr(0<Pr<1). Mothers are selected in the same way as the fathers.The algorithm of selecting parents is shown in Algorithm2.Algorithm 2:Select parents.
( 1 ) Select m fathers from the pool as follows: (1.1) Replace the last generation’s fathers by other individuals in pool that have better fitness. (1.2) Replace the remaining last generation’s fathers again with a probability Pr (0<Pr<1). ( 2 ) Select m mothers from pool as follows: (2.1) Replace the last generation’s mothers by other individuals in pool that have better fitness. (2.2) Replace the remaining last generation’s mothers again with a probability Pr.(iv) Crossover and Mutation. The selected parents are randomly paired to exchange information to generate new two individuals. In the crossover, the tails of a pair of chromosomes (individuals) are swapped at a random point along the gene sequence with a crossover probability Pc(0.7<Pc<1). After the crossover, the gene in sequence is mutated. This means the offspring are obtained by randomly setting approximately Pm of the dimensions of the individual within the search space, where Pm denotes mutation probability (0<Pm<0.2).Thus, a new hybrid FWA is proposed by introducing the GA into the mutation process of AFWA. The flowchart of the proposed optimization algorithm is shown in Figure7.Figure 7
The flowchart of hybrid FWA.Herexmax and xmin are the upper and the lower bounds of the search space, respectively.


## 4.3. Hybrid FWA-Based Parameter Optimization

The proposed hybrid FWA-based parameter optimization method combines the advantages of SRA and the hybrid FWA. By the SRA, the cost function is given to evaluate the HV tracking performance under uncertainties. Subsequently, the hybrid FWA is used to determine the optimal design parameters to satisfy the tracking performance requirements of HV with high probability. The flowchart of the proposed hybrid FWA-based parameter optimization method is given in Figure8.Figure 8
Flowchart of the hybrid FWA-based parameter optimization method.To illustrate the search process of the proposed hybrid FWA-based parameter optimization method in detail, the following steps are given:(1)
Generate several solutions by the hybrid FWA search process.
(a)
Randomly initialize a population ofn fireworks in the search space.(b)
For each firework, generateNi explosion sparks within the explosion amplitude Ai, and subsequently the positions xj of the explosion sparks are obtained.(c)
Encode all individuals as chromosomes.(d)
Select2m parents from all the chromosomes, and 2m diverse individuals are generated via the genetic operator.(2)
Evaluate the solution’s fitness with the cost function in SRA.
(a)
With the stochastic robustness indices listed in Table1 and the indicator function I· as defined in (19), calculate the indicator function value for the corresponding index.(b)
By the Monte Carlo simulation,N samples under uncertainties are generated to estimate the probability P^d in which the design requirements of HV control are not satisfied.(c)
For all the solutions generated in the search process, calculate the cost functionJ^d in (20).(3)
Prepare for the next step searching.
(a)
After the evaluation of all solution’s fitness, the optimal solution is selected as a firework in the next generation. Its adaptive amplitudeA∗ is calculated based on (28).(b)
Randomly selectn-1 fireworks among all the individuals.(4)
Check if the stop criterion is fulfilled.(5)
The optimal HOSM parameters are obtained, and an excellent HV tracking performance under uncertainties is achieved.


## 5. Simulation Study

### 5.1. Computational Intelligence Algorithm Test Cases

In this section, typical nonlinear benchmark functions in [25] are employed to test the effectiveness of the proposed hybrid FWA. For the comparison, the GA, PSO, AFWA, and proposed hybrid FWA are run on the benchmarks for 300000 evaluations per function. Each experiment for testing algorithm is repeated 50 times.In the testing, the parameters settings of algorithms are listed in Table2.Table 2
Parameters settings in algorithms.
Algorithms Algorithm coefficients GA Population size: 200. Binary coded chromosome length: 10. Crossover probability: 0.7. Mutation probability: 0.015. PSO Particle number: 30. Inertia weight:ωmax=0.9ωmin=0.4. Learning factor: c1=c2=1.5. AFWA Total sparks number: 200. Other parameters are the same as in [20]. Hybrid FWA Genetic operator parameters:m=15, Pr=0.5, Pc=0.95, and Pm=0.025. Other parameters are the same as AFWA.The first function is the Bent Cigar function and is described as follows:(29)f1x=x12+106∑i=2Dxi2,where xi∈-100,100, D=4. The Bent Cigar function is a unimodal function and is smooth. However, it has a narrow ridge. It has the global minimum f1∗=0 when xi=0, i=1,…,D. The second function is the Rosenbrock function that is described as follows:(30)f2x=∑i=1D-1100xi2-xi+12+xi-12,where xi∈-100,100, D=30. The Rosenbrock function is a nonconvex function in which the global minimum is inside a long, narrow, and parabolic shaped flat valley. It has the global minimum f2∗=0 when xi=1, i=1,…,D. The third function is the Griewank function described as follows:(31)f3x=∑i=1Dxi24000-∏i=1Dcos⁡xii+1,where xi∈-100,100, D=30. The Griewank function is a multimodal function. It has the global minimum f3∗=0 when xi=0, i=1,…,D. The fourth function is the Alpine function described as follows:(32)f4x=∑i=1D110xisin⁡xi10+xi,where xi∈-100,100, D=30. The Alpine function is a multimodal function. It has the global minimum f4∗=0 when xi=0, i=1,…,D. The fifth function is the Rastrigin function that is described in (33)f5x=∑i=1Dxi2-10cos⁡2πxi+10,where xi∈-100,100, D=30. The Rastrigin function is a multimodal function, which has huge number of local optima. It has the global minimum f5∗=0 when xi=0, i=1,…,D. The last function is the expanded Schaffer F6 function described in (34)gx,y=0.5+sin2x2+y2-0.51+0.001x2+y22,f6x=gx1,x2+⋯+gxD-1,xD+gxD,x1,where xi∈-100,100, D=30. The expanded Schaffer F6 function is a multimodal function. It has the global minimum f6∗=0 when xi=0, i=1,…,D.The testing results are given in Table3.Table 3
Algorithm testing results.
Function ID Metric GA PSO AFWA Hybrid FWA 1 Mean 2.0817 e + 4 63.8216 44.4774 36.6983 Std. 1.3336 e + 4 59.7943 40.1607 39.9258 2 Mean 234.6298 628.8852 58.3773 56.1416 Std. 123.0631 477.8930 31.7385 23.7506 3 Mean 1.2007 1.7338 0.7527 0.7125 Std. 0.2004 0.5626 0.2188 0.1904 4 Mean 0.2287 0.6958 0.0708 0.0658 Std. 0.0687 0.4486 0.0259 0.0250 5 Mean 22.9280 55.2201 19.6178 16.8546 Std. 24.7651 24.7651 15.2610 12.0650 6 Mean 4.8939 3.7720 3.8627 3.6126 Std. 0.9717 1.2440 1.0428 0.8973As shown in Table3, the proposed hybrid FWA presents the means closest to global minimum. Therefore, the testing results indicate that the hybrid FWA proposed in this study exhibits better search efficiency, when compared to the GA, PSO, and AFWA.


### 5.2. Algorithm Analysis in Parameter Optimization

In order to analyze the parameter searching efficiency of algorithms, the GA, PSO, AFWA, and proposed hybrid FWA are used to search for the optimal design parameters of the HOSM controller of HV. In the search, for all the algorithms, the number of individuals is 32, and the number of iterations is 15. For the AFWA, the number of fireworks is 5, the total number of sparks is 32, and the number of mutation sparks is 4. For the hybrid FWA, the number of fireworks and total sparks is the same as in AFWA,m=3, and Pm=0.1. The other parameters of algorithms are set the same as shown in Table 2.The ranges of the uncertainties in HV are as follows:(35)Δm≤0.5%,Δρ≤10%,ΔCL≤10%,ΔCD≤10%,ΔCT≤10%,ΔCM≤10%.The search space of the HV controller parameters is given in Table4.Table 4
Search space of HV controller parameters.
Controller parameters Bound α ξ v [ 5,50 ] β ξ v 1 [ 0.01,10 ] β ξ v 2 [ 0.01,10 ] α ξ h [ 10,100 ] β ξ h 1 [ 0.01,10 ] β ξ h 2 [ 0.01,10 ] β ξ h 3 [ 0.01,10 ]As given in Table1, the cost function Jd in SRA is a weighted sum of 10 probabilities of the design requirements to guide the search of the HOSM controller parameters. The parameters specified for the indicator function I· are as follows: TTPh,1=5, TTPh,2=1, TTPv,1=0.05, TTPv,2=0.02, STPh,1=3, STPh,2=1, STPv,1=0.05, STPv,2=0.005, FCPh,1=1, FCPh,2=0.2, FCPv,1=1, FCPv,2=0.2, CEh,1=2, CEh,2=1, CEv,1=2, and CEv,2=1. The duration of interested transient stage is defined by the parameters ah=0.15 and av=0.4.The results of the HV performance optimization using various optimization algorithms are shown in Figure9. The x-axis of the figure shows the number of iterations, and the y-axis shows the optimal value of the cost function J^d. The comparative result indicates that the proposed hybrid FWA exhibits better global search ability for the optimal HV control parameters than that of the GA, PSO, and AFWA.Figure 9
HV performance optimization using various algorithms.


### 5.3. Results of Optimal HOSM Controller Design

With the proposed hybrid FWA-based parameter optimization method, we shall examine the performance of the optimal HOSM controller in the trajectory tracking of HV. Initially, the cruising flight conditions are as follows: Mach numberMa=15, V=15060 ft/s, h=110000 ft, γ=0deg, and q=0 deg/s. At the cruising flight conditions, the aerodynamic parameters CL, CD, CT, CMα, CMα,δe, and CMα,q are given as follows:(36)CL=0.6203α,CD=0.6450α2+0.0043378α+0.003772,CT=0.02576β,β≤1,0.02241+0.15β,β>1,CMα=-0.035α2+0.036617α+5.3261×10-6,CMα,δe=ceδe-α,ce=0.0292,CMα,q=c-q2V-6.796α2+0.3015α-0.2289.After 15 search iterations by the proposed hybrid FWA-based parameter optimization algorithm, the optimal quasi-continuous HOSM controller parameters are determined as follows:dH∗=28.50,1.74,2.90,56.66,0.39,1.22,6.85. Using AFWA, the optimal controller parameters are determined as follows: dA∗={40.65,1.36,8,35,0.40,1.15,4}. For comparison purposes, the other two sets of design parameters are given: The quasi-continuous HOSM controller parameters (not optimized) in [8] are as follows: dN={30,1,2,10,0.5,1,3}. The HOSM controller parameters determined by the improved PSO in [11] are as follows: dP∗={107.33,0.67,2.83,22.29,125.38, 0.54,1.06,5.50,26.64}.In order to demonstrate the tracking performance of HV under uncertainties, the command tracking trajectories using four sets of controller parameters are given in Figure10. In the simulation, the reference command is generated to control the HV to climb 800 ft at constant velocity in about 15 s. The parametric uncertainties are set as follows: Δm=0.1%, Δρ=2%, ΔCL=2%, ΔCD=2%, ΔCT=2%, and ΔCM=2%, which are within the range given in (35).Figure 10
Tracking trajectories using various controller parameters.In Figure10, the trajectories of altitude h, velocity V, angle of attack α, and throttle setting β are depicted by the solid lines, and the reference command is shown by the dotted line. The simulation results demonstrate that the optimal controller parameters dH∗ determined by hybrid FWA provide a stable and high-accuracy tracking of the reference command in the presence of uncertainties. The command tracking error of the HV control system using the parameters dH∗ remains the smallest, when compared to the controller parameters dA∗, dP∗, and dN. Besides, a faster dynamic response as well as lower fuel consumption is achieved using the parameters determined by the proposed hybrid FWA.Next, with randomly generated uncertainties, the command tracking trajectories using four sets of controller parameters are demonstrated in Figure11. The uncertain parameters are assumed to be uniformly distributed within the bounds given in (35). The results indicate that the optimal controller parameters dH∗ determined by the proposed hybrid FWA not only guarantee the HV system stability, but also exhibit a better tracking performance under bounded uncertainties.Figure 11
Random responses using various controller parameters.Therefore, the simulation results demonstrate that the HV controller designed by the proposed hybrid FWA-based parameter optimization method achieves an excellent tracking performance in the presence of uncertainties.


## 5.1. Computational Intelligence Algorithm Test Cases

In this section, typical nonlinear benchmark functions in [25] are employed to test the effectiveness of the proposed hybrid FWA. For the comparison, the GA, PSO, AFWA, and proposed hybrid FWA are run on the benchmarks for 300000 evaluations per function. Each experiment for testing algorithm is repeated 50 times.In the testing, the parameters settings of algorithms are listed in Table2.Table 2
Parameters settings in algorithms.
Algorithms Algorithm coefficients GA Population size: 200. Binary coded chromosome length: 10. Crossover probability: 0.7. Mutation probability: 0.015. PSO Particle number: 30. Inertia weight:ωmax=0.9ωmin=0.4. Learning factor: c1=c2=1.5. AFWA Total sparks number: 200. Other parameters are the same as in [20]. Hybrid FWA Genetic operator parameters:m=15, Pr=0.5, Pc=0.95, and Pm=0.025. Other parameters are the same as AFWA.The first function is the Bent Cigar function and is described as follows:(29)f1x=x12+106∑i=2Dxi2,where xi∈-100,100, D=4. The Bent Cigar function is a unimodal function and is smooth. However, it has a narrow ridge. It has the global minimum f1∗=0 when xi=0, i=1,…,D. The second function is the Rosenbrock function that is described as follows:(30)f2x=∑i=1D-1100xi2-xi+12+xi-12,where xi∈-100,100, D=30. The Rosenbrock function is a nonconvex function in which the global minimum is inside a long, narrow, and parabolic shaped flat valley. It has the global minimum f2∗=0 when xi=1, i=1,…,D. The third function is the Griewank function described as follows:(31)f3x=∑i=1Dxi24000-∏i=1Dcos⁡xii+1,where xi∈-100,100, D=30. The Griewank function is a multimodal function. It has the global minimum f3∗=0 when xi=0, i=1,…,D. The fourth function is the Alpine function described as follows:(32)f4x=∑i=1D110xisin⁡xi10+xi,where xi∈-100,100, D=30. The Alpine function is a multimodal function. It has the global minimum f4∗=0 when xi=0, i=1,…,D. The fifth function is the Rastrigin function that is described in (33)f5x=∑i=1Dxi2-10cos⁡2πxi+10,where xi∈-100,100, D=30. The Rastrigin function is a multimodal function, which has huge number of local optima. It has the global minimum f5∗=0 when xi=0, i=1,…,D. The last function is the expanded Schaffer F6 function described in (34)gx,y=0.5+sin2x2+y2-0.51+0.001x2+y22,f6x=gx1,x2+⋯+gxD-1,xD+gxD,x1,where xi∈-100,100, D=30. The expanded Schaffer F6 function is a multimodal function. It has the global minimum f6∗=0 when xi=0, i=1,…,D.The testing results are given in Table3.Table 3
Algorithm testing results.
Function ID Metric GA PSO AFWA Hybrid FWA 1 Mean 2.0817 e + 4 63.8216 44.4774 36.6983 Std. 1.3336 e + 4 59.7943 40.1607 39.9258 2 Mean 234.6298 628.8852 58.3773 56.1416 Std. 123.0631 477.8930 31.7385 23.7506 3 Mean 1.2007 1.7338 0.7527 0.7125 Std. 0.2004 0.5626 0.2188 0.1904 4 Mean 0.2287 0.6958 0.0708 0.0658 Std. 0.0687 0.4486 0.0259 0.0250 5 Mean 22.9280 55.2201 19.6178 16.8546 Std. 24.7651 24.7651 15.2610 12.0650 6 Mean 4.8939 3.7720 3.8627 3.6126 Std. 0.9717 1.2440 1.0428 0.8973As shown in Table3, the proposed hybrid FWA presents the means closest to global minimum. Therefore, the testing results indicate that the hybrid FWA proposed in this study exhibits better search efficiency, when compared to the GA, PSO, and AFWA.


## 5.2. Algorithm Analysis in Parameter Optimization

In order to analyze the parameter searching efficiency of algorithms, the GA, PSO, AFWA, and proposed hybrid FWA are used to search for the optimal design parameters of the HOSM controller of HV. In the search, for all the algorithms, the number of individuals is 32, and the number of iterations is 15. For the AFWA, the number of fireworks is 5, the total number of sparks is 32, and the number of mutation sparks is 4. For the hybrid FWA, the number of fireworks and total sparks is the same as in AFWA,m=3, and Pm=0.1. The other parameters of algorithms are set the same as shown in Table 2.The ranges of the uncertainties in HV are as follows:(35)Δm≤0.5%,Δρ≤10%,ΔCL≤10%,ΔCD≤10%,ΔCT≤10%,ΔCM≤10%.The search space of the HV controller parameters is given in Table4.Table 4
Search space of HV controller parameters.
Controller parameters Bound α ξ v [ 5,50 ] β ξ v 1 [ 0.01,10 ] β ξ v 2 [ 0.01,10 ] α ξ h [ 10,100 ] β ξ h 1 [ 0.01,10 ] β ξ h 2 [ 0.01,10 ] β ξ h 3 [ 0.01,10 ]As given in Table1, the cost function Jd in SRA is a weighted sum of 10 probabilities of the design requirements to guide the search of the HOSM controller parameters. The parameters specified for the indicator function I· are as follows: TTPh,1=5, TTPh,2=1, TTPv,1=0.05, TTPv,2=0.02, STPh,1=3, STPh,2=1, STPv,1=0.05, STPv,2=0.005, FCPh,1=1, FCPh,2=0.2, FCPv,1=1, FCPv,2=0.2, CEh,1=2, CEh,2=1, CEv,1=2, and CEv,2=1. The duration of interested transient stage is defined by the parameters ah=0.15 and av=0.4.The results of the HV performance optimization using various optimization algorithms are shown in Figure9. The x-axis of the figure shows the number of iterations, and the y-axis shows the optimal value of the cost function J^d. The comparative result indicates that the proposed hybrid FWA exhibits better global search ability for the optimal HV control parameters than that of the GA, PSO, and AFWA.Figure 9
HV performance optimization using various algorithms.


## 5.3. Results of Optimal HOSM Controller Design

With the proposed hybrid FWA-based parameter optimization method, we shall examine the performance of the optimal HOSM controller in the trajectory tracking of HV. Initially, the cruising flight conditions are as follows: Mach numberMa=15, V=15060 ft/s, h=110000 ft, γ=0deg, and q=0 deg/s. At the cruising flight conditions, the aerodynamic parameters CL, CD, CT, CMα, CMα,δe, and CMα,q are given as follows:(36)CL=0.6203α,CD=0.6450α2+0.0043378α+0.003772,CT=0.02576β,β≤1,0.02241+0.15β,β>1,CMα=-0.035α2+0.036617α+5.3261×10-6,CMα,δe=ceδe-α,ce=0.0292,CMα,q=c-q2V-6.796α2+0.3015α-0.2289.After 15 search iterations by the proposed hybrid FWA-based parameter optimization algorithm, the optimal quasi-continuous HOSM controller parameters are determined as follows:dH∗=28.50,1.74,2.90,56.66,0.39,1.22,6.85. Using AFWA, the optimal controller parameters are determined as follows: dA∗={40.65,1.36,8,35,0.40,1.15,4}. For comparison purposes, the other two sets of design parameters are given: The quasi-continuous HOSM controller parameters (not optimized) in [8] are as follows: dN={30,1,2,10,0.5,1,3}. The HOSM controller parameters determined by the improved PSO in [11] are as follows: dP∗={107.33,0.67,2.83,22.29,125.38, 0.54,1.06,5.50,26.64}.In order to demonstrate the tracking performance of HV under uncertainties, the command tracking trajectories using four sets of controller parameters are given in Figure10. In the simulation, the reference command is generated to control the HV to climb 800 ft at constant velocity in about 15 s. The parametric uncertainties are set as follows: Δm=0.1%, Δρ=2%, ΔCL=2%, ΔCD=2%, ΔCT=2%, and ΔCM=2%, which are within the range given in (35).Figure 10
Tracking trajectories using various controller parameters.In Figure10, the trajectories of altitude h, velocity V, angle of attack α, and throttle setting β are depicted by the solid lines, and the reference command is shown by the dotted line. The simulation results demonstrate that the optimal controller parameters dH∗ determined by hybrid FWA provide a stable and high-accuracy tracking of the reference command in the presence of uncertainties. The command tracking error of the HV control system using the parameters dH∗ remains the smallest, when compared to the controller parameters dA∗, dP∗, and dN. Besides, a faster dynamic response as well as lower fuel consumption is achieved using the parameters determined by the proposed hybrid FWA.Next, with randomly generated uncertainties, the command tracking trajectories using four sets of controller parameters are demonstrated in Figure11. The uncertain parameters are assumed to be uniformly distributed within the bounds given in (35). The results indicate that the optimal controller parameters dH∗ determined by the proposed hybrid FWA not only guarantee the HV system stability, but also exhibit a better tracking performance under bounded uncertainties.Figure 11
Random responses using various controller parameters.Therefore, the simulation results demonstrate that the HV controller designed by the proposed hybrid FWA-based parameter optimization method achieves an excellent tracking performance in the presence of uncertainties.


## 6. Conclusion

In this study, we propose an improved hybrid FWA-based parameter optimization method for nonlinear HV control under uncertainties. An approach of searching for the optimal design parameters is developed by using two processes. The first process involves modeling the relation between the design parameters and the cost function that evaluates the likelihood of system instability and design requirement violation by using SRA. Subsequently, the cost function is minimized by the improved hybrid FWA to achieve a satisfactory tracking performance for the HV system with uncertainties. The proposed method makes it easier and more efficient to solve the optimization problem of satisfying all the HV design requirements with high probability. When compared with other algorithms, the hybrid FWA exhibits better efficiency in solving the HV parameter optimization problem with respect to uncertainties. Moreover, it is also efficient in solving other complex optimization problems.



---

*Source: 9098151-2018-03-04.xml*