# Wind Identification of Spinning Projectile Using Improved Multiobjective Differential Evolution Algorithm

Authors: Dingye Zhang; Jun Guan; Wenjun Yi; Youran Xia; Kangjian Wang
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

Wind velocity has an important influence on the ballistic characteristics of uncontrolled projectiles. It is difficult to precisely determine the projectile’s impact position if wind velocity information from the projectile’s flight process cannot be collected. A wind velocity identification technique of spinning projectile based on the multiobjective chaotic adaptive differential evolution algorithm is suggested to increase the estimation accuracy of wind velocity and ballistic prediction accuracy during projectile flight. The variation law of projectile aerodynamic characteristics under no wind situation is calculated using the 4D kinematic model of spinning projectile. Three aerodynamic parameter coefficients are chosen as reference variables, three objective functions are defined by mean square error, and the two components of wind velocity along the ground coordinate system are used as decision variables in the identification process. The study identifies the two components, which is based on the multiobjective chaotic adaptive differential evolution algorithm. Several groups of wind velocity identification under constant and variable wind conditions are numerically simulated. The results show that the suggested method can estimate wind velocity effectively and precisely throughout the flight of spinning projectile.


---

## Body

## 1. Introduction

Wind velocity, being a significant meteorological factor, has a significant influence on projectile flight ballistic. It is one of the most direct and significant influences on projectile aerodynamic properties, particularly for conventional spinning projectile and in the uncontrolled flight phase of ballistic correction projectile. Wind velocity has two characteristics: speed and direction. The word “wind direction” refers to the direction in which the wind is blowing. The azimuth of the wind direction is the angle formed by rotating clockwise from due north. The magnitude and direction of the pressure gradient, the Coriolis inertia force, the air viscosity coefficient, and the roughness of the surface are all factors that influence wind direction and speed [1]. As a result, describing the dispersion of wind using a single mathematical formula is challenging.In order to facilitate the study of the wind velocity in the process of projectile flight, the changing wind is generally averaged over a period of time, during which the wind speed vector is constant. The average wind is the name given to this average number. The gust is the difference between the instantaneous wind and this average value. The average wind throughout the ballistic causes each projectile’s ballistic to have the same deviation as a windless ballistic, such as range and direction deviation. However, the gust causes each projectile to have a distinct deviation, which is one of the major causes of ballistic dispersion [1]. To facilitate the analysis of the problem, this study assumes that the wind direction is parallel to the ground and ignores the impact of the vertical component of wind speed. Wind velocity is characterized by longitudinal wind and cross wind. The wind velocity component parallel to the shooting plane is called longitudinal wind and positive downwind and is denoted by Wx. Longitudinal wind causes not only range change but also directional deviation. The wind velocity component perpendicular to the shooting plane is called cross wind, which is positive from left to right when viewed in the shooting direction. It is expressed by Wz. The cross wind not only causes side deflection but also affects the range.Although the Beidou navigation system and inertial navigation system carried by the spinning projectile provide better information about the velocity and position of the projectile relative to the navigation coordinate system, they do not provide information about wind velocity. In actual flight tests, the information on wind velocity is generally measured by sounding the balloon. The parameters such as air pressure, air temperature, humidity, and wind direction and speed can be measured from top to bottom after the balloon is released. Similarly, information about wind velocity can also be obtained in the wind tunnel experiment. Figures1–3 simulate the dynamic process of a spinning projectile in actual flight.Figure 1
The state when the rudder is just opened.Figure 2
The state that the rudder is fully opened but not yet ruddered.Figure 3
Sine ruddering state.Figure1 shows the state when the rudder of the projectile is just opened, Figure 2 shows the state when the rudder is fully opened but not yet ruddered, and Figure 3 shows the state when the rudder is ruddered sinusoidally. At the later stage of sine ruddering, the rudder generates a hinge moment, so that the projectile receives an upward force in the direction of the combined force, at which time the projectile “heads up” to achieve the purpose of later gliding range. During the experiment, information about the wind velocity can be obtained through the wind velocity measurement sensor.In terms of simulation tests, in the 1990s, Costello [2] conducted a detailed study of the ballistic simulation model for duck-rudder guided mortars under windy conditions, which improved the accuracy of the ballistic drop point under conditions that the wind speed was known. In 2008, Zhao et al. [3] investigated the effect of the variation of longitudinal wind and cross wind on the range of the end-guided projectile considering the presence of longitudinal wind and cross wind. However, all of the above methods require the determination of wind velocity before shooting. If wind velocity is not measured, wind velocity needs to be identified as a way to obtain wind velocity information.Aiming at the identification of wind velocity, Zhang et al. [4] proposed a method to identify the longitudinal wind speed online based on the measured value of the onboard sensor and the theoretical value of the projectile drag coefficient, but the method is only applicable to the identification of the longitudinal wind speed and cannot estimate the cross wind speed value. Sun and Yang [5] proposed a method of wind speed identification assisted by overload control. This method uses the function of overload angular rate control loop to control the angle of attack to zero and then carries out online identification of wind speed according to the motion law of missile, which is applicable to guided munitions with strong overload control ability. The problem is that it is more restricted in use and cannot be applied to the uncontrolled projectile. For uncontrolled projectile, Gao et al. [6] proposed a method to identify longitudinal wind and cross wind by using sequential quadratic programming based on ballistic measurement data, but the wind speed estimation error of this method is large when the ballistic height is large, which exceeds ±10% of the actual wind speed.To address the problem of wind velocity prediction, Mok et al. [7] proposed an adaptive Kalman filter-based storm surge prediction model and applied it to the Macau inner harbor for validation. The storm surge was successfully predicted for the next 3 hours, but its predicted wind velocity magnitude had a maximum peak of 18%, which exceeded the actual wind value by ±10%.In order to improve the situation that there is no available air data sensor and the flight stability and impact point accuracy of spinning projectile considering the effects of longitudinal wind and cross wind, this paper proposes a wind velocity identification method for spinning projectile based on the multiobjective chaotic adaptive differential evolution algorithm. The remainder of this paper is structured as follows. In Section2, the kinematic model and aerodynamic model of spinning projectiles are used to solve the drag coefficient, the derivative of the lift coefficient, and the joint partial derivatives of the Magnus force coefficient of spinning projectile under the conditions of no wind, constant wind, and variable wind. Section 3 details how the multiobjective optimization method is applied to wind velocity estimation and the steps of the proposed new method. In Section 4, the multiobjective chaotic adaptive differential evolution algorithm is used to identify the wind speed over a period of time at a given wind speed to obtain longitudinal and cross wind values. The simulation results show that the proposed algorithm can effectively identify the longitudinal and cross wind values under constant and variable wind conditions and still has a high identification accuracy under variable wind conditions, which can satisfy the identification of wind velocity without available air data sensors to a certain extent.


## 2. Aerodynamic Model of Spinning Projectile

### 2.1. 4D Kinematic Model of Spinning Projectile

One of the most important assumptions of the mass-point ballistic equation of the spinning projectile is that the angle of attackδ=0, i.e., the projectile axis is always in line with the velocity direction. However, the actual situation is not ideal. The projectile axis will be subject to various perturbations (initial perturbations, ballistic bending, asymmetric factors, etc.) and will rotate around the center of mass to form an angle of attack. The induced drag, lift, and Magnus force generated by the angle of attack will make the actual ballistic deviate from the ideal ballistic [1]. Therefore, the use of the mass-point ballistic equation to calculate the projectile aerodynamic parameters must have a large error. In order to reduce the error, the angle of attack of the projectile after perturbation should be considered, so it is necessary to use the complete six degrees of freedom (6D) or its simplified equation to solve it, but the six degrees of freedom equation is more complicated and the computational workload is large. In this paper, the drag, lift, and Magnus force generated by the angle of attack are added to the mass-point ballistic equation, while the angle of attack itself is directly calculated by the analytical solution of the angular motion equation. The improved mass-point ballistic equation is called the modified mass-point ballistic equation (4D). The modified mass-point ballistic equation is generally only suitable for calculating the full ballistic with an angle of shoot less than 60°.In studying the 4D model, it is assumed that the moment acting on the projectile is in equilibrium, at which time the projectile axis is the dynamic equilibrium axis and the angle of attack is the dynamic equilibrium angle. The key to modify the mass-point ballistic equation is to determine the dynamic equilibrium angle vectorαe. Projecting it onto the ground coordinate system O−xyz, the calculation formula is as follows.
(1)αe=αex,αey,αez,(2)αex=x1abgx2VrxVry−x1Vrzx13+x22x1Vr2,αey=−x1x2abgVrx2+Vrz2x13+x22x1Vr2,αez=x1abgx1Vrx+x2VryVrzx13+x22x1Vr2,Among them,(3)x1=1−ρSd2maa−abCz″γ̇Vr2,x2=ρS2maa−abCy′Vr2,(4)aa=2mldγ̇my″ρSlmz′Cy′Vr4+ρSld2Cz″my″γ̇2Vr2,ab=2Cγ̇Cy′ρSlmz′Cy′Vr4+ρSld2Cz″my″γ̇2Vr2,(5)Vr=Vx−Wx,Vy,Vz−Wz.Furthermore, the modified mass-point ballistic equation established in the ground coordinate systemO−xyz is as follows.
(6)dVxdt=−ρS2mVrVx−WxCx0+Cx2⋅αe2+ρS2mVr2Cy′αex+ρSd2mγ̇Cz″αeyVz−Wz−αezVy,dVydt=−ρS2mVrVyCx0+Cx2⋅αe2+ρS2mVr2Cy′αey+ρSd2mγ̇Cz″αezVx−Wx−αexVz−Wz−g,dVzdt=−ρS2mVrVz−WzCx0+Cx2⋅αe2+ρS2mVr2Cy′αez+ρSd2mγ̇Cz″αexVy−αeyVx−Wx,dxdt=Vx,dydt=Vy,dzdt=Vz,dγ̇dt=−ρSld2CVrγ̇mxz′,dγdt=γ̇.In Equations (1)-(6), g is the strength of the gravitational field near the earth’s surface, m is the total mass of the projectile, S is the characteristic area of the projectile, which is usually selected as the maximum cross-sectional area of the projectile, d is the maximum cross-sectional diameter of the projectile, l is the characteristic length of the projectile, and ρ is the atmospheric density. γ is the roll angle and γ̇ is the roll angle speed. Vr is the velocity of the projectile relative to the air, Vx, Vy, and Vz are the three components of the projectile velocity relative to the ground coordinate system O−xyz, Wx and Wz are the Ox axis and Oz axis components of the wind speed relative to the ground coordinate system O−xyz, x, y, and z are the three position coordinates of the projectile relative to the ground coordinate system O−xyz, and αex, αey, αez is the dynamic balance angle vector αe’s three components relative to the ground coordinate system O−xyz; C is the axial moment of inertia of the projectile. Cx0 and Cx2 are the zero lift drag coefficient and induced drag coefficient, respectively, Cx=Cx0+Cx2⋅αe2 is the drag coefficient of the projectile, Cy′ is the derivative of the lift coefficient of the projectile, and Cz″ is the joint partial derivative of the Magnus force coefficient of the projectile; mxz′ is the derivative of extreme damping moment coefficient, mz′ is the derivative of overturning moment coefficient, and my″ is the joint partial derivative of Magnus moment coefficient.


### 2.2. 4D Aerodynamic Model of Spinning Projectile

Based on the 4D kinematic model of spinning projectile, the calculation formulas of the drag coefficient, the derivative of the lift coefficient, and the joint partial derivatives of the Magnus force coefficient are derived.(7)CxCy′Cz″=A1A2A3B1B2B3C1C2C3−12mV̇x2mV̇y+g2mV̇z.Among them,(8)A1=−ρSVrVx−Wx,B1=−ρSVrVy,C1=−ρSVrVz−Wz,(9)A2=ρSVr2αex,B2=ρSVr2αey,C2=ρSVr2αez,(10)A3=ρSdγ̇αeyVz−Wz−αezVy,B3=ρSdγ̇αezVx−Wx−αexVz−Wz,C3=ρSdγ̇αexVy−αeyVx−Wx.Therefore, the aerodynamic coefficient parameter identification and the wind velocity model identification are based on Equation (7) and with the help of the dynamic model. In Section 2.3, the variation law of these three aerodynamic parameters under no wind, constant wind, and variable wind conditions is solved.


### 2.3. Calculation of Aerodynamic Parameters of Spinning Projectile

The drag coefficient of the projectile with time in the no wind (i.e.,Wx=0 and Wz=0), constant wind (e.g., when Wx=5 and Wz=6), and variable wind (when Wx=5+sin2π/10t and Wz=6+cos2π/10t) conditions is shown in Figure 4.Figure 4
Variation curve of the drag coefficient under three conditions.The derivative of the lift coefficient of the projectile as a function of time in the no wind (i.e.,Wx=0 and Wz=0), constant wind (e.g., when Wx=5 and Wz=6), and variable wind (when Wx=5+sin2π/10t and Wz=6+cos2π/10t) conditions is shown in Figure 5.Figure 5
Variation curve of the derivative of the lift coefficient under three conditions.The joint partial derivatives of the Magnus force coefficient of the projectile with time in the no wind (i.e.,Wx=0 and Wz=0), constant wind (e.g., when Wx=5 and Wz=6), and variable wind (when Wx=5+sin2π/10t and Wz=6+cos2π/10t) conditions is shown in Figure 6.Figure 6
Variation curve of the joint partial derivatives of the Magnus force coefficient under three conditions.As can be seen from Figures4–6, compared with the windless state, the derivative of the lift coefficient and the joint partial derivatives of the Magnus force coefficient of the projectile change significantly after the addition of constant wind, while the drag coefficient remains basically unchanged. Compared with the constant wind state, the overall change trend of these three aerodynamic parameters is basically the same after adding trigonometric pulsation perturbation, but fluctuations occur at certain segmental positions, and the fluctuations have a more obvious effect on the derivative of lift coefficient.


## 2.1. 4D Kinematic Model of Spinning Projectile

One of the most important assumptions of the mass-point ballistic equation of the spinning projectile is that the angle of attackδ=0, i.e., the projectile axis is always in line with the velocity direction. However, the actual situation is not ideal. The projectile axis will be subject to various perturbations (initial perturbations, ballistic bending, asymmetric factors, etc.) and will rotate around the center of mass to form an angle of attack. The induced drag, lift, and Magnus force generated by the angle of attack will make the actual ballistic deviate from the ideal ballistic [1]. Therefore, the use of the mass-point ballistic equation to calculate the projectile aerodynamic parameters must have a large error. In order to reduce the error, the angle of attack of the projectile after perturbation should be considered, so it is necessary to use the complete six degrees of freedom (6D) or its simplified equation to solve it, but the six degrees of freedom equation is more complicated and the computational workload is large. In this paper, the drag, lift, and Magnus force generated by the angle of attack are added to the mass-point ballistic equation, while the angle of attack itself is directly calculated by the analytical solution of the angular motion equation. The improved mass-point ballistic equation is called the modified mass-point ballistic equation (4D). The modified mass-point ballistic equation is generally only suitable for calculating the full ballistic with an angle of shoot less than 60°.In studying the 4D model, it is assumed that the moment acting on the projectile is in equilibrium, at which time the projectile axis is the dynamic equilibrium axis and the angle of attack is the dynamic equilibrium angle. The key to modify the mass-point ballistic equation is to determine the dynamic equilibrium angle vectorαe. Projecting it onto the ground coordinate system O−xyz, the calculation formula is as follows.
(1)αe=αex,αey,αez,(2)αex=x1abgx2VrxVry−x1Vrzx13+x22x1Vr2,αey=−x1x2abgVrx2+Vrz2x13+x22x1Vr2,αez=x1abgx1Vrx+x2VryVrzx13+x22x1Vr2,Among them,(3)x1=1−ρSd2maa−abCz″γ̇Vr2,x2=ρS2maa−abCy′Vr2,(4)aa=2mldγ̇my″ρSlmz′Cy′Vr4+ρSld2Cz″my″γ̇2Vr2,ab=2Cγ̇Cy′ρSlmz′Cy′Vr4+ρSld2Cz″my″γ̇2Vr2,(5)Vr=Vx−Wx,Vy,Vz−Wz.Furthermore, the modified mass-point ballistic equation established in the ground coordinate systemO−xyz is as follows.
(6)dVxdt=−ρS2mVrVx−WxCx0+Cx2⋅αe2+ρS2mVr2Cy′αex+ρSd2mγ̇Cz″αeyVz−Wz−αezVy,dVydt=−ρS2mVrVyCx0+Cx2⋅αe2+ρS2mVr2Cy′αey+ρSd2mγ̇Cz″αezVx−Wx−αexVz−Wz−g,dVzdt=−ρS2mVrVz−WzCx0+Cx2⋅αe2+ρS2mVr2Cy′αez+ρSd2mγ̇Cz″αexVy−αeyVx−Wx,dxdt=Vx,dydt=Vy,dzdt=Vz,dγ̇dt=−ρSld2CVrγ̇mxz′,dγdt=γ̇.In Equations (1)-(6), g is the strength of the gravitational field near the earth’s surface, m is the total mass of the projectile, S is the characteristic area of the projectile, which is usually selected as the maximum cross-sectional area of the projectile, d is the maximum cross-sectional diameter of the projectile, l is the characteristic length of the projectile, and ρ is the atmospheric density. γ is the roll angle and γ̇ is the roll angle speed. Vr is the velocity of the projectile relative to the air, Vx, Vy, and Vz are the three components of the projectile velocity relative to the ground coordinate system O−xyz, Wx and Wz are the Ox axis and Oz axis components of the wind speed relative to the ground coordinate system O−xyz, x, y, and z are the three position coordinates of the projectile relative to the ground coordinate system O−xyz, and αex, αey, αez is the dynamic balance angle vector αe’s three components relative to the ground coordinate system O−xyz; C is the axial moment of inertia of the projectile. Cx0 and Cx2 are the zero lift drag coefficient and induced drag coefficient, respectively, Cx=Cx0+Cx2⋅αe2 is the drag coefficient of the projectile, Cy′ is the derivative of the lift coefficient of the projectile, and Cz″ is the joint partial derivative of the Magnus force coefficient of the projectile; mxz′ is the derivative of extreme damping moment coefficient, mz′ is the derivative of overturning moment coefficient, and my″ is the joint partial derivative of Magnus moment coefficient.


## 2.2. 4D Aerodynamic Model of Spinning Projectile

Based on the 4D kinematic model of spinning projectile, the calculation formulas of the drag coefficient, the derivative of the lift coefficient, and the joint partial derivatives of the Magnus force coefficient are derived.(7)CxCy′Cz″=A1A2A3B1B2B3C1C2C3−12mV̇x2mV̇y+g2mV̇z.Among them,(8)A1=−ρSVrVx−Wx,B1=−ρSVrVy,C1=−ρSVrVz−Wz,(9)A2=ρSVr2αex,B2=ρSVr2αey,C2=ρSVr2αez,(10)A3=ρSdγ̇αeyVz−Wz−αezVy,B3=ρSdγ̇αezVx−Wx−αexVz−Wz,C3=ρSdγ̇αexVy−αeyVx−Wx.Therefore, the aerodynamic coefficient parameter identification and the wind velocity model identification are based on Equation (7) and with the help of the dynamic model. In Section 2.3, the variation law of these three aerodynamic parameters under no wind, constant wind, and variable wind conditions is solved.


## 2.3. Calculation of Aerodynamic Parameters of Spinning Projectile

The drag coefficient of the projectile with time in the no wind (i.e.,Wx=0 and Wz=0), constant wind (e.g., when Wx=5 and Wz=6), and variable wind (when Wx=5+sin2π/10t and Wz=6+cos2π/10t) conditions is shown in Figure 4.Figure 4
Variation curve of the drag coefficient under three conditions.The derivative of the lift coefficient of the projectile as a function of time in the no wind (i.e.,Wx=0 and Wz=0), constant wind (e.g., when Wx=5 and Wz=6), and variable wind (when Wx=5+sin2π/10t and Wz=6+cos2π/10t) conditions is shown in Figure 5.Figure 5
Variation curve of the derivative of the lift coefficient under three conditions.The joint partial derivatives of the Magnus force coefficient of the projectile with time in the no wind (i.e.,Wx=0 and Wz=0), constant wind (e.g., when Wx=5 and Wz=6), and variable wind (when Wx=5+sin2π/10t and Wz=6+cos2π/10t) conditions is shown in Figure 6.Figure 6
Variation curve of the joint partial derivatives of the Magnus force coefficient under three conditions.As can be seen from Figures4–6, compared with the windless state, the derivative of the lift coefficient and the joint partial derivatives of the Magnus force coefficient of the projectile change significantly after the addition of constant wind, while the drag coefficient remains basically unchanged. Compared with the constant wind state, the overall change trend of these three aerodynamic parameters is basically the same after adding trigonometric pulsation perturbation, but fluctuations occur at certain segmental positions, and the fluctuations have a more obvious effect on the derivative of lift coefficient.


## 3. Wind Estimation Method

### 3.1. Multiobjective Optimization Design

It can be seen from Section2.3 that when the values of Wx and Wz are changed, it causes a change in the three aerodynamic parameters. This means that three models can be proposed to search the solutions of two wind components to minimize the error of a particular aerodynamic model. Therefore, the problem is a specific three-objective optimization problem for which a unique solution should be obtained if the actual wind can be found. In the actual simulation experiments, a solution set close to the real wind will be obtained after optimization.It is a very common problem in engineering problems to deal with situations where multiple objectives (including physical constraints, operational constraints, and nonlinear constraints) need to be optimized. Therefore, solving these problems from the classical single-objective optimization perspective is difficult to achieve. Thus, multiobjective optimization methods are introduced. The multiobjective optimization problem (MOP) can be expressed as(11)minμ∈Rfμ=f1μ,f2μ,⋯,fmμ,where μ is the solution that simultaneously minimizes m cost functions fii=1,2,⋯,m. Generally, it is impossible to find solutions that satisfy all requirements at once, so the optimizer must simultaneously provide the set of solutions of all objectives that are not improved by any other set. This set of solution sets is called Pareto set Θ∗, and their values in the objective space are called Pareto front fp∗.Definition 1 (Pareto optimality [8]).
An objective vectorfμ1 is Pareto optimal if there exist no other objective vectors fμ2 such that for all i∈1,2,⋯,m there is fiμ2≤fiμ1, and for at least one j∈1,2,⋯,m there is fjμ2<fjμ1.Definition 2 (strict dominance [8]).
An objective vectorfμ1 is dominated by another objective vector fμ2 if fiμ2<fiμ1 for all i∈1,2,⋯,m.Definition 3 (dominance [8]).
An objective vectorfμ1 dominates another vector fμ2 if fμ1 is no worse than fμ2 in all objectives and is better in at least one objective, i.e., fμ1≺fμ2.Figure7 illustrates how to adopt multiobjective optimization as the design method [9]. The process consists of three stages: MOP definition, algorithm optimization, and analysis and decision-making [10]. The design must be considered as a holistic effort in which each stage has equal importance for the successful execution of the design process [11]. Therefore, the objectives, decision space, and constraints must be selected carefully during the MOP definition phase in order to optimize the correct problems during the optimization process. Finally, once the approximation of Pareto front is obtained, an in-depth analysis should be performed in order to find the most appropriate solution in the decision-making phase. The same topology is followed in the wind estimation process.Figure 7
Multiobjective optimization design flow chart.


### 3.2. Definition of Multiobjective Optimization Problem

The 4D aerodynamic model of spinning projectile in Section2.2 is used as the identification model from which the values of the two wind velocity components are identified. If the solutions that minimize three model errors simultaneously are found in the two components of the wind velocity, the solutions are likely to be the actual wind experienced by the projectile during flight. Three objectives are defined, and each aerodynamic model uses the mean square error (MSE) as the performance index of the identification process; then, three cost functions can be defined for each simulation experiment.
(12)f1Wx,z=1N∑h=1NCxth,Wx,Wz−C∧xth,Wx,Wz2,(13)f2Wx,z=1N∑h=1NCy′th,Wx,Wz−C∧y′th,Wx,Wz2,(14)f3Wx,z=1N∑h=1NCz″th,Wx,Wz−C∧z″th,Wx,Wz2.In Equations (12)–(14), C^x, C^y′, and C^z″ are the calculated values provided by the aerodynamic model for Cx, Cy′, and Cy″, respectively. These three cost functions constitute the objective space, and the two wind velocity components Wx and Wz define the decision space. Finally, constraints can be included in objectives and decision variables. In this simulation experiment, constraints are introduced only in the decision space in order to reduce the space of possible solutions.


### 3.3. Multiobjective Optimization Process

In the last two decades, mainstream multiobjective optimization methods usually use evolutionary algorithms, because they have the advantage of being less susceptible to the essential characteristics of the decision space and objective space of multiobjective optimization problems and the ability to obtain multiple Pareto optimal solutions by running the algorithm once solving multiobjective optimization problems. According to the evolutionary mechanism, multiobjective evolutionary algorithm (MOEA) can be divided into three categories [12]: decomposition-based MOEA, dominance relationship-based MOEA, and indicator-based MOEA. Among them, the Pareto dominance relationship-based MOEA methods are more common.Decomposition-based MOEA is the conversion of a multiobjective optimization problem into a single-objective optimization problem by combining or aggregating all subobjectives being optimized into a single objective. The indicator-based MOEA guides the search process and the solution selection process based on performance evaluation indicators. And the basic idea of Pareto dominance relationship-based MOEA is to identify all nondominated individuals from the current evolutionary population using a Pareto-based fitness assignment strategy. This approach was first proposed by Goldberg [13] in 1989. Later, multiobjective genetic algorithm (MOGA [14]), nondominated sorting genetic algorithm (NSGA [15]), and niche Pareto genetic algorithm (NPGA [16]) were proposed, which gradually formed the representative first generation of multiobjective evolutionary optimization algorithms. Since the end of the 20th century, the research trend of multiobjective evolutionary optimization algorithms has changed dramatically, strength Pareto evolutionary algorithm (SPEA [17]), Pareto archived evolutionary strategy (PAES [18]), the Pareto envelope-based selection algorithm (PESA [19], and PESA-II [20]). After that, improved algorithms such as SPEA2 [21], NPGA2 [22], NSGA-II [23], orthogonal multiobjective evolutionary algorithm (OMOEA [24]), and entropy-based multiobjective evolutionary algorithm (EMOEA [25]) were proposed one after another, forming the second generation of multiobjective evolutionary optimization algorithms. The algorithms in this period are mainly characterized by elite retention strategies, and some better methods for maintaining population diversity were proposed, such as crowding distance-based methods, spatial supergrid-based methods, and clustering-based methods.Many evolutionary algorithms based on biological incentive mechanisms have been rapidly developed. These algorithms are classified as natural computing, such as particle swarm optimization algorithm, immune evolutionary algorithm, ant colony optimization algorithm, fish swarm algorithm, and artificial bee colony algorithm. The common feature is that they are optimization algorithms designed by drawing on some mechanism, mechanism or rule, etc. from natural phenomena. No evolutionary algorithm is better than other evolutionary algorithms because they have their own advantages and disadvantages. Among many evolutionary algorithms, differential evolutionary algorithm has the advantages of fewer control parameters and faster convergence. It has been successfully applied to single-objective and multiobjective optimization problems in engineering.In this paper, a multiobjective optimization algorithm based on chaotic adaptive differential evolution (multiobjective chaotic adaptive differential evolution (MOCADE)) is proposed to optimize multiple objectives. The algorithm mainly consists of four steps: initializing the population with chaos, performing adaptive differential evolution operation on the population, calculating the Pareto values of population individuals, and filtering population individuals based on Pareto dominance criterion. The flow chart of the algorithm is shown in Figure8.Figure 8
Wind velocity estimation process.


### 3.4. Chaotic Adaptive Differential Evolution Process

#### 3.4.1. Chaos Initializing

Differential evolution algorithm usesXpop real-valued parameters vectors with Nvar dimension to take them as the population of each generation, and each individual is denoted as follows.
(15)xxpopGxpop=1,2,⋯,Xpop,where xpop denotes the sequence of individuals in the population. G denotes evolutionary algebra. Xpop denotes the population size. In this paper, Xpop=200 is taken, and Xpop remains unchanged in the process of minimization. Nvar denotes the number of decision variables, which is Nvar=2 in this paper.In DE algorithm, a common method of population initializing is to randomly select from the values within the given boundary constraints. The convergence speed and efficiency of the algorithm’s search process are affected by the degree of scattering of the initial positions and the uniformity of their positions in the search space. When the solution space is large and the number of individuals in the population is not large enough, it may cause uneven distribution of individuals in initializing population. And chaos is a kind of irregular motion, which comes from nonlinear dynamical system, and its motion state has good ergodicity and randomness. In this paper, the population is initialized by a one-dimensional logistic chaotic map, which has extremely complex dynamical behavior. Its mathematical expression is as follows.(16)xi+1=μ⋅xi⋅1−xi.Among them,μ∈0,4 is called logistic parameter. When x∈0,1, the logistic map is in a chaotic state.Figure9 depicts the possible values of x obtained after iteration for different values of μ for a certain value of the initial condition x0.Figure 9
Value range ofx.As can be seen in Figure9, it can be seen that the closer μ is to 4, the closer the value range of x is to the average distribution over the whole 0,1 interval; therefore, the logistic control parameter μ=4 is selected in this paper.The steps of initializing the population of logistic chaotic map are shown in Algorithm1. Firstly, an Nvar-dimensional vector x0 is randomly generated, and then Xpop vectors xxpop are generated by Equation (16). Finally, each component of xxpop is carried to the value interval xnvarL,xnvarU of the initialization boundary. The mathematical formulation of the carrier operation to generate the parent population is as follows.
(17)xParent0=xnvarL+xxpop⋅xnvarU−xnvarL.Algorithm 1: Logistic chaotic map initialization population.
Input: Nvar, Xpop, μ, xnvarL, xnvarUOutput: xParent01 randomly generate aNvar-dimensional x02 use Equation (16) generate xxpop3 use Equation (17) generate xParent0


#### 3.4.2. Adaptive Mutating

After the initialization operation, the algorithm enters the mutation operation phase. For each target vectorxxpopGxpop=1,2,⋯,Xpop, the mutation strategy of the basic differential evolution algorithm can be described as “DE/rand/1,” and its mutation vector is generated by the following formula.
(18)vxpopG+1=xr1G+F⋅xr2G−xr3G.Among them,xr1G, xr2G, and xr3G are randomly selected mutually dissimilar individuals in the G-th generation population. r1, r2, and r3 are mutually dissimilar integers randomly selected in the interval 1,Xpop, and all are different from xpop. F is a real constant factor that controls the scaling of the deviation variable, usually F∈0,2.DE/rand/1 has many advantages such as simple structure, convenient, and fast using and has very powerful global search capability. The basic differential evolution algorithm takes real constants for the mutation operator in the search process, and the mutation operator is more difficult to determine in the simulation experiment. If the mutation rate is too large, the algorithm is inefficient in searching and the accuracy of the global optimal solution is low. If the variation rate is too small, the population diversity is reduced and the phenomenon of “early maturity” is easy to occur [26]. Therefore, a differential evolution algorithm with adaptive mutation operator is used in this paper.
(19)λ=e1−Gm/Gm+1−G,F=F0×2λ.Among them,F0 denotes the initial mutation operator, and F0=0.4 is taken in this paper. Gm denotes the maximum evolutionary algebra. In this paper, Gm=600. G denotes the current evolutionary algebra. When the algorithm starts searching, the mutation operator is 2F0, which has a large value, and the population maintains good individual diversity. As the algorithm searches, the mutation operator decreases gradually, and the mutation operator approaches F0 at the later stage, so that the population retains good information and thus avoids the optimal solution being destroyed.


#### 3.4.3. Hybrid Crossover Operation

For each individualxxpopG of the population, the basic crossover operation applies binomial crossover to create the test vector uxpopG.
(20)uxpop,nvarG+1=vxpop,nvarG+1,ifRnvar≤CRornvar=nvarrand,xxpop,nvarG+1,ifRnvar>CRornvar≠nvarrand,xpop=1,2,⋯,Xpop;nvar=1,2,⋯,Nvar.Among them,Rnvar denotes a random number uniformly distributed between 0,1. nvarrand denotes a random integer in the range 1,Nvar, which guarantees that at least one gene is inherited from the parent individual. CR is the crossover probability, usually in the interval 0,1.The crossover operation of Equation (20) is expressed as “bin,” and its combination with the mutation operation of Equation (18) is expressed as “DE/rand/1/bin.” This algorithm is very popular in global optimization [27]. A random range crossover operator is used in this paper.
(21)CR=0.5×1+rand0,1.This way the mean value of the crossover operator is kept at 0.75, taking into account the possible random variation in the amplification of the difference vector, which is conducive to maintaining the population diversity.Simulated binary crossover (SBX) was designed by Deb and Agrawal [28], which is used to simulate the search pattern of single point binary crossover operator in real coded representation. It is an important recombination operator, which is applied to various multiobjective optimization algorithms. The mathematical expression is as follows.
(22)ur1G+1=0.5×1+β⋅xr1G+1−β⋅xr2G,ur2G+1=0.5×1−β⋅xr1G+1+β⋅xr2G,wherein ur1G+1 and ur2G+1 are child nodes generated by two parent nodes xr1G and xr2G. β follows the polynomial probability distribution and is calculated as follows.
(23)β=2×rand1/1+η,ifrand≤0.5,12−2×rand1/1+η,otherwise.Among them, the distribution indexη is a predefined nonnegative real number, and in general, the value of η in SBX is usually taken as 15 or 20. In this paper, η takes 20. At this time, SBX has a strong local search capability around the parent [29]. The steps of SBX operation are shown in Algorithm 2.Algorithm 2: SBX operator for the i-th gene bit.
Input: xr1G,xr2G,ηOutput: ur1G+1,ur2G+11 use Equation (23) to generate β2 use Equation (22) to generate ur1G+1, ur2G+1Considering that the binomial crossover strategy has a good global search ability and SBX has a strong parent-centered local search ability, this paper mixes binomial crossover with simulated binary crossover and sets a value ofχ that keeps changing with the number of population iterations to balance the global search ability and local search ability of the algorithm. Let χ=1−0.6×G/Gm. After the algorithm starts to search, before the iteration proceeds to Gm/2, because the value of χ is larger, the crossover strategy of Equation (20) is chosen at this time, and the algorithm has a strong global search ability. Then, the algorithm iteration proceeds to the second half, because the value of χ decreases, the crossover strategy of Equation (22) is chosen, and the algorithm has a strong local search ability. The hybrid crossover strategy that combines the binomial crossover strategy and the simulated binary crossover strategy is expressed as “bin-SBX.” The pseudocode of “bin-SBX” hybrid crossover strategy is shown in Algorithm 3.Algorithm 3: Bin-SBX.
Input: xr1G,xr2G,η,vxpop,nvarG+1,xxpop,nvarG+1Output: xChildG+11 use Equation (23) to generate β2 use Equation (22) to generate u′xpop,nvarG+13 ifχ≥0.74  use Equation (20) to generate xChildG+1=uxpop,nvarG+15  else ifχ<0.76xChildG+1=u′xpop,nvarG+17  end8 end


#### 3.4.4. Selection Operation

Using the basic selection strategy of differential evolution, the test vectors are compared with the target vectors in the current population according to the greedy criterion, and the vectors with smaller objective function values are retained. All individuals in the next generation are better or at least as good as their corresponding individuals in the current population.


## 3.1. Multiobjective Optimization Design

It can be seen from Section2.3 that when the values of Wx and Wz are changed, it causes a change in the three aerodynamic parameters. This means that three models can be proposed to search the solutions of two wind components to minimize the error of a particular aerodynamic model. Therefore, the problem is a specific three-objective optimization problem for which a unique solution should be obtained if the actual wind can be found. In the actual simulation experiments, a solution set close to the real wind will be obtained after optimization.It is a very common problem in engineering problems to deal with situations where multiple objectives (including physical constraints, operational constraints, and nonlinear constraints) need to be optimized. Therefore, solving these problems from the classical single-objective optimization perspective is difficult to achieve. Thus, multiobjective optimization methods are introduced. The multiobjective optimization problem (MOP) can be expressed as(11)minμ∈Rfμ=f1μ,f2μ,⋯,fmμ,where μ is the solution that simultaneously minimizes m cost functions fii=1,2,⋯,m. Generally, it is impossible to find solutions that satisfy all requirements at once, so the optimizer must simultaneously provide the set of solutions of all objectives that are not improved by any other set. This set of solution sets is called Pareto set Θ∗, and their values in the objective space are called Pareto front fp∗.Definition 1 (Pareto optimality [8]).
An objective vectorfμ1 is Pareto optimal if there exist no other objective vectors fμ2 such that for all i∈1,2,⋯,m there is fiμ2≤fiμ1, and for at least one j∈1,2,⋯,m there is fjμ2<fjμ1.Definition 2 (strict dominance [8]).
An objective vectorfμ1 is dominated by another objective vector fμ2 if fiμ2<fiμ1 for all i∈1,2,⋯,m.Definition 3 (dominance [8]).
An objective vectorfμ1 dominates another vector fμ2 if fμ1 is no worse than fμ2 in all objectives and is better in at least one objective, i.e., fμ1≺fμ2.Figure7 illustrates how to adopt multiobjective optimization as the design method [9]. The process consists of three stages: MOP definition, algorithm optimization, and analysis and decision-making [10]. The design must be considered as a holistic effort in which each stage has equal importance for the successful execution of the design process [11]. Therefore, the objectives, decision space, and constraints must be selected carefully during the MOP definition phase in order to optimize the correct problems during the optimization process. Finally, once the approximation of Pareto front is obtained, an in-depth analysis should be performed in order to find the most appropriate solution in the decision-making phase. The same topology is followed in the wind estimation process.Figure 7
Multiobjective optimization design flow chart.


## 3.2. Definition of Multiobjective Optimization Problem

The 4D aerodynamic model of spinning projectile in Section2.2 is used as the identification model from which the values of the two wind velocity components are identified. If the solutions that minimize three model errors simultaneously are found in the two components of the wind velocity, the solutions are likely to be the actual wind experienced by the projectile during flight. Three objectives are defined, and each aerodynamic model uses the mean square error (MSE) as the performance index of the identification process; then, three cost functions can be defined for each simulation experiment.
(12)f1Wx,z=1N∑h=1NCxth,Wx,Wz−C∧xth,Wx,Wz2,(13)f2Wx,z=1N∑h=1NCy′th,Wx,Wz−C∧y′th,Wx,Wz2,(14)f3Wx,z=1N∑h=1NCz″th,Wx,Wz−C∧z″th,Wx,Wz2.In Equations (12)–(14), C^x, C^y′, and C^z″ are the calculated values provided by the aerodynamic model for Cx, Cy′, and Cy″, respectively. These three cost functions constitute the objective space, and the two wind velocity components Wx and Wz define the decision space. Finally, constraints can be included in objectives and decision variables. In this simulation experiment, constraints are introduced only in the decision space in order to reduce the space of possible solutions.


## 3.3. Multiobjective Optimization Process

In the last two decades, mainstream multiobjective optimization methods usually use evolutionary algorithms, because they have the advantage of being less susceptible to the essential characteristics of the decision space and objective space of multiobjective optimization problems and the ability to obtain multiple Pareto optimal solutions by running the algorithm once solving multiobjective optimization problems. According to the evolutionary mechanism, multiobjective evolutionary algorithm (MOEA) can be divided into three categories [12]: decomposition-based MOEA, dominance relationship-based MOEA, and indicator-based MOEA. Among them, the Pareto dominance relationship-based MOEA methods are more common.Decomposition-based MOEA is the conversion of a multiobjective optimization problem into a single-objective optimization problem by combining or aggregating all subobjectives being optimized into a single objective. The indicator-based MOEA guides the search process and the solution selection process based on performance evaluation indicators. And the basic idea of Pareto dominance relationship-based MOEA is to identify all nondominated individuals from the current evolutionary population using a Pareto-based fitness assignment strategy. This approach was first proposed by Goldberg [13] in 1989. Later, multiobjective genetic algorithm (MOGA [14]), nondominated sorting genetic algorithm (NSGA [15]), and niche Pareto genetic algorithm (NPGA [16]) were proposed, which gradually formed the representative first generation of multiobjective evolutionary optimization algorithms. Since the end of the 20th century, the research trend of multiobjective evolutionary optimization algorithms has changed dramatically, strength Pareto evolutionary algorithm (SPEA [17]), Pareto archived evolutionary strategy (PAES [18]), the Pareto envelope-based selection algorithm (PESA [19], and PESA-II [20]). After that, improved algorithms such as SPEA2 [21], NPGA2 [22], NSGA-II [23], orthogonal multiobjective evolutionary algorithm (OMOEA [24]), and entropy-based multiobjective evolutionary algorithm (EMOEA [25]) were proposed one after another, forming the second generation of multiobjective evolutionary optimization algorithms. The algorithms in this period are mainly characterized by elite retention strategies, and some better methods for maintaining population diversity were proposed, such as crowding distance-based methods, spatial supergrid-based methods, and clustering-based methods.Many evolutionary algorithms based on biological incentive mechanisms have been rapidly developed. These algorithms are classified as natural computing, such as particle swarm optimization algorithm, immune evolutionary algorithm, ant colony optimization algorithm, fish swarm algorithm, and artificial bee colony algorithm. The common feature is that they are optimization algorithms designed by drawing on some mechanism, mechanism or rule, etc. from natural phenomena. No evolutionary algorithm is better than other evolutionary algorithms because they have their own advantages and disadvantages. Among many evolutionary algorithms, differential evolutionary algorithm has the advantages of fewer control parameters and faster convergence. It has been successfully applied to single-objective and multiobjective optimization problems in engineering.In this paper, a multiobjective optimization algorithm based on chaotic adaptive differential evolution (multiobjective chaotic adaptive differential evolution (MOCADE)) is proposed to optimize multiple objectives. The algorithm mainly consists of four steps: initializing the population with chaos, performing adaptive differential evolution operation on the population, calculating the Pareto values of population individuals, and filtering population individuals based on Pareto dominance criterion. The flow chart of the algorithm is shown in Figure8.Figure 8
Wind velocity estimation process.


## 3.4. Chaotic Adaptive Differential Evolution Process

### 3.4.1. Chaos Initializing

Differential evolution algorithm usesXpop real-valued parameters vectors with Nvar dimension to take them as the population of each generation, and each individual is denoted as follows.
(15)xxpopGxpop=1,2,⋯,Xpop,where xpop denotes the sequence of individuals in the population. G denotes evolutionary algebra. Xpop denotes the population size. In this paper, Xpop=200 is taken, and Xpop remains unchanged in the process of minimization. Nvar denotes the number of decision variables, which is Nvar=2 in this paper.In DE algorithm, a common method of population initializing is to randomly select from the values within the given boundary constraints. The convergence speed and efficiency of the algorithm’s search process are affected by the degree of scattering of the initial positions and the uniformity of their positions in the search space. When the solution space is large and the number of individuals in the population is not large enough, it may cause uneven distribution of individuals in initializing population. And chaos is a kind of irregular motion, which comes from nonlinear dynamical system, and its motion state has good ergodicity and randomness. In this paper, the population is initialized by a one-dimensional logistic chaotic map, which has extremely complex dynamical behavior. Its mathematical expression is as follows.(16)xi+1=μ⋅xi⋅1−xi.Among them,μ∈0,4 is called logistic parameter. When x∈0,1, the logistic map is in a chaotic state.Figure9 depicts the possible values of x obtained after iteration for different values of μ for a certain value of the initial condition x0.Figure 9
Value range ofx.As can be seen in Figure9, it can be seen that the closer μ is to 4, the closer the value range of x is to the average distribution over the whole 0,1 interval; therefore, the logistic control parameter μ=4 is selected in this paper.The steps of initializing the population of logistic chaotic map are shown in Algorithm1. Firstly, an Nvar-dimensional vector x0 is randomly generated, and then Xpop vectors xxpop are generated by Equation (16). Finally, each component of xxpop is carried to the value interval xnvarL,xnvarU of the initialization boundary. The mathematical formulation of the carrier operation to generate the parent population is as follows.
(17)xParent0=xnvarL+xxpop⋅xnvarU−xnvarL.Algorithm 1: Logistic chaotic map initialization population.
Input: Nvar, Xpop, μ, xnvarL, xnvarUOutput: xParent01 randomly generate aNvar-dimensional x02 use Equation (16) generate xxpop3 use Equation (17) generate xParent0


### 3.4.2. Adaptive Mutating

After the initialization operation, the algorithm enters the mutation operation phase. For each target vectorxxpopGxpop=1,2,⋯,Xpop, the mutation strategy of the basic differential evolution algorithm can be described as “DE/rand/1,” and its mutation vector is generated by the following formula.
(18)vxpopG+1=xr1G+F⋅xr2G−xr3G.Among them,xr1G, xr2G, and xr3G are randomly selected mutually dissimilar individuals in the G-th generation population. r1, r2, and r3 are mutually dissimilar integers randomly selected in the interval 1,Xpop, and all are different from xpop. F is a real constant factor that controls the scaling of the deviation variable, usually F∈0,2.DE/rand/1 has many advantages such as simple structure, convenient, and fast using and has very powerful global search capability. The basic differential evolution algorithm takes real constants for the mutation operator in the search process, and the mutation operator is more difficult to determine in the simulation experiment. If the mutation rate is too large, the algorithm is inefficient in searching and the accuracy of the global optimal solution is low. If the variation rate is too small, the population diversity is reduced and the phenomenon of “early maturity” is easy to occur [26]. Therefore, a differential evolution algorithm with adaptive mutation operator is used in this paper.
(19)λ=e1−Gm/Gm+1−G,F=F0×2λ.Among them,F0 denotes the initial mutation operator, and F0=0.4 is taken in this paper. Gm denotes the maximum evolutionary algebra. In this paper, Gm=600. G denotes the current evolutionary algebra. When the algorithm starts searching, the mutation operator is 2F0, which has a large value, and the population maintains good individual diversity. As the algorithm searches, the mutation operator decreases gradually, and the mutation operator approaches F0 at the later stage, so that the population retains good information and thus avoids the optimal solution being destroyed.


### 3.4.3. Hybrid Crossover Operation

For each individualxxpopG of the population, the basic crossover operation applies binomial crossover to create the test vector uxpopG.
(20)uxpop,nvarG+1=vxpop,nvarG+1,ifRnvar≤CRornvar=nvarrand,xxpop,nvarG+1,ifRnvar>CRornvar≠nvarrand,xpop=1,2,⋯,Xpop;nvar=1,2,⋯,Nvar.Among them,Rnvar denotes a random number uniformly distributed between 0,1. nvarrand denotes a random integer in the range 1,Nvar, which guarantees that at least one gene is inherited from the parent individual. CR is the crossover probability, usually in the interval 0,1.The crossover operation of Equation (20) is expressed as “bin,” and its combination with the mutation operation of Equation (18) is expressed as “DE/rand/1/bin.” This algorithm is very popular in global optimization [27]. A random range crossover operator is used in this paper.
(21)CR=0.5×1+rand0,1.This way the mean value of the crossover operator is kept at 0.75, taking into account the possible random variation in the amplification of the difference vector, which is conducive to maintaining the population diversity.Simulated binary crossover (SBX) was designed by Deb and Agrawal [28], which is used to simulate the search pattern of single point binary crossover operator in real coded representation. It is an important recombination operator, which is applied to various multiobjective optimization algorithms. The mathematical expression is as follows.
(22)ur1G+1=0.5×1+β⋅xr1G+1−β⋅xr2G,ur2G+1=0.5×1−β⋅xr1G+1+β⋅xr2G,wherein ur1G+1 and ur2G+1 are child nodes generated by two parent nodes xr1G and xr2G. β follows the polynomial probability distribution and is calculated as follows.
(23)β=2×rand1/1+η,ifrand≤0.5,12−2×rand1/1+η,otherwise.Among them, the distribution indexη is a predefined nonnegative real number, and in general, the value of η in SBX is usually taken as 15 or 20. In this paper, η takes 20. At this time, SBX has a strong local search capability around the parent [29]. The steps of SBX operation are shown in Algorithm 2.Algorithm 2: SBX operator for the i-th gene bit.
Input: xr1G,xr2G,ηOutput: ur1G+1,ur2G+11 use Equation (23) to generate β2 use Equation (22) to generate ur1G+1, ur2G+1Considering that the binomial crossover strategy has a good global search ability and SBX has a strong parent-centered local search ability, this paper mixes binomial crossover with simulated binary crossover and sets a value ofχ that keeps changing with the number of population iterations to balance the global search ability and local search ability of the algorithm. Let χ=1−0.6×G/Gm. After the algorithm starts to search, before the iteration proceeds to Gm/2, because the value of χ is larger, the crossover strategy of Equation (20) is chosen at this time, and the algorithm has a strong global search ability. Then, the algorithm iteration proceeds to the second half, because the value of χ decreases, the crossover strategy of Equation (22) is chosen, and the algorithm has a strong local search ability. The hybrid crossover strategy that combines the binomial crossover strategy and the simulated binary crossover strategy is expressed as “bin-SBX.” The pseudocode of “bin-SBX” hybrid crossover strategy is shown in Algorithm 3.Algorithm 3: Bin-SBX.
Input: xr1G,xr2G,η,vxpop,nvarG+1,xxpop,nvarG+1Output: xChildG+11 use Equation (23) to generate β2 use Equation (22) to generate u′xpop,nvarG+13 ifχ≥0.74  use Equation (20) to generate xChildG+1=uxpop,nvarG+15  else ifχ<0.76xChildG+1=u′xpop,nvarG+17  end8 end


### 3.4.4. Selection Operation

Using the basic selection strategy of differential evolution, the test vectors are compared with the target vectors in the current population according to the greedy criterion, and the vectors with smaller objective function values are retained. All individuals in the next generation are better or at least as good as their corresponding individuals in the current population.


## 3.4.1. Chaos Initializing

Differential evolution algorithm usesXpop real-valued parameters vectors with Nvar dimension to take them as the population of each generation, and each individual is denoted as follows.
(15)xxpopGxpop=1,2,⋯,Xpop,where xpop denotes the sequence of individuals in the population. G denotes evolutionary algebra. Xpop denotes the population size. In this paper, Xpop=200 is taken, and Xpop remains unchanged in the process of minimization. Nvar denotes the number of decision variables, which is Nvar=2 in this paper.In DE algorithm, a common method of population initializing is to randomly select from the values within the given boundary constraints. The convergence speed and efficiency of the algorithm’s search process are affected by the degree of scattering of the initial positions and the uniformity of their positions in the search space. When the solution space is large and the number of individuals in the population is not large enough, it may cause uneven distribution of individuals in initializing population. And chaos is a kind of irregular motion, which comes from nonlinear dynamical system, and its motion state has good ergodicity and randomness. In this paper, the population is initialized by a one-dimensional logistic chaotic map, which has extremely complex dynamical behavior. Its mathematical expression is as follows.(16)xi+1=μ⋅xi⋅1−xi.Among them,μ∈0,4 is called logistic parameter. When x∈0,1, the logistic map is in a chaotic state.Figure9 depicts the possible values of x obtained after iteration for different values of μ for a certain value of the initial condition x0.Figure 9
Value range ofx.As can be seen in Figure9, it can be seen that the closer μ is to 4, the closer the value range of x is to the average distribution over the whole 0,1 interval; therefore, the logistic control parameter μ=4 is selected in this paper.The steps of initializing the population of logistic chaotic map are shown in Algorithm1. Firstly, an Nvar-dimensional vector x0 is randomly generated, and then Xpop vectors xxpop are generated by Equation (16). Finally, each component of xxpop is carried to the value interval xnvarL,xnvarU of the initialization boundary. The mathematical formulation of the carrier operation to generate the parent population is as follows.
(17)xParent0=xnvarL+xxpop⋅xnvarU−xnvarL.Algorithm 1: Logistic chaotic map initialization population.
Input: Nvar, Xpop, μ, xnvarL, xnvarUOutput: xParent01 randomly generate aNvar-dimensional x02 use Equation (16) generate xxpop3 use Equation (17) generate xParent0


## 3.4.2. Adaptive Mutating

After the initialization operation, the algorithm enters the mutation operation phase. For each target vectorxxpopGxpop=1,2,⋯,Xpop, the mutation strategy of the basic differential evolution algorithm can be described as “DE/rand/1,” and its mutation vector is generated by the following formula.
(18)vxpopG+1=xr1G+F⋅xr2G−xr3G.Among them,xr1G, xr2G, and xr3G are randomly selected mutually dissimilar individuals in the G-th generation population. r1, r2, and r3 are mutually dissimilar integers randomly selected in the interval 1,Xpop, and all are different from xpop. F is a real constant factor that controls the scaling of the deviation variable, usually F∈0,2.DE/rand/1 has many advantages such as simple structure, convenient, and fast using and has very powerful global search capability. The basic differential evolution algorithm takes real constants for the mutation operator in the search process, and the mutation operator is more difficult to determine in the simulation experiment. If the mutation rate is too large, the algorithm is inefficient in searching and the accuracy of the global optimal solution is low. If the variation rate is too small, the population diversity is reduced and the phenomenon of “early maturity” is easy to occur [26]. Therefore, a differential evolution algorithm with adaptive mutation operator is used in this paper.
(19)λ=e1−Gm/Gm+1−G,F=F0×2λ.Among them,F0 denotes the initial mutation operator, and F0=0.4 is taken in this paper. Gm denotes the maximum evolutionary algebra. In this paper, Gm=600. G denotes the current evolutionary algebra. When the algorithm starts searching, the mutation operator is 2F0, which has a large value, and the population maintains good individual diversity. As the algorithm searches, the mutation operator decreases gradually, and the mutation operator approaches F0 at the later stage, so that the population retains good information and thus avoids the optimal solution being destroyed.


## 3.4.3. Hybrid Crossover Operation

For each individualxxpopG of the population, the basic crossover operation applies binomial crossover to create the test vector uxpopG.
(20)uxpop,nvarG+1=vxpop,nvarG+1,ifRnvar≤CRornvar=nvarrand,xxpop,nvarG+1,ifRnvar>CRornvar≠nvarrand,xpop=1,2,⋯,Xpop;nvar=1,2,⋯,Nvar.Among them,Rnvar denotes a random number uniformly distributed between 0,1. nvarrand denotes a random integer in the range 1,Nvar, which guarantees that at least one gene is inherited from the parent individual. CR is the crossover probability, usually in the interval 0,1.The crossover operation of Equation (20) is expressed as “bin,” and its combination with the mutation operation of Equation (18) is expressed as “DE/rand/1/bin.” This algorithm is very popular in global optimization [27]. A random range crossover operator is used in this paper.
(21)CR=0.5×1+rand0,1.This way the mean value of the crossover operator is kept at 0.75, taking into account the possible random variation in the amplification of the difference vector, which is conducive to maintaining the population diversity.Simulated binary crossover (SBX) was designed by Deb and Agrawal [28], which is used to simulate the search pattern of single point binary crossover operator in real coded representation. It is an important recombination operator, which is applied to various multiobjective optimization algorithms. The mathematical expression is as follows.
(22)ur1G+1=0.5×1+β⋅xr1G+1−β⋅xr2G,ur2G+1=0.5×1−β⋅xr1G+1+β⋅xr2G,wherein ur1G+1 and ur2G+1 are child nodes generated by two parent nodes xr1G and xr2G. β follows the polynomial probability distribution and is calculated as follows.
(23)β=2×rand1/1+η,ifrand≤0.5,12−2×rand1/1+η,otherwise.Among them, the distribution indexη is a predefined nonnegative real number, and in general, the value of η in SBX is usually taken as 15 or 20. In this paper, η takes 20. At this time, SBX has a strong local search capability around the parent [29]. The steps of SBX operation are shown in Algorithm 2.Algorithm 2: SBX operator for the i-th gene bit.
Input: xr1G,xr2G,ηOutput: ur1G+1,ur2G+11 use Equation (23) to generate β2 use Equation (22) to generate ur1G+1, ur2G+1Considering that the binomial crossover strategy has a good global search ability and SBX has a strong parent-centered local search ability, this paper mixes binomial crossover with simulated binary crossover and sets a value ofχ that keeps changing with the number of population iterations to balance the global search ability and local search ability of the algorithm. Let χ=1−0.6×G/Gm. After the algorithm starts to search, before the iteration proceeds to Gm/2, because the value of χ is larger, the crossover strategy of Equation (20) is chosen at this time, and the algorithm has a strong global search ability. Then, the algorithm iteration proceeds to the second half, because the value of χ decreases, the crossover strategy of Equation (22) is chosen, and the algorithm has a strong local search ability. The hybrid crossover strategy that combines the binomial crossover strategy and the simulated binary crossover strategy is expressed as “bin-SBX.” The pseudocode of “bin-SBX” hybrid crossover strategy is shown in Algorithm 3.Algorithm 3: Bin-SBX.
Input: xr1G,xr2G,η,vxpop,nvarG+1,xxpop,nvarG+1Output: xChildG+11 use Equation (23) to generate β2 use Equation (22) to generate u′xpop,nvarG+13 ifχ≥0.74  use Equation (20) to generate xChildG+1=uxpop,nvarG+15  else ifχ<0.76xChildG+1=u′xpop,nvarG+17  end8 end


## 3.4.4. Selection Operation

Using the basic selection strategy of differential evolution, the test vectors are compared with the target vectors in the current population according to the greedy criterion, and the vectors with smaller objective function values are retained. All individuals in the next generation are better or at least as good as their corresponding individuals in the current population.


## 4. Simulation and Analysis

### 4.1. Evaluation Index

The multiobjective evolutionary algorithm has two performance evaluation criteria: (1) the degree to which the obtained solution set approaches the Pareto optimal solution set and (2) whether the obtained solution set has a good diversity distribution. The two criteria proposed based on this are [30] the convergence criterion r and the diversity criterion Δ. Assume that the corresponding set of the solution set obtained by the multiobjective evolution algorithm in the target space is F, and the theoretically optimal Pareto front is F′.Definition 4 (convergence criterionr).
The convergence criterionrF,F′ is used to measure the degree to which F approaches F′.
(24)rF,F′=1F∑f∈Fminf−f′,f′∈F′.
In this paper,F=f1,f2,f3, and F′=f1′,f2′,f3′. f1=f11,f12,⋯,f1nT, f2=f21,f22,⋯,f2nT, and f3=f31,f32,⋯,f3nT, where n is the number of Pareto front obtained by the multiobjective evolutionary algorithm, also known as the number of Pareto set. f1′=0, f2′=0, and f3′=0.Obviously, the smaller the value ofrF,F′, the closer F is to F′. In particular, when F⊆F′, rF,F′=0.Definition 5 (diversity criterionΔ).
The diversity criterionΔF,F′ is used to measure the degree of diversity distribution of F on F′.
(25)ΔF,F′=df−dl+∑i=1F−1di−d¯df+dl+F−1d¯,where di is the Euclidean distance between two consecutive solutions in the objective space, d¯ is the arithmetic mean of all distances dii=1,2,⋯,F−1, and df and dl denote the Euclidean distance between the theoretically optimal extremal solution f1′,f2′,f3′ and the boundary solutions f11,f21,f31 and f1n,f2n,f3n of the obtained nondominated set. Obviously, the smaller the value of ΔF,F′, the better the diversity distribution of F. When F is perfectly uniformly distributed and contains the extreme target vectors in F′, ΔF,F′=0.


### 4.2. Simulation Result

This paper creates a simulation environment as a verification tool in which the spinning projectile model can withstand different winds. Two hypotheses are made. Hypothesis 1 is that the projectile is subjected to a constant wind velocity magnitude throughout the flight. Hypothesis 2 is that the projectile is subjected to a pulsating wind velocity that varies uniformly with time throughout the flight. In the actual flight test, the wind velocity is estimated by other methods without wind velocity measurement, and the approximate range of wind speed is derived. In the simulation experiment, this idea is reflected in the selection of optimization boundary. Although the magnitude of these wind velocity is known in the simulation test, it can be assumed that the size of wind velocity can be identified by selecting the appropriate optimization boundary on the premise of unknown wind velocity magnitude.


#### 4.2.1. Constant Wind Simulation

When the projectile is subjected to a constant wind velocity throughout the flight, that is, given constant longitudinal and cross wind values, the two components are identified. In this paper, three groups of constant wind simulation tests are done, and each group of tests is simulated three times with different optimization boundaries. After the multiobjective optimization completes the maximum number of iterations, it provides a set of solution sets that approximate the Pareto setΘ∗. Table 1 shows some of the values extracted from the solution sets. The second column of the table gives the different optimization boundaries set for each test, and the third column gives the number of solutions obtained that approximate the Pareto set. The fourth column gives the average value of the solutions in the solution set obtained for each test, which is a dichotomous vector expressed as the wind velocity W=Wx,Wz, where both Wx and Wz have units of m/s. The values in the fifth column indicate the standard deviation of the whole set. The sixth and seventh columns give the absolute and relative errors of each component of the wind estimation.Table 1
Identification results of constant wind velocity components.
Wind valueOptimization boundaryΘ∗Average valueStandard deviationAbsolute errorRelative errorWx=51,10134.9800795.9993180.4707190.6519270.0199210.0006820.0039840.0001141,15135.1775135.8492370.8645500.4747730.1775130.1507630.0355030.025127Wz=61,20105.3050365.8851570.7088900.5164410.3050360.1148430.0610070.019141Wx=71,1057.3677442.9800410.3901030.1774900.3677440.0199590.0525350.0066531,1526.9705503.1655210.0436410.2077090.0294500.1655210.0042070.055174Wz=31,2056.9725992.8328620.4430000.4849010.0274010.1671380.0039140.055713Wx=95,1548.85040912.1674060.1843690.3332370.1495910.1674060.0166210.0139515,2039.41738012.0448031.1832650.0072120.4173800.0448030.0463760.003734Wz=125,2579.48925412.1671071.0277960.5019320.4892540.1671070.0543620.013926Table1 shows that the identification process converges to the actual wind with a small error. Empirically, if the obtained identification results are within ±10% error of the set wind values, it can be proved that this method is effective for the estimation of constant wind. The results in the fourth column of the table show that the average values obtained from these three sets of tests are within the acceptable error of the set wind values; then, the validity of the method is preliminarily verified. Usually, even for constant wind and known model structure, the estimated winds are not unique, but very close to the point set of real winds. There are two reasons that may lead to this situation. First, the optimization process has not yet fully converged to the optimal value. Second, the size of the optimization boundary interval of the search space also affects the identification results. It can be seen from the second and third group of tests that when the expanded optimization boundary interval reaches a certain value, the standard deviation of the identified point set may increase, which means that there are some points in the point set after expanding the optimization boundary that deviate greatly from the real values.Considering the convergence and diversity criteria of the method for identifying constant wind, the values ofrF,F′ and ΔF,F′ for three groups of simulation tests are calculated according to Equations (24) and (25) as shown in Table 2.Table 2
Convergence and diversity of Pareto front obtained from constant wind simulation.
Wind valueOptimization boundaryΘ∗rF,F′ΔF,F′Wx=51,10133.21×10−110.6319751,15131.10×10−100.820344Wz=61,20108.56×10−110.927726Wx=71,1054.15×10−110.0127411,1522.62×10−130.367500Wz=31,2052.31×10−110.511674Wx=95,1547.63×10−120.1668035,2031.55×10−100.900160Wz=125,2571.85×10−100.637612Table2 shows that in terms of convergence, in general, the smaller the optimization boundary, the better the convergence of Pareto front. In particular, from the first group and the second group of experiments, Pareto front has better convergence when a larger optimization boundary optimization obtains fewer Pareto set. In terms of diversity, in general, the smaller the optimization boundary, the better the Pareto front diversity distribution. In particular, from the third group of experiments, when the more Pareto set are optimized, the better the Pareto front diversity distribution. Therefore, the convergence and diversity distribution of Pareto front are related to the size of optimization boundary and the number of Pareto set.


#### 4.2.2. Variable Wind Simulation

When the projectile is subjected to a pulsating wind velocity that varies uniformly with time throughout the flight, the two components are identified in this paper given the longitudinal and cross wind values that vary according to the triangular function. The wind is modeled as follows.(26)Wx=5+sin2π/10t,Wz=6+cos2π/10t,where the units of Wx and Wz are also m/s.According to Section2.3, it can be seen that the whole flight simulation test lasts about 90 s. In the tests of this section, the variation period of wind is T=10s. The tests are taken for three time periods for simulation, namely, 0−1s (T/10), 0−2.5s (T/4), and 0−5s (T/2), and it is assumed in this paper that the wind velocity magnitude is uniformly constant in a short time period. The same optimization boundaries are taken as in Section 4.2.1 for Wx=5 and Wz=6 for nine simulations. Table 3 shows some values extracted from the solution set. Table 4 shows the average values of the three time periods under this pulsating wind velocity.Table 3
Identification results of the pulsating wind velocity component.
Time periodOptimization boundaryΘ∗Average valueStandard deviationAbsolute errorRelative error0−1s1,1055.5142996.9295870.6084351.3676920.2074020.0049440.0390820.0007131,1575.5233116.7969770.6029522.9006570.2164140.1375540.0407800.0198361,2035.0436117.0436441.8985690.4540540.2632860.1091130.0496120.0157351,1015.5874766.839677000.0491440.2030570.0087190.0305960−2.5s1,1525.3168666.5199990.2674490.5353680.3197540.1166210.0567280.0175721,2025.4750846.3565390.1406660.7976170.1615360.2800810.0286580.0422021,1045.6405715.5928111.6345902.0476030.0039510.4071890.0007010.0678650−5s1,1555.8470305.9995540.6888612.7926120.2104100.0004460.0373290.0000741,2085.5044575.7641770.7020403.5128660.1321630.2358230.0234470.039304Table 4
Mean values of the pulsating wind velocity function.
Time periodMean values0−1s5.3068976.9345310−2.5s5.6366206.6366200−5s5.6366206.000000The comparison between Tables3 and 4 shows that the algorithm still has a good identification accuracy (±10%) when the optimization boundary becomes larger and larger. Of course, it is worth mentioning that, similar to the constant wind simulation, the standard deviation of the Pareto set identified by the algorithm becomes larger with the expansion of the optimization boundary, which means that the accuracy of the algorithm in searching for the optimal solution is declining, but it can still meet the accuracy requirement of the identification of the average wind velocity magnitude.


## 4.1. Evaluation Index

The multiobjective evolutionary algorithm has two performance evaluation criteria: (1) the degree to which the obtained solution set approaches the Pareto optimal solution set and (2) whether the obtained solution set has a good diversity distribution. The two criteria proposed based on this are [30] the convergence criterion r and the diversity criterion Δ. Assume that the corresponding set of the solution set obtained by the multiobjective evolution algorithm in the target space is F, and the theoretically optimal Pareto front is F′.Definition 4 (convergence criterionr).
The convergence criterionrF,F′ is used to measure the degree to which F approaches F′.
(24)rF,F′=1F∑f∈Fminf−f′,f′∈F′.
In this paper,F=f1,f2,f3, and F′=f1′,f2′,f3′. f1=f11,f12,⋯,f1nT, f2=f21,f22,⋯,f2nT, and f3=f31,f32,⋯,f3nT, where n is the number of Pareto front obtained by the multiobjective evolutionary algorithm, also known as the number of Pareto set. f1′=0, f2′=0, and f3′=0.Obviously, the smaller the value ofrF,F′, the closer F is to F′. In particular, when F⊆F′, rF,F′=0.Definition 5 (diversity criterionΔ).
The diversity criterionΔF,F′ is used to measure the degree of diversity distribution of F on F′.
(25)ΔF,F′=df−dl+∑i=1F−1di−d¯df+dl+F−1d¯,where di is the Euclidean distance between two consecutive solutions in the objective space, d¯ is the arithmetic mean of all distances dii=1,2,⋯,F−1, and df and dl denote the Euclidean distance between the theoretically optimal extremal solution f1′,f2′,f3′ and the boundary solutions f11,f21,f31 and f1n,f2n,f3n of the obtained nondominated set. Obviously, the smaller the value of ΔF,F′, the better the diversity distribution of F. When F is perfectly uniformly distributed and contains the extreme target vectors in F′, ΔF,F′=0.


## 4.2. Simulation Result

This paper creates a simulation environment as a verification tool in which the spinning projectile model can withstand different winds. Two hypotheses are made. Hypothesis 1 is that the projectile is subjected to a constant wind velocity magnitude throughout the flight. Hypothesis 2 is that the projectile is subjected to a pulsating wind velocity that varies uniformly with time throughout the flight. In the actual flight test, the wind velocity is estimated by other methods without wind velocity measurement, and the approximate range of wind speed is derived. In the simulation experiment, this idea is reflected in the selection of optimization boundary. Although the magnitude of these wind velocity is known in the simulation test, it can be assumed that the size of wind velocity can be identified by selecting the appropriate optimization boundary on the premise of unknown wind velocity magnitude.


### 4.2.1. Constant Wind Simulation

When the projectile is subjected to a constant wind velocity throughout the flight, that is, given constant longitudinal and cross wind values, the two components are identified. In this paper, three groups of constant wind simulation tests are done, and each group of tests is simulated three times with different optimization boundaries. After the multiobjective optimization completes the maximum number of iterations, it provides a set of solution sets that approximate the Pareto setΘ∗. Table 1 shows some of the values extracted from the solution sets. The second column of the table gives the different optimization boundaries set for each test, and the third column gives the number of solutions obtained that approximate the Pareto set. The fourth column gives the average value of the solutions in the solution set obtained for each test, which is a dichotomous vector expressed as the wind velocity W=Wx,Wz, where both Wx and Wz have units of m/s. The values in the fifth column indicate the standard deviation of the whole set. The sixth and seventh columns give the absolute and relative errors of each component of the wind estimation.Table 1
Identification results of constant wind velocity components.
Wind valueOptimization boundaryΘ∗Average valueStandard deviationAbsolute errorRelative errorWx=51,10134.9800795.9993180.4707190.6519270.0199210.0006820.0039840.0001141,15135.1775135.8492370.8645500.4747730.1775130.1507630.0355030.025127Wz=61,20105.3050365.8851570.7088900.5164410.3050360.1148430.0610070.019141Wx=71,1057.3677442.9800410.3901030.1774900.3677440.0199590.0525350.0066531,1526.9705503.1655210.0436410.2077090.0294500.1655210.0042070.055174Wz=31,2056.9725992.8328620.4430000.4849010.0274010.1671380.0039140.055713Wx=95,1548.85040912.1674060.1843690.3332370.1495910.1674060.0166210.0139515,2039.41738012.0448031.1832650.0072120.4173800.0448030.0463760.003734Wz=125,2579.48925412.1671071.0277960.5019320.4892540.1671070.0543620.013926Table1 shows that the identification process converges to the actual wind with a small error. Empirically, if the obtained identification results are within ±10% error of the set wind values, it can be proved that this method is effective for the estimation of constant wind. The results in the fourth column of the table show that the average values obtained from these three sets of tests are within the acceptable error of the set wind values; then, the validity of the method is preliminarily verified. Usually, even for constant wind and known model structure, the estimated winds are not unique, but very close to the point set of real winds. There are two reasons that may lead to this situation. First, the optimization process has not yet fully converged to the optimal value. Second, the size of the optimization boundary interval of the search space also affects the identification results. It can be seen from the second and third group of tests that when the expanded optimization boundary interval reaches a certain value, the standard deviation of the identified point set may increase, which means that there are some points in the point set after expanding the optimization boundary that deviate greatly from the real values.Considering the convergence and diversity criteria of the method for identifying constant wind, the values ofrF,F′ and ΔF,F′ for three groups of simulation tests are calculated according to Equations (24) and (25) as shown in Table 2.Table 2
Convergence and diversity of Pareto front obtained from constant wind simulation.
Wind valueOptimization boundaryΘ∗rF,F′ΔF,F′Wx=51,10133.21×10−110.6319751,15131.10×10−100.820344Wz=61,20108.56×10−110.927726Wx=71,1054.15×10−110.0127411,1522.62×10−130.367500Wz=31,2052.31×10−110.511674Wx=95,1547.63×10−120.1668035,2031.55×10−100.900160Wz=125,2571.85×10−100.637612Table2 shows that in terms of convergence, in general, the smaller the optimization boundary, the better the convergence of Pareto front. In particular, from the first group and the second group of experiments, Pareto front has better convergence when a larger optimization boundary optimization obtains fewer Pareto set. In terms of diversity, in general, the smaller the optimization boundary, the better the Pareto front diversity distribution. In particular, from the third group of experiments, when the more Pareto set are optimized, the better the Pareto front diversity distribution. Therefore, the convergence and diversity distribution of Pareto front are related to the size of optimization boundary and the number of Pareto set.


### 4.2.2. Variable Wind Simulation

When the projectile is subjected to a pulsating wind velocity that varies uniformly with time throughout the flight, the two components are identified in this paper given the longitudinal and cross wind values that vary according to the triangular function. The wind is modeled as follows.(26)Wx=5+sin2π/10t,Wz=6+cos2π/10t,where the units of Wx and Wz are also m/s.According to Section2.3, it can be seen that the whole flight simulation test lasts about 90 s. In the tests of this section, the variation period of wind is T=10s. The tests are taken for three time periods for simulation, namely, 0−1s (T/10), 0−2.5s (T/4), and 0−5s (T/2), and it is assumed in this paper that the wind velocity magnitude is uniformly constant in a short time period. The same optimization boundaries are taken as in Section 4.2.1 for Wx=5 and Wz=6 for nine simulations. Table 3 shows some values extracted from the solution set. Table 4 shows the average values of the three time periods under this pulsating wind velocity.Table 3
Identification results of the pulsating wind velocity component.
Time periodOptimization boundaryΘ∗Average valueStandard deviationAbsolute errorRelative error0−1s1,1055.5142996.9295870.6084351.3676920.2074020.0049440.0390820.0007131,1575.5233116.7969770.6029522.9006570.2164140.1375540.0407800.0198361,2035.0436117.0436441.8985690.4540540.2632860.1091130.0496120.0157351,1015.5874766.839677000.0491440.2030570.0087190.0305960−2.5s1,1525.3168666.5199990.2674490.5353680.3197540.1166210.0567280.0175721,2025.4750846.3565390.1406660.7976170.1615360.2800810.0286580.0422021,1045.6405715.5928111.6345902.0476030.0039510.4071890.0007010.0678650−5s1,1555.8470305.9995540.6888612.7926120.2104100.0004460.0373290.0000741,2085.5044575.7641770.7020403.5128660.1321630.2358230.0234470.039304Table 4
Mean values of the pulsating wind velocity function.
Time periodMean values0−1s5.3068976.9345310−2.5s5.6366206.6366200−5s5.6366206.000000The comparison between Tables3 and 4 shows that the algorithm still has a good identification accuracy (±10%) when the optimization boundary becomes larger and larger. Of course, it is worth mentioning that, similar to the constant wind simulation, the standard deviation of the Pareto set identified by the algorithm becomes larger with the expansion of the optimization boundary, which means that the accuracy of the algorithm in searching for the optimal solution is declining, but it can still meet the accuracy requirement of the identification of the average wind velocity magnitude.


## 4.2.1. Constant Wind Simulation

When the projectile is subjected to a constant wind velocity throughout the flight, that is, given constant longitudinal and cross wind values, the two components are identified. In this paper, three groups of constant wind simulation tests are done, and each group of tests is simulated three times with different optimization boundaries. After the multiobjective optimization completes the maximum number of iterations, it provides a set of solution sets that approximate the Pareto setΘ∗. Table 1 shows some of the values extracted from the solution sets. The second column of the table gives the different optimization boundaries set for each test, and the third column gives the number of solutions obtained that approximate the Pareto set. The fourth column gives the average value of the solutions in the solution set obtained for each test, which is a dichotomous vector expressed as the wind velocity W=Wx,Wz, where both Wx and Wz have units of m/s. The values in the fifth column indicate the standard deviation of the whole set. The sixth and seventh columns give the absolute and relative errors of each component of the wind estimation.Table 1
Identification results of constant wind velocity components.
Wind valueOptimization boundaryΘ∗Average valueStandard deviationAbsolute errorRelative errorWx=51,10134.9800795.9993180.4707190.6519270.0199210.0006820.0039840.0001141,15135.1775135.8492370.8645500.4747730.1775130.1507630.0355030.025127Wz=61,20105.3050365.8851570.7088900.5164410.3050360.1148430.0610070.019141Wx=71,1057.3677442.9800410.3901030.1774900.3677440.0199590.0525350.0066531,1526.9705503.1655210.0436410.2077090.0294500.1655210.0042070.055174Wz=31,2056.9725992.8328620.4430000.4849010.0274010.1671380.0039140.055713Wx=95,1548.85040912.1674060.1843690.3332370.1495910.1674060.0166210.0139515,2039.41738012.0448031.1832650.0072120.4173800.0448030.0463760.003734Wz=125,2579.48925412.1671071.0277960.5019320.4892540.1671070.0543620.013926Table1 shows that the identification process converges to the actual wind with a small error. Empirically, if the obtained identification results are within ±10% error of the set wind values, it can be proved that this method is effective for the estimation of constant wind. The results in the fourth column of the table show that the average values obtained from these three sets of tests are within the acceptable error of the set wind values; then, the validity of the method is preliminarily verified. Usually, even for constant wind and known model structure, the estimated winds are not unique, but very close to the point set of real winds. There are two reasons that may lead to this situation. First, the optimization process has not yet fully converged to the optimal value. Second, the size of the optimization boundary interval of the search space also affects the identification results. It can be seen from the second and third group of tests that when the expanded optimization boundary interval reaches a certain value, the standard deviation of the identified point set may increase, which means that there are some points in the point set after expanding the optimization boundary that deviate greatly from the real values.Considering the convergence and diversity criteria of the method for identifying constant wind, the values ofrF,F′ and ΔF,F′ for three groups of simulation tests are calculated according to Equations (24) and (25) as shown in Table 2.Table 2
Convergence and diversity of Pareto front obtained from constant wind simulation.
Wind valueOptimization boundaryΘ∗rF,F′ΔF,F′Wx=51,10133.21×10−110.6319751,15131.10×10−100.820344Wz=61,20108.56×10−110.927726Wx=71,1054.15×10−110.0127411,1522.62×10−130.367500Wz=31,2052.31×10−110.511674Wx=95,1547.63×10−120.1668035,2031.55×10−100.900160Wz=125,2571.85×10−100.637612Table2 shows that in terms of convergence, in general, the smaller the optimization boundary, the better the convergence of Pareto front. In particular, from the first group and the second group of experiments, Pareto front has better convergence when a larger optimization boundary optimization obtains fewer Pareto set. In terms of diversity, in general, the smaller the optimization boundary, the better the Pareto front diversity distribution. In particular, from the third group of experiments, when the more Pareto set are optimized, the better the Pareto front diversity distribution. Therefore, the convergence and diversity distribution of Pareto front are related to the size of optimization boundary and the number of Pareto set.


## 4.2.2. Variable Wind Simulation

When the projectile is subjected to a pulsating wind velocity that varies uniformly with time throughout the flight, the two components are identified in this paper given the longitudinal and cross wind values that vary according to the triangular function. The wind is modeled as follows.(26)Wx=5+sin2π/10t,Wz=6+cos2π/10t,where the units of Wx and Wz are also m/s.According to Section2.3, it can be seen that the whole flight simulation test lasts about 90 s. In the tests of this section, the variation period of wind is T=10s. The tests are taken for three time periods for simulation, namely, 0−1s (T/10), 0−2.5s (T/4), and 0−5s (T/2), and it is assumed in this paper that the wind velocity magnitude is uniformly constant in a short time period. The same optimization boundaries are taken as in Section 4.2.1 for Wx=5 and Wz=6 for nine simulations. Table 3 shows some values extracted from the solution set. Table 4 shows the average values of the three time periods under this pulsating wind velocity.Table 3
Identification results of the pulsating wind velocity component.
Time periodOptimization boundaryΘ∗Average valueStandard deviationAbsolute errorRelative error0−1s1,1055.5142996.9295870.6084351.3676920.2074020.0049440.0390820.0007131,1575.5233116.7969770.6029522.9006570.2164140.1375540.0407800.0198361,2035.0436117.0436441.8985690.4540540.2632860.1091130.0496120.0157351,1015.5874766.839677000.0491440.2030570.0087190.0305960−2.5s1,1525.3168666.5199990.2674490.5353680.3197540.1166210.0567280.0175721,2025.4750846.3565390.1406660.7976170.1615360.2800810.0286580.0422021,1045.6405715.5928111.6345902.0476030.0039510.4071890.0007010.0678650−5s1,1555.8470305.9995540.6888612.7926120.2104100.0004460.0373290.0000741,2085.5044575.7641770.7020403.5128660.1321630.2358230.0234470.039304Table 4
Mean values of the pulsating wind velocity function.
Time periodMean values0−1s5.3068976.9345310−2.5s5.6366206.6366200−5s5.6366206.000000The comparison between Tables3 and 4 shows that the algorithm still has a good identification accuracy (±10%) when the optimization boundary becomes larger and larger. Of course, it is worth mentioning that, similar to the constant wind simulation, the standard deviation of the Pareto set identified by the algorithm becomes larger with the expansion of the optimization boundary, which means that the accuracy of the algorithm in searching for the optimal solution is declining, but it can still meet the accuracy requirement of the identification of the average wind velocity magnitude.


## 5. Conclusions

In this paper, a multiobjective chaotic adaptive differential evolutionary algorithm is proposed to identify the longitudinal and cross wind values during the flight of the projectile based on the 4D kinematic model of spinning projectile. After conducting several sets of constant wind simulation tests and variable wind simulation tests, the errors of the identification results are within±10%, which proves the effectiveness of the algorithm. At the same time, the convergence and diversity indexes in constant wind simulation are calculated, and the convergence and diversity indexes of test results are analyzed qualitatively and quantitatively. In addition, the successful application of the algorithm in the 4D kinematic model of spinning projectile provides a basis for the identification of aerodynamic parameters to be performed later.



---

*Source: 9726595-2022-09-12.xml*