# Study on the Application Scheme of Aerodynamic Coefficient Identification Based on the Differential Evolution Algorithm

Authors: Pengfei Liu; Hongsong Cao; Shunshan Feng; Hengzhu Liu; Ye Du
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The present paper studies the application schemes of aerodynamic coefficient identification based on the differential evolution (DE) algorithm from free flight data. Two identification schemes utilizing the DE algorithm are proposed, including the whole discrete point (WDP) scheme and the differential interval linear (DIL) scheme. Several comparative tests are conducted to study the performances of the two schemes by using the noisy simulated flight data and the measured radar data. The results show satisfactory performance in estimating the aerodynamic coefficients for the two schemes. The WDP scheme provides more accurate and robust coefficients estimates in exchange of more than two times greater computational effort compared with the DIL scheme.


---

## Body

## 1. Introduction

The identification of aerodynamic coefficients, based on free flight data, is one of the important technologies in the field of projectile aerodynamics [1–5]. Compared with wind tunnel tests or numerical methods, aerodynamic coefficient identification can use actual flight data and full-scale projectiles. This could ensure higher fidelity and enable it to become a critical approach for determining the aerodynamic coefficients, especially in the verification phase of aerodynamic configuration design, ballistic flight control design, and fire table generation [1, 2, 6–11].In recent years, with the improvements in computer performance, modern intelligence optimization algorithms have been developed and applied to identify aerodynamic coefficients. Compared with the traditional gradient-based identification method [12–21], intelligence optimization algorithms have extreme superiority in obtaining the global optimum and handling discontinuous nonconvex objectives [22]. Due to the great value in academics and engineering, a lot of work has been devoted to the research of intelligence optimization algorithms for aerodynamic coefficient identification. The genetic algorithm (GA) was used to estimate projectile aerodynamic coefficients from free flight range data [21]. The adaptive GA was used for nonlinear aerodynamic parameter identification of an axis-symmetrical tactical missile [23]. The partial swarm optimization (PSO) algorithm was applied to the aerodynamic coefficient identification of a spinning symmetric projectile [24]. The artificial bee colony (ABC) algorithm was applied to identify the high angle of attack aerodynamic parameters [25].The differential evolution (DE) algorithm proposed by Storn and Price is a differential operator-based evolutionary optimization technique [26–28]. This algorithm requires very few control parameters and adopts real number coding to add diversity to the difference distributions [27–29]. These unique features make it accurate, robust, easy to use, and converge fast compared with previous intelligence optimization algorithms. The DE algorithm has already been successfully used in the fields of aerodynamic design, mechanical design, and digital filter design [26, 30–32].In the aerodynamic coefficient identification from the whole trajectory data, it is well known that a key process is to split the measurements into multiple small intervals, and then, the coefficients can be identified in the intervals. In the previous research [1, 2, 12, 13, 15–19, 21, 23–25], the coefficients in an interval are usually modeled as constants. In the work of Zaikang and Lyster, the drag coefficients in each interval are modeled as a cubic spline [17]. The identification results are in good agreement with the measured data. However, due to the limited measurements of the small interval and the constant assumption of the coefficient, the unbalanced noise distribution would lead to overfitting of the data, and the coupling effect between the data intervals would be ignored. These drawbacks reduced the accuracy and robustness.In this article, two identification schemes based on the DE algorithm are proposed to improve the performance of identifying the aerodynamic coefficients from trajectory data, including the whole discrete point (WDP) scheme and the differential interval linear (DIL) scheme. Several comparative tests are conducted, and the characteristics of these two schemes with regard to the computational time, accuracy, and robustness are discussed.


## 2. Projectile Flight Dynamic Model

In this study, the system model of aerodynamic coefficient identification is a standard six-degree-of-freedom (6-DOF) model typically used in flight dynamic modeling of projectiles. Figure1 gives two standard reference frames used to develop this model. The 6-DOF model consists of the inertial position components of the projectile mass center x,y,z from an inertial frame and the standard aerospace sequence Euler angles γ,ϕ2,ϕa.Figure 1
Schematic of the projectile’s position coordinates.The kinematic and dynamic differential equations for describing the motion of the projectile in flight are given as follows [1, 33]:(1)x˙y˙z˙=cosφ2cosφa−sinφa−sinφ2cosφacosφ2sinφacosφa−sinφ2sinφasinφ20cosφ2uvw,(2)γ˙φ2˙φa˙=10−tanφ20−10001cosφ2pqr,(3)u˙v˙w˙=XmYmZm+0−rqr0−rtanφ2−qrtanφ20uvw,(4)p˙q˙r˙=IA000IC000IC−1LMN+0−rqr0−rtanφ2−qrtanφ20IApICqICr.The forces acting on the projectile in (3) consist of gravity G and the aerodynamic force. The aerodynamic force is split into a standard aerodynamic force A and Magnus aerodynamic force M. Gliding flight is assumed in this study. The combination of forces is expressed in the following equation:(5)XYZ=XGYGZG+XAYAZA+XMYMZM,(6)XGYGZG=−mgcosφ2sinφacosφa−sinφ2sinφa,(7)XAYAZA=12ρVCG2S−CD−cosαcosβsinαcosαcosβ+Clαsin2β+cos2βsin2α−sinβcosαcosβ−cos2βsinαcosα,(8)XMYMZM=12ρVCG2SCypapdVCG0−sinαcosβsinβ.In (7), the drag coefficient CD is estimated by considering the zero-yaw drag coefficient CD0 and the yaw drag coefficient CD2 [2]:(9)CD=CD0+CD2=CD0+Clα+0.5CD0α2+β2.It is known thatCD has the most significant effect on the trajectory, and CD0 represents the main part of CD [34].The moment acting on the projectile in (4) consists of the moment due to A, the moment due to M, and the unsteady aerodynamic moments UA, as shown in the following equation:(10)LMN=LAMANA+LMMMNN+LUAMUANUA.The moments in (10) are expressed as follows:(11)LAMANA=12ρSlVCG2Cmα0sinαcosβ−sinβ,(12)LMMMNN=12ρSlVCG2−CnpαpdVCG0−sinβ−sinαcosβ,(13)LUAMUANUA=12ρSlVCG2ClDD+ClppdVCGCmqqdVCGCmqrdVCG.The aerodynamic angles of attackα and β are then defined by the following:(14)tanα=wu,(15)tanβ=vu2+w2.The above equations are the major equations of the projectile flight dynamic model. More details regarding the projectile flight model can be found in [33]. All aerodynamic coefficients CD0,CD2,Clα,Cypα,Cmα,Cnpα,ClDD,Clp,andCmq depend upon the local Mach number and are obtained during simulation by linear interpolation. To solve this model, the fourth-order Runge–Kutta algorithm is employed in this article. In addition, all calculations are performed with a real record of meteorological data.


## 3. DE Algorithm-Based Aerodynamic Coefficient Identification

As mentioned above, the aerodynamic coefficient identification problem can be posed as an optimization problem; the basic task is to seek the set of coefficientsθ that best fit the measured data, namely, to minimize the cost function Oθ at different Mach numbers. In this article, the maximum likelihood criterion is chosen to formulate the object cost function, and thus, Oθ is defined as follows:(16)Oθ=∑i=1nεTiR−1εi+lnR,where θ is the vector of unknown coefficients and εi is the deviation of the measured state vector from the calculated state vector, which is defined as follows:(17)εi=HM−HC,where HM and HC are the measured and calculated state vectors, respectively, which are obtained from the 6-DOF dynamic simulation with the identified coefficients in one step. R is the covariance matrix of measurement noise. As the statistical property of the measurement noise is unknown, the optimum estimation R^ usually replaces R:(18)R≈R^=1n∑i=1nεiεTi.The detailed implementation steps for the identification of aerodynamic coefficients using the DE algorithm can be described as follows.Step 1 (initialization operation): generateNPi=1⋯NP individuals randomly in the whole search space. The initial population is given by the following:(19)θki,G=θLi,G+r×θLi,G−θUi,G,wherer is a uniformly distributed random number in the interval 0,1 and θki,G is the kth component of the individual θi in the G generation. θLi,G and θUi,G are the lower and upper bounds of θki,G, respectively.Step 2: evaluate the fitness of each individual in the population according to (16).Step 3: create a new population as follows:(a)
Mutation operation: a mutant vector is generated according to the following:(20)Vki,G=θkr1,G+Fθkr2,G−θkr3,G.(b)
Crossover operation: a new trail vector is generated by mixing the mutant vector and the target vector.(21)Uki,G=Vki,G,ifrk≤CRi,Gork=krand,θki,G,otherwise,whererk is a random number taken from the uniform distribution 0,1.(c)
Selection operation: selection determines whether the target or the trial vector survives to the next iteration. For each mutant vector and its corresponding target vector in the current population, the object function values are compared using a greedy selection scheme. The selection operation is described as follows:(22)θki,G=Uki,G+1,ifOUki,G+1<OUki,G,θki,G,otherwise.Step 4: repeat steps 2-3 until the termination criteria are satisfied. In this article, the termination criteria are defined as follows:(23)G>Gmax.The flowchart of the implemented DE algorithm is shown in Figure2.Figure 2
Flowchart of the implemented DE algorithm.


## 4. Identification Schemes

The first scheme, which is the whole discrete point (WDP) scheme, as shown in Figure3, treats the measured trajectory as a whole. The trajectory would not be divided into several intervals, and all of the coefficients corresponding to different Mach numbers in an aerodynamic coefficients table would be identified simultaneously. In this scheme, the design variable is defined as follows:(24)θ=CM1,CM2,…,CMj,…CMJ−1,CMJT,where CMj represents the coefficient corresponding to the jth Mach number in the coefficients table, j=1,2,…,J, and J is the number of Mach numbers. The dimension of an individual in population k is equal to the number of Mach numbers J.Figure 3
Schematic diagram of the identification schemes.In the second scheme, which is the differential interval linear (DIL) scheme, as shown in Figure3, the whole trajectory would be divided into several intervals. The relationship between the coefficients and the Mach numbers in an interval is treated as a linear function, where the slope does not have to be zero. A constant factor KMn in each interval is introduced in this scheme, which is defined as follows:(25)KMn=CrealCprior.As seen from (25), KMn represents the gap between the prior value and the real value of aerodynamic coefficients, and Creal can be determined by identifying KMn in each interval successively from the free flight data. Namely, KMn becomes the design variable in this optimization problem, which is defined as follows:(26)θ=KMn,where n=1,2,…,N, and N is the number of intervals. It is easy to find that the dimension of an individual in population k in this scheme is equal to the number of the correction factors in each interval.In addition, for comparison research, the conventional differential interval constant (DIC) scheme is used in this paper. This scheme divides the whole trajectory into several intervals as in the DIL scheme; yet, the coefficients in each interval are modeled as constants and identified successively. In this scheme, the design variable is defined as follows:(27)θ=CMi,where CMi represents the coefficient in the ith interval, i=1,2,…,I, and I is the number of intervals. The dimension of an individual in population k in this scheme is equal to the number of coefficients in each interval.Additionally, the prior information, namely, the general trend or characteristics, is required for all of the schemes mentioned above. This prior information could be used to divide the Mach number interval, to set lower and upper bounds of the design variable, to obtain the feasibility of identification, etc. In practical engineering, the prior information could be obtained from wind tunnel tests, numerical methods, or engineering experience.


## 5. Results and Discussion

Several comparative tests are conducted in two different scenarios to study the performances of the identification schemes presented in this paper. A fin-stabilized projectile is selected as the test bed, and the physical characteristics of this projectile are listed in Table1. In order to simplify the problem, the present study focuses on the identification of the drag coefficient, and the methods presented in this paper could be easily used for other aerodynamic coefficients. Other aerodynamic coefficients are input as known values.Table 1
The projectile’s physical properties.
m (kg)IA (kg-m2)IC (kg-m2)CG from nose (m)d (m)45.3400.11531.6271.1900.122In the first scenario, the simulated measurement data is generated from a 6-DOF projectile dynamic simulation. To model the random noise, which inevitably exists in the trajectory data measured by radar, the Gaussian noise is added into the simulated position and velocity. Three standard deviations are used equal to 0.2, 1.0, and 5.0 times the nominal standard deviation are considered, and the nominal standard deviationσ is 3 m and 0.01 m/s for the position and velocity, respectively. The identification is conducted 50 times for each scheme, and the final identification results are the average values of the 50 identification results. The meteorological data and the initial launch conditions used in this scenario are chosen to be the same as the second scenario, which will be described in the next paragraph.In the second scenario, phased array radar data from a free flight test is used as the measurement data to identify the aerodynamic drag coefficient. In this test, available measurements are the position data and radial velocity data, and the meteorological data were recorded with both ground and balloon-carried instruments throughout the test events. In Figure4, the recorded data for the pressure, temperature, wind speed, and wind direction as a function of the altitude are shown. For the present projectile, the meteorological data below 12 km is mainly used. In Table 2, the initial launch conditions of the test are listed.Figure 4
Meteorological data versus altitude. (a) Temperature and pressure versus altitude. (b) Wind speed and wind direction versus altitude.
(a)(b)Table 2
Initial conditions for the two scenarios.
StateFirst scenario/second scenarioLaunching elevation angle (degrees)42.7Launching azimuth angle (degrees)358.953Velocity (m/s)41.317Spin rate (rad/s)25.133Table3 describes the initialization of the parameters implemented in the DE algorithm. The main control parameters that need to be set are the number of individuals NP, the weight factor F, and the crossover rate CR. Every control parameter is tuned through appropriate tradeoff studies for the optimal process. Detailed information about the principle of tuning the control parameters can be found in [27, 28, 30, 35]. In addition, to improve the efficiency of identification processing, the multiprocessing technique is used in the program.Table 3
Parameter initialization for the three schemes.
ParametersIdentification schemesDILWDPDICNumber of individualsNP8016080Weight factorF0.8Crossover rateCR0.9Maximum number of generationsGmax300Mutation strategyDE/current-to-best/1/binSelection strategyBest-parent-childTermination criteriaGmaxNumber of measurement data in interval25All25Number of parallel processes15In Figure5, the reconstructions of the trajectory with the identification results from two scenarios are plotted, including the altitude versus range and radial velocity versus range. In the first scenario (Figures 5(a)–5(c)), as can be observed, the match error enlarged with the increase of the noise. In Figure 5(c), relatively significant error can be found, especially in the apex of the trajectory. There is a possibility of high sensitivity of the drag coefficient in the transonic region. In Figure 5(d), in the descent stage of the trajectory, the error appears to be enlarged. It is possible that due to the limit of the performance of the radar system, the noise existing in the measured radar trajectory data enlarged with the increase of the range, which reduced the accuracy of identification at the descent stage.Figure 5
The trajectory matching for different schemes from the two scenarios. (a)σnoise=0.2σ. (b) σnoise=1.0σ. (c) σnoise=5.0σ. (d) Radar data.
(a)(b)(c)(d)In Figure6, the identification results from the first scenario at varying levels of noise are plotted. It is used for closer inspection of the performance of the identification schemes introduced in the present paper. It can be seen that the WDP scheme gives the most accurate result. This is consistent with the results presented in Figure 5. As can be observed from Figure 5, the reconstruction of the trajectory of the WDP schemes matches the true data very well whether in the apex region or the descent stage of the two scenarios. This is because more measurements are considered in the WDP scheme, which enable the sampled noise to approach a normal distribution, and then leads to a more consistent nominal cost. This advantage can avoid the overfitting of data and increase the accuracy of the identification. Additionally, because the WDP scheme identifies all of the coefficients at different Mach numbers simultaneously, the error generated in the early stage of the trajectory would not be accumulated and spread to the later stage. In contrast, since the DIC scheme and the DIL scheme identified the coefficients along the trajectory interval successively, the coupling error would gradually accumulate and reduce the identification accuracy. Additionally, in each interval of the aerodynamic coefficients table, the coefficients are obtained by linear interpolation. It is actually that the coefficients are treated as linear function of Mach number, and the WDP scheme identifies both the slope and intercept. Namely, the WDP scheme can use the first-order information of the model of aerodynamic coefficients.Figure 6
Identification results from simulated radar data for three schemes. (a)σnoise=0.2σ. (b) σnoise=1.0σ. (c) σnoise=5.0σ.
(a)(b)(c)For the DIL scheme, it is intended to improve the identification accuracy. However, as can be observed from Figures5 and 6, the accuracy of the DIL scheme is much closer to that of the conventional DIC scheme. It is found that the sums of the absolute error of the range from the second scenario for the DIL scheme and the DIC scheme are 1.5E + 5 and 2.15E + 5, respectively, whereas it is 3.1E + 4 for the WDP scheme. It is possible that the prior information is not accurate enough. Although the DIL scheme models the coefficients in an interval with a linear function instead of a constant function as in the DIC scheme, the specific form of the linear function, namely, the prior information, may not be close enough to the real value. Moreover, the introduced correction factor would simultaneously affect the slope and intercept, which would lose some accuracy. In the future work, it is suggested to choose both the slope and intercept as design variables to be identified in the WDP scheme. Moreover, in Figure 6(c), significant fluctuation from the results of the DIC scheme can be observed, especially in the transonic region. This is because the DIC scheme uses constants to model the aerodynamic coefficients, which means that the converging order is zero order. In the field of optimization, it also means that the robustness of the DIC scheme is worse than that of the WDP scheme and the DIL scheme.In Figure7, the computation CPU times using different schemes with simulated data and radar data are plotted. In all of the tests, the computation CPU times for the DIC scheme and the DIL scheme are very close to each other, which varied from 75 s to 81 s. For the WDP scheme, the computation time is much longer than that of the other two schemes. As seen from Figure 7, the computation CPU time for the WDP scheme is up to 211 s.Figure 7
Mean computation CPU times.A direct reason for this situation is that the WDP scheme has larger computational budgets in each iteration of the optimization process since the dimension of an individual in the population is much higher than that in the DIC scheme and the DIL scheme. In the present study, we choose the drag coefficient as the identified parameter. According to (24), (26), and (27), the dimension of the identification problem for the DIC scheme and the DIL scheme is one, and the dimension for the WDP scheme depends upon the number of Mach numbers in the aerodynamic coefficients table. In the present study, the number of Mach numbers is 9.Although the WDP scheme has advantages with regard to accuracy and robustness compared with the other two schemes, it is in exchange of greater computational effort. However, with the development of computer science and technology, an acceptable level of the computational time could be achieved by utilizing high-performance computing technology, such as the parallel computation used in the present paper. Additionally, one can improve the efficiency of the identification by using multiple schemes. For example, in the different Mach number regions, one can apply the different schemes depending upon the degree of fluctuation of the aerodynamic coefficients. For the WDP scheme, one can divide the trajectory into two or three intervals to improve the efficiency while maintaining the accuracy in cases, where the number of Mach numbers in the aerodynamic coefficients table is much more than 9, as in the present paper.


## 6. Conclusions

Two identification schemes based on DE algorithm-based aerodynamic coefficient identification are proposed. Several comparative tests are conducted by identifying the aerodynamic coefficients from the noisy simulated flight data and radar data. The results show that both schemes presented in this paper can accurately estimate the aerodynamic coefficients. The WDP scheme and the DIL scheme improve the robustness of the identification based on the DE algorithm. Since the WDP scheme treats the trajectory as a whole, which effectively reduces the coupling error, the WDP scheme gives relatively better identification results compared with the DIL scheme and the DIC scheme. However, the computation times are more than two times greater compared with the other two schemes.For application of the DE algorithm-based projectile aerodynamic coefficient identification in practical engineering, it is suggested that future work should focus on the following objectives. First, the fusion of multiple schemes or a new scheme might be developed to improve the overall identification performance with limited measurements. Second, the efficiency of the DE algorithm should be further improved, especially in dealing with multidimensional problems.



---

*Source: 2545878-2020-08-05.xml*