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.
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 . 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 . 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. Figures 1–3 simulate the dynamic process of a spinning projectile in actual flight.
Figure 1 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 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 , which exceeded the actual wind value by .
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 Section 2, 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 , 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 .
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 . Projecting it onto the ground coordinate system , the calculation formula is as follows.
Among them,
Furthermore, the modified mass-point ballistic equation established in the ground coordinate system is as follows.
In Equations (1)-(6), is the strength of the gravitational field near the earth’s surface, is the total mass of the projectile, is the characteristic area of the projectile, which is usually selected as the maximum cross-sectional area of the projectile, is the maximum cross-sectional diameter of the projectile, is the characteristic length of the projectile, and is the atmospheric density. is the roll angle and is the roll angle speed. is the velocity of the projectile relative to the air, , , and are the three components of the projectile velocity relative to the ground coordinate system , and are the axis and axis components of the wind speed relative to the ground coordinate system , , , and are the three position coordinates of the projectile relative to the ground coordinate system , and , , is the dynamic balance angle vector ’s three components relative to the ground coordinate system ; is the axial moment of inertia of the projectile. and are the zero lift drag coefficient and induced drag coefficient, respectively, is the drag coefficient of the projectile, is the derivative of the lift coefficient of the projectile, and is the joint partial derivative of the Magnus force coefficient of the projectile; is the derivative of extreme damping moment coefficient, is the derivative of overturning moment coefficient, and 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.
Among them,
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., and ), constant wind (e.g., when and ), and variable wind (when and ) conditions is shown in Figure 4.
The derivative of the lift coefficient of the projectile as a function of time in the no wind (i.e., and ), constant wind (e.g., when and ), and variable wind (when and ) conditions is shown in Figure 5.
The joint partial derivatives of the Magnus force coefficient of the projectile with time in the no wind (i.e., and ), constant wind (e.g., when and ), and variable wind (when and ) conditions is shown in Figure 6.
As can be seen from Figures 4–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 Section 2.3 that when the values of and 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 where is the solution that simultaneously minimizes cost functions . 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 .
Definition 1 (Pareto optimality [8]). An objective vector is Pareto optimal if there exist no other objective vectors such that for all there is , and for at least one there is .
Definition 2 (strict dominance [8]). An objective vector is dominated by another objective vector if for all .
Definition 3 (dominance [8]). An objective vector dominates another vector if is no worse than in all objectives and is better in at least one objective, i.e., .
Figure 7 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.
3.2. Definition of Multiobjective Optimization Problem
The 4D aerodynamic model of spinning projectile in Section 2.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.
In Equations (12)–(14), , , and are the calculated values provided by the aerodynamic model for , , and , respectively. These three cost functions constitute the objective space, and the two wind velocity components and 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 Figure 8.
3.4. Chaotic Adaptive Differential Evolution Process
3.4.1. Chaos Initializing
Differential evolution algorithm uses real-valued parameters vectors with dimension to take them as the population of each generation, and each individual is denoted as follows. where denotes the sequence of individuals in the population. denotes evolutionary algebra. denotes the population size. In this paper, is taken, and remains unchanged in the process of minimization. denotes the number of decision variables, which is 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.
Among them, is called logistic parameter. When , the logistic map is in a chaotic state.
Figure 9 depicts the possible values of obtained after iteration for different values of for a certain value of the initial condition .
As can be seen in Figure 9, it can be seen that the closer is to , the closer the value range of is to the average distribution over the whole interval; therefore, the logistic control parameter is selected in this paper.
The steps of initializing the population of logistic chaotic map are shown in Algorithm 1. Firstly, an -dimensional vector is randomly generated, and then vectors are generated by Equation (16). Finally, each component of is carried to the value interval of the initialization boundary. The mathematical formulation of the carrier operation to generate the parent population is as follows.
3.4.2. Adaptive Mutating
After the initialization operation, the algorithm enters the mutation operation phase. For each target vector , 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.
Among them, , , and are randomly selected mutually dissimilar individuals in the -th generation population. , , and are mutually dissimilar integers randomly selected in the interval , and all are different from . is a real constant factor that controls the scaling of the deviation variable, usually .
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.
Among them, denotes the initial mutation operator, and is taken in this paper. denotes the maximum evolutionary algebra. In this paper, . denotes the current evolutionary algebra. When the algorithm starts searching, the mutation operator is , 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 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 individual of the population, the basic crossover operation applies binomial crossover to create the test vector .
Among them, denotes a random number uniformly distributed between . denotes a random integer in the range , which guarantees that at least one gene is inherited from the parent individual. is the crossover probability, usually in the interval .
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.
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. wherein and are child nodes generated by two parent nodes and . follows the polynomial probability distribution and is calculated as follows.
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.
Considering 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 . After the algorithm starts to search, before the iteration proceeds to , 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.
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 and the diversity criterion . Assume that the corresponding set of the solution set obtained by the multiobjective evolution algorithm in the target space is , and the theoretically optimal Pareto front is .
Definition 4 (convergence criterion ). The convergence criterion is used to measure the degree to which approaches . In this paper, , and . , , and , where is the number of Pareto front obtained by the multiobjective evolutionary algorithm, also known as the number of Pareto set. , , and .
Obviously, the smaller the value of , the closer is to . In particular, when , .
Definition 5 (diversity criterion ). The diversity criterion is used to measure the degree of diversity distribution of on . where is the Euclidean distance between two consecutive solutions in the objective space, is the arithmetic mean of all distances , and and denote the Euclidean distance between the theoretically optimal extremal solution and the boundary solutions and of the obtained nondominated set. Obviously, the smaller the value of , the better the diversity distribution of . When is perfectly uniformly distributed and contains the extreme target vectors in , .
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 , where both and 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 shows that the identification process converges to the actual wind with a small error. Empirically, if the obtained identification results are within 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 of and for three groups of simulation tests are calculated according to Equations (24) and (25) as shown in Table 2.
Table 2 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. where the units of and are also m/s.
According to Section 2.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 . The tests are taken for three time periods for simulation, namely, (), (), and (), 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 and 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.
The comparison between Tables 3 and 4 shows that the algorithm still has a good identification accuracy () 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 , 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.
Data Availability
All data are from simulation experiments.
Conflicts of Interest
There are no conflicts of interest to declare.
Acknowledgments
The authors acknowledge the financial support provided by the National Natural Science Foundation of China (Project No. 11472136).