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.
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. Figure 1 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 from an inertial frame and the standard aerospace sequence Euler angles .
The kinematic and dynamic differential equations for describing the motion of the projectile in flight are given as follows [1, 33]:
The forces acting on the projectile in (3) consist of gravity and the aerodynamic force. The aerodynamic force is split into a standard aerodynamic force and Magnus aerodynamic force . Gliding flight is assumed in this study. The combination of forces is expressed in the following equation:
In (7), the drag coefficient is estimated by considering the zero-yaw drag coefficient and the yaw drag coefficient [2]:
It is known that has the most significant effect on the trajectory, and represents the main part of [34].
The moment acting on the projectile in (4) consists of the moment due to , the moment due to , and the unsteady aerodynamic moments , as shown in the following equation:
The moments in (10) are expressed as follows:
The aerodynamic angles of attack and are then defined by the following:
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 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 at different Mach numbers. In this article, the maximum likelihood criterion is chosen to formulate the object cost function, and thus, is defined as follows:where is the vector of unknown coefficients and is the deviation of the measured state vector from the calculated state vector, which is defined as follows:where and are the measured and calculated state vectors, respectively, which are obtained from the 6-DOF dynamic simulation with the identified coefficients in one step. is the covariance matrix of measurement noise. As the statistical property of the measurement noise is unknown, the optimum estimation usually replaces :
The detailed implementation steps for the identification of aerodynamic coefficients using the DE algorithm can be described as follows. Step 1 (initialization operation): generate individuals randomly in the whole search space. The initial population is given by the following: where is a uniformly distributed random number in the interval and is the component of the individual in the generation. and are the lower and upper bounds of , 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:(b)Crossover operation: a new trail vector is generated by mixing the mutant vector and the target vector. where is a random number taken from the uniform distribution .(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: Step 4: repeat steps 2-3 until the termination criteria are satisfied. In this article, the termination criteria are defined as follows:
The flowchart of the implemented DE algorithm is shown in Figure 2.
4. Identification Schemes
The first scheme, which is the whole discrete point (WDP) scheme, as shown in Figure 3, 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:where represents the coefficient corresponding to the Mach number in the coefficients table, , and is the number of Mach numbers. The dimension of an individual in population is equal to the number of Mach numbers .
In the second scheme, which is the differential interval linear (DIL) scheme, as shown in Figure 3, 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 in each interval is introduced in this scheme, which is defined as follows:
As seen from (25), represents the gap between the prior value and the real value of aerodynamic coefficients, and can be determined by identifying in each interval successively from the free flight data. Namely, becomes the design variable in this optimization problem, which is defined as follows:where , and is the number of intervals. It is easy to find that the dimension of an individual in population 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:where represents the coefficient in the interval, , and is the number of intervals. The dimension of an individual in population 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 Table 1. 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.
In 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 Figure 4, 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.
(a)
(b)
Table 3 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 , the weight factor , and the crossover rate . 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.
In Figure 5, 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.
(a)
(b)
(c)
(d)
In Figure 6, 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.
(a)
(b)
(c)
For the DIL scheme, it is intended to improve the identification accuracy. However, as can be observed from Figures 5 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 Figure 7, 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.
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.
Nomenclature
Projectile Dynamic Model| : | Projectile inertia components |
| : | Projectile angular velocity vector components |
| : | Projectile linear velocity vector components |
| : | Projectile location vector components |
| : | Projectile aerodynamic angles of attack |
| : | Projectile Euler angles |
| : | Aerodynamic force vector components |
| : | Aerodynamic moment vector components |
| : | Atmospheric density |
| : | Projectile reference diameter |
| : | Projectile reference length |
| : | Projectile reference area |
| : | Projectile center of gravity location from the nose |
| : | Projectile mass |
| : | Magnitude of the mass center velocity |
| : | Projectile radial velocity |
| : | Drag coefficient |
| : | Zero-yaw drag coefficient |
| : | Yaw drag coefficient |
| : | Lift force coefficient derivative with respect to the total yaw angle |
| : | Magnus force aerodynamic coefficient |
| : | Pitching moment aerodynamic coefficient |
| : | Magnus moment aerodynamic coefficient |
| : | Roll moment aerodynamic coefficient due to fin cant |
| : | Roll rate damping moment aerodynamic coefficient |
| : | Pitch rate damping moment aerodynamic coefficient . |
| : | State vector of the ballistic model |
| : | Deviations of the measured state from the calculated state |
| : | Object function |
| : | Number of the measured points |
| : | Mach number |
| : | Mach number interval |
| : | , vector to be perturbed; , number of difference vectors; , type of crossover |
| : | Number of individuals |
| : | Number of generations |
| : | Parameters within the vectors |
| : | Mutant vector |
| : | Target vector |
| : | Weight factor |
| : | Crossover rate. |
| : | Earth frame |
| : | Aerodynamic drag force vector components |
| : | Calculated and measured states of the projectile |
| : | Lower and upper limits of the search space. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the Natural Science Foundation of North University of China (Grant no. XJJ201806).