ZHENG Xu,GAO Yejun,JING Wuxing,and WANG Yongsheng

1.Department of Aerospace Engineering,Harbin Institute of Technology,Harbin 150001,China;2.The 28th Research Institute of China Electronics Technology Group Corporation,Nanjing 210007,China;3.Beijing Aerospace Long-march Flying Vehicle Institute,Beijing 100076,China;4.Beijing Aerospace Technology Institute,Beijing 100074,China

Abstract:In the case of the given design variables and constraint functions,this paper is concerned with the rapid overall parameters design of trajectory, propulsion and aerodynamics for long-range ballistic missiles based on the index of the minimum take-off mass.In contrast to the traditional subsystem independent design, this paper adopts the research idea of the combination of the subsystem independent design and the multisystem integration design.Firstly, the trajectory, propulsion and aerodynamics of the subsystem are separately designed by the engineering design, including the design of the minimum energy trajectory, the computation of propulsion system parameters,and the calculation of aerodynamic coefficient and dynamic derivative of the missile by employing the software of missile DATCOM. Then, the uniform design method is used to simplify the constraint conditions and the design variables through the integration design, and the accurate design of the optimized variables would be accomplished by adopting the uniform particle swarm optimization (PSO) algorithm. Finally, the automation design software is written for the three-stage solid ballistic missile. The take-off mass of 29 850 kg is derived by the subsystem independent design,and 20 constraints are reduced by employing the uniform design on the basis of 29 design variables and 32 constraints, and the take-off mass is dropped by 1 850 kg by applying the combination of the uniform design and PSO. The simulation results demonstrate the effectiveness and feasibility of the proposed hybrid optimization technique.

Keywords: ballistic missile, independent design, multidisciplinary integrated design, uniform design, particle swarm optimization(PSO).

1.Introduction

Under the given design variables and constraint functions,the overall design parameters of propulsion,trajectory and aerodynamics for the long-range ballistic missile need to be rapidly derived to ensure that the take-off mass reaches the minimum so as to satisfy the requirements of overall optimization for each subsystem.The design process of the long-range ballistic missile is composed of the concept design,preliminary design and detailed design[1,2].The research contents in this paper belong to the concept design.Tasks in the concept design need to consider the overall scheme selection, subsystems scheme design and overall parameters optimization, which determines 75% expense in the whole design process[3,4].

In the traditional design, the independent management of each subsystem sequentially completes the tasks and passes design results to the next procedure after giving the requirements of the overall system [5,6]. The overall design scheme satisfying tactical index requirements would be acquired by the repeated iteration.The traditional independent design of subsystems captures the essence factors that design variables impact on the design index and simplifies the design process.However,this design ideas have several drawbacks [7,8]. For example, subsystem design results ignore the coupling among each subsystem due to the intentional coordination,which conversely affects the feasibility and optimality of the overall scheme.

To overcome the defects of independent design of the subsystem, multidisciplinary design optimization (MDO)[9–11] such as considering trajectory and propulsion is brought to the forefront with the development of computer and optimization, whose design flow is shown in Fig. 1.The arrows in Fig.1 show the internal relationships during the integrated design process. By calculating engine parameters, mass parameters and aerodynamic parameters,the trajectory is optimized and overall parameters satisfying the constraints are obtained through continuous iteration. The integrated design is able to deal with mutual coupling among subsystems, get over relative independence among disciplines,shorten design cycle,reduce development funds, and improve the performance of the overall system[12].Using a diversity of the intelligent optimization algorithm, some scholars get successful works with respect to multidisciplinary integrated design problems of flight vehicles. In [13,14], the integrated design was conducted through the combination of the propulsion system, guidance control system, trajectory, and aerodynamic shape; 24 design variables of multidiscipline were optimized and a satisfactory result was obtained that the range of the ballistic missile increased 11% more than that of sole trajectory optimization. Zeeshan et al. [15]proposed an MDO strategy for the conceptual design of a multistage interceptor by applying a hybrid optimization method composed of the genetic algorithm(GA)and sequential quadratic programming(SQP),which provides designers with an original and impressive approach.Zafar[16,17] developed a hybrid algorithm with a switch to a local search algorithm for multidisciplinary design of the solid launch vehicle at the conceptual level.In[18,19],an effort was presented to optimize the multidisciplinary design of an entire space launch vehicle to the low earth orbit,consisting of multiple stages by using the GA with the goal of minimizing vehicle weight.Sun et al.[20]set up a mutual coupling optimization model involving the propulsion system of multistage solid rocket, exterior trajectory and aerodynamic characteristics, and optimized 15 variables by employing a variety of the intelligent optimization algorithm.

Fig.1 Flow chart of integrated design

It can be seen that the integrated design of multisystem displays favorable results when coping with multidisciplinary coupling problems. Nevertheless, the integrated design itself does not make full use of already existed independent design criteria, thus the initial values have a large deviation from independent design parameters, which easily plunges into the local optimum and largely differs from practical engineering [21]. To overcome this deficiency,a set of superior initial values in this paper will be obtained by the independent design to conduct the further developed integrated design algorithm.

Therefore,combining the advantages of the independent design of subsystems and the integrated design of multisystem, an alliance of subsystem independent design and multidisciplinary integrated design for the long-range ballistic missile along with applying the uniform design and particle swarm optimization(PSO)is conducted.The elaborate selection for design variables,constraints,initial values and initial population will remarkably improve the efficiency of the PSO. The organization of this paper is shown as follows: firstly, overall preliminary design parameters are separately acquired by the application of the trajectory model, propulsion system model and aerodynamic model; then, design variables of each subsystem would be optimized by adopting the combination of the uniform design and PSO based on the index of the minimum launch mass; finally, simulation results will be derived to validate the effectiveness of the proposed optimization method.

2.Independent design of trajectory,propulsion,and aerodynamic shape

Once warhead mass, launch point and target point are determined, the main design parameters of the overall scheme for ballistic missiles can be preliminarily obtained by the engineering design method.

2.1 Minimum energy trajectory design

The minimum energy trajectory(MET)is the flight trajectory where the shutdown speed of the ballistic missile is minimal.Under this situation,the carried fuel is the least,which guarantees the minimum take-off mass. The initial design of MET is completed within the model of the twobody problem,whose flight path is shown in Fig.2.

Fig.2 Schematic diagram of MET

The shape of the elliptical trajectory can be determined after knowing semi-major axis aband eccentricity e of ellipse.According to the two-body problem,the speed of the shutdown time is

where μ is a geocentric gravitational constant,which is the radius of the earth at the launch point.

The distance between point M in the trajectory and the center of the earth is

where E is eccentric anomaly of point M.

Therefore,when a missile is located at the launch point,we have

Thus, the semi-major axis abof the ellipse can be expressed as

The relationship between the eccentric anomaly E and true anomaly θ is,

The true anomaly of the launch point can be derived from the range angle,that is

Substituting (4), (5) and (6) into (1), we can get that the shutdown velocity is the function of the eccentricity e.Thus, VKwould be the minimum when the derivative of VKto e equals zero. In this case, the eccentricity of the MET is

The shutdown velocity VKof the MET can be acquired by substituting eminVKand other relative parameters into(1).

As a result of ignoring some factors such as rotation of the earth and atmospheric drag, the derived shutdown velocity has a deviation from the accurate value. Therefore,the accurate shutdown velocity can be determined through the iteration method[22].

2.2 Design of propulsion system

2.2.1 Calculation of propellant mass of each stage

The total engine velocity VIcorrelates with VK,that is,

where ΔVgiis the gravity loss velocity in each stage,ΔVD1is the drag loss velocity in the first stage, ΔVT1is the nozzle pressure loss velocity in the first stage,V0is the velocity component in the trajectory plane caused by the rotation of the earth.

VIcan be expressed as

where Isiis the vacuum specific impulse of the engine in the ith stage, μkiis the propellant mass ratio in the ith stage, namely, the ratio of propellant mass mpiin the ith sub-stage to initial mass m0iin the ith stage.

The take-off mass of the ballistic missile[23]is

where mmis the mass of warhead,μFiis the engine mass ratio in the ith stage,namely the ratio of engine mass mFiin the ith sub-stage to initial mass m0iin the ith stage,Niis the structure coefficient in the ith stage,namely the ratio of mass of tail and shell in the ith stage to the initial mass m0iin the ith stage.The relationship of these parameters in each stage is shown in Fig.3.

Fig.3 Schematic diagram of N stages of ballistic missile

In the case of the known warhead mass, the following function needs to be maximized in order to minimize the take-off mass of the ballistic missile[23],namely

The maximal value of function U with the constraint of(9) can be obtained by applying the Lagrange multiplier method,whose expression is shown as

The take-off mass can be determined by the propellant mass ratio in the ith stage in the case of the given engine mass ratio and vacuum specific impulse.The extremum of function F would be worked out by the following function,that is,

Equation(13)contains n+1 unknown quantities of μki(i = 1,...,n)and λ,and has n+1 equations,which can solve the optimal propellant mass ratio μkiof each stage.

The recursive formula of the initial mass of each stage can be expressed as

Thus,the propellant mass of each stage is

2.2.2 Design of engine parameters

The parameters of an engine principally include the engine diameter,the length of the cylindrical section,the average pressure of the combustion chamber,the expansion ratio of the nozzle,the expansion semi-angle of the nozzle and the working-time of the engine.

The preliminary design value of engine diameter Dciand the length of cylindrical section Lciin each stage can refer to parameters of the same type ballistic missile. The design of the average pressure of the combustion chamber needs to ensure normal combustion of the propellant as well as the maximum impulse provided by the unit mass of the engine.In general,the average pressure of the combustion chamber in the first stage is optional for 6–10 MPa,and it would be lower in the second and upper stage with the optional value of 4–5 MPa.

The expansion ratio of the nozzle is computed[24]as

where peis the exit pressure of the nozzle with the value of 0.045–0.075 MPa in the first stage and 0.015–0.03 MPa in the second and upper stage,and k is the ratio of specific heat of gas.

The working-time of the engine can be written as

within which

where C?is the characteristic velocity,CFis the thrust coefficient, ηcis the impulse coefficient of the combustion chamber, ηnis the impulse coefficient of the nozzle, αMis the initial expansion semi-angle of the nozzle,αeis the exit expansion semi-angle of the nozzle, N01is the ratio of thrust to weight in the first stage, g0is the gravity acceleration of the standard sea level,ζ is the compensation coefficient.

Because of the significant influence of the objective function induced by the engine diameter,the length of the cylindrical section, and the exit expansion semi-angle of the nozzle,the parameters in this paper would be selected as the optimization variables of engine parameters, and other parameters are computed according to relative formulas.

2.3 Ascent dynamics and flight program

Throughout this paper the notation is such that a boldface letter means a vector, and the plain form of the same letter indicates the magnitude of this vector. The equations of motion of the ballistic missile in an inertial coordinate system can be expressed as

where r and v ∈R3are the inertial position and velocity vectors,g is the gravitational acceleration,Tvacis the full vacuum thrust magnitude,η > 0 is the engine throttle, T is the current thrust magnitude including effects of throttle modulation and thrust loss due to the back pressure.In this formulation, the total engine thrust is assumed to be aligned with the body longitudinal axis, and is not gimbaled independently.The vectors A and N are the aerodynamic forces in the body longitudinal and normal direction,respectively;Ibis the unit vector defining the vehicle body longitudinal axis;m(t)is the mass of the vehicle at the current time t. The specific impulse of the engine is Isp, and g0represents the gravitational acceleration magnitude on the surface of the earth.For better numerical conditioning,the following nondimensionalization is used.

(i)The distances are normalized by R0,the radius of the earth at the equator.

(iii)The velocities are normalized bythe circular velocity around the earth at R0.

The gravity is modeled by the Newtonian central gravity field.With some abuse of notation,we use the same names hereafter for the dimensionless variables.The dimensionless equations of motion from(21)are then

where the prime indicates differentiation with respect to the dimensionless time. Now A and N are the aerodynamic accelerations in g0in the body longitudinal and normal direction,respectively,and T is the magnitude of the thrust acceleration in g0.

Define the dimensionless wind-relative velocity

where ωeis the earth angular rotation rate vector normalized byand vwis the dimensionless wind velocity.

The magnitudes of the dimensionless aerodynamic and thrust accelerations are given by

The flight process in the ascent phrase generally consists of atmospheric flight phrase and vacuum flight phrase, of which the atmospheric flight phrase comprises the vertical ascending phrase(0?t1),subsonic turning phrase(t1?t2)and trajectory turning phrase(t2?tk1).

The ascent flight program of the pitch angle is designed as

where ? is the trajectory inclination angle, α is the angle of attack(AOA),ωezis the component of axis z of angular velocity of the earth,?cx0iis the pitch angle of the initial moment in the ith stage(i 2),is the gradient of the pitch angle in the ith stage(i 2),tki?1is the shutdown time of engine in the(i ?1)th stage.

The AOA in the subsonic turning phrase can be devised as

where αmis the maximal absolute value of the AOA,a is the turning constant coefficient.

2.4 Aerodynamic shape design of warhead

The shape of warhead of the ballistic missile as shown in Fig. 4 is generally a small blunt cone. In order to reduce the resistance of warhead during flight,the ratio μrbmaxof bluntness radius rbto the bottom diameter Dc3of warhead is generally less than 0.1, and the half cone angle can be chosen from 8?to 15?.

Fig.4 Cone shape of warhead

The expressions of bluntness radius rband the length of warhead lmare

In this paper, we calculate aerodynamic parameters of the ballistic missile by adopting the aerodynamic engineering software named“Missile DATCOM”developed by US air force.The manipulation procedure is listed as follows.Firstly, the datum is generated in the light of the aerodynamic shape parameters of missile and warhead.Then,AOA, Mach number, altitude, reference area, reference length,control instructions and other related data would be written into a “dat”file to be read for the subprogram.Finally,the executable program “misdat.exe”is employed to compute aerodynamic parameters,whose results are saved in a file that is to be read at the time of trajectory calculation.

3.Description of integrated design problems

After the independent design of the propulsion system,trajectory and aerodynamic shape,the multidisciplinary integrated design needs to be implemented for achieving the overall optimal performance of the ballistic missile.Above all, design variables, constraints functions and objective functions of the optimization model are illustrated in the following elaboration.

3.1 Analysis of design variables

In the problem of MDO,the optimized design variables of the propulsion system,trajectory design and aerodynamic shape are selected as follows.

Variables of the propulsion system consist of engine diameter Dci, length Lciof the cylindrical section of each stage, exit expansion semi-angle αeiof the nozzle, average working pressure pciof the combustor,expansion ratio εAiof the nozzle, engine working time tai, initial expansion semi-angle αMiof the nozzle,and the exit expansion semi-angle of the nozzle αei

Variables of aerodynamic shape consist of nose bluntness radius rband length lmof warhead.

Variables of the trajectory design consist of launch azimuth A0, vertical ascending time t1, turning constant coefficient a, maximal absolute value αmof the AOA and gradient of the pitch angle in the ith stage

The above optimization design variables can be expressed as the form of vector,that is,

3.2 Analysis of constraint functions

According to the requirements of the overall design and subsystem design,the design of the ballistic missile needs to be satisfied under the constraints of the propulsion system,trajectory and aerodynamic parameters.

The constraints of the propulsion system include that the engine mass ratio meets the requirement of μFminμFat each stage,the difference between the initial and the exit expansion semi-angle of the nozzle meets the requirement of ΔαΔαmaxat each stage,the ratio of the nozzle exit diameter to the body diameter meets the requirement of μdDμdDmaxat each stage,the slenderness ratio of the fuel meets the requirement of kldminkldkldmax,the take-off ratio of the thrust to weight meets the requirement of N01minN01N01max, the ratio of the nozzle exit pressure to the external pressure meets the requirement of

The constraints of the aerodynamic shape include that the ratio of the bluntness radius to the bottom diameter meets the requirement of μbDμbDmax,the half cone angle meets the requirement of θdminθθdmax,the missile slenderness ratio meets the requirement of kLDmin

The constraints of the trajectory include that the axis overload of the missile meets the requirement of nxnxmax, the normal overload of the missile meets the requirement of nynymax, the absolute value of the AOA in the subsonic turning phrase meets the requirement of αmαmmax,the ascent dynamic pressure meets the requirement of qqmax, the separation height between the first and second stages meets the requirement of H12minH12,the trajectory inclination angle at reentry time meets the requirement of ?emin?eθemax, the longitudinal and lateral deviations of the hit point meet the requirements of ΔLΔLmaxand ΔZΔZmax.

The boundary value of the constraint function is shown in Table 1.The constraints of each subsystem can be uniformly expressed as the form of gi0, such as g1=μFmin? μF0. The key and constraint functions are acquired and illustrated as follows:

Table 1 Boundary value of the constraint function

3.3 Analysis of objective function

In engineering, we generally emphasize on the total cost,take-off mass,range of the ballistic missile or other parameters designed as the objective function[25].In the case of the predetermined parameters of warhead mass and missile range,the take-off mass of the ballistic missile is proportional to the cost of the weapon system. As a consequence,the take-off mass is chosen as the evaluation criteria of the propulsion system,trajectory design and aerodynamic configuration to guarantee its minimal value under the condition of multi-constraints.The objective function,constraints and boundary of variables are depicted as

where Xminand Xmaxare respectively the lower and upper bounds of optimization design variables X.

In the optimization problem abstracted from(33),with the requirements of boundary conditions and constraint vectors, the minimum take-off mass is achieved by optimizing the design variables according to the multidisciplinary parameters.

4.Simplification of design variables and constraint functions

Too many design variables and constraints will increase the calculation workload, then reducing the efficiency of the optimization.It has been known that each design variable has different effects on the objective functions and constraints, and each constraint plays different importance in the solution space.Therefore,screening significant design variables and constraints is helpful to improve the computation efficiency. In the case of many design variables,fewer and effective test plans will be implemented by the uniform design,which can be obtained by the relationship between test factors and optimization indexes. Based on above advantages,the uniform design test is employed to simplify the design variables and constraints.

4.1 Simplification of design variables

In the uniform design test, because of the simple principle and easy realization, the grid point method is used to construct the uniform design table as follows, which can generate a series of numbers with n rows and m columns[26].

(i)For a given set of test parameters n, the positive integer h whose value is less than n would be searched to constitute the column vector h=[h1h2··· hm]Twith the requirement that the greatest common divisor between n and h is one.

(ii) The element of the column i in the uniform design table is

where i and j are respectively i = 1,...,n and j =1,...,m.

If the number of experiment factors is s,it haskinds of experiment schemes.Because the uniformity of the selection scheme affects experiment results, the minimum deviation criterion is generally utilized to select the experiment scheme with better uniformity. The centralized L2 deviation is commonly used to evaluate the uniformity criterion,whose formula can be expressed as

where xijis the transform of element uijof the uniform design table,namely,

After a series of design points are obtained by the uniform design table, the partial derivatives (sensitivity) of each design variable to the objective function can be calculated at each point,namely ?f/?Xi.If the mean value of?f/?Xiis large, while the standard deviation is small, it shows that the design variable has a significant influence on the objective function in the design space. Otherwise, the design variable has little influence on the objective function. ?f/?Xiis considered as a normal distribution variable, and the credibility level of its absolute value greater than parameter b is

where b is a given constant value,μiis the mean value of?f/?Xi,σiis the standard deviation of ?f/?Xi,Df is a simplified representation of ?f/?Xi.The greater value of probability P(|?f/?Xi|>b)is,and the greater impact of design variable on the objective function will be.

For the form of constraint function gi0,the concept of cumulative constraints is introduced to describe the effect of design variables on constraint functions.The definition of the cumulative constraint is

where p is a parameter with a value of p = 1.The partial derivatives ?C/?Xiof the cumulative constraint to each design variable at each test point can be computed by using the centralized L2 deviation method. The credibility level of its absolute value greater than parameter c is

where DC is a simplified representation of ?C/?Xi.

Similarly, the greater value of probability P (|?C/?Xi|>c) is, and the greater impact of the design variable on constraint function will be.

4.2 Simplification of constraint functions

From a statistical point of view, if the average value of a constraint function is greater than zero at all test points and its standard deviation is smaller,it shows that this constraint plays an important role in the solution space. On the contrary, if the mean value of a constraint function is negative and its standard deviation is relatively small,this constraint function is regarded as an unimportant one.Constraint function giis considered as a normal distribution and the probability of gi>0 is shown as follows:

where μiis the mean value at the ith constraint function and σiis the standard deviation at the ith constraint function.

Similarly,the greater value of probability P(gi>0)is,the greater importance the constraint will be of.

5.Integrated optimization of design variables

The PSO, as an intelligent optimization algorithm, can deal with multiple particles simultaneously and carry out the search in parallel, which is widely explored as a result of distinct advantages of simple algorithm and easy implementation[27].Distinct with the GA and simulated annealing, the PSO can eliminate the operation of the crossover and mutation and reveals fast convergence speed and strong search balance capability both globally and locally. However, the PSO displays several shortcomings[28], such as low search accuracy, easily to fall into the local optimum, and large error in solving the optimization problem concerned with discrete variables. In order to overcome these disadvantages of the PSO, this paper pledges the global optimum instead of the local optimum in search by using the uniform design to initialize the particle population[29].

5.1 Initialization of particle population using the uniform design

As a random search algorithm,two search directions will be conducted in the PSO.The first one is to inherit the existing information,e.g.,the optimal particle position of the current population, which guarantees the global convergence of the population.The other one is the random disturbance,e.g.,the existing random weighting coefficients ξ and η in(40),which guarantees the exploration of the population.The initial particle population is randomly generated in most PSO strategies.This initialization is easy and effective in most cases.However,if the random initial particle population cannot fully represent the characteristics in the design space, the global optimum information will be missing, which will lead to the premature of particle population and converge into local optimum.

Therefore, in order to improve the coverage of the initial population in the design space,a series of design points obtained by the uniform design table will be employed to initialize the particles before optimization.Fig.5 shows the comparison between uniform design points and the random point set,which reflects that uniform design whose points are dispersed evenly can better indicate the inherent characteristics of design space and is conducive to containing global optimal information.

Fig.5 Comparison between uniform design and random point set

5.2 Precise design of optimization variables

The left particles in Fig. 5 distribute evenly and orderly,covering the design space nicely and avoiding falling into the local optimum effectively in the iterative optimization process. Each particle in the graph represents each multidisciplinary vector, gradually converging to the optimal point satisfying the constraint index in the process of continuous optimization.

A set of optimization variablesXiin (31)corresponding to each particleiin the particle swarm is a set of parameters in the solution space.Each particle in the particle swarm composed ofnparticles is computed to update its flight velocity and position in each time of search.Thek+1 iteration formula[30]is

whereis the velocity of particleifor variablej,is the position of particleifor variablej,is the best position found by particleiso far for variablej,is the best position in the swarm for variablej,ωis inertia weight to maintain the original velocity,c1is particle weight coefficient to hold the historical optimal position,c2is the particle weight coefficient to track the optimal position of population,ξandηare independent random weighting coefficients between 0 and 1,anddis a constraint factor.

The inertia weightωcontrols the exploration properties of the algorithm,with larger values facilitating a more global behavior and smaller values facilitating a more local behavior.The global convergence performance of the PSO will be improved prominently as the value ofωdecreases with the increase of iteration numbers.Theωin this paper is designed[31]as

whereωmaxandωminare respectively the maximum and minimum of the inertia weight,iteis the current iteration number,anditemaxis the total iteration number.

6.Simulations

6.1 Independent design results of each subsystem

Launch point parameters are set as 116.4?E, 39.3?N,0 m elevation. Target point parameters are set as– 155.5?E, 19.7?N, 0 m elevation.The warhead mass is 900 kg. Through the independent design of each subsystem, the parameters will reach a balance after multiple iterations.The ratios of thrust to weight of each stageN01,N02,N03are respectively derived as 2.5,2.7,3.0,the structure ratios of each stageN1,N2,N3are respectively derived as 0.011, 0.011, 0.012, and the calculated take-off mass of the ballistic missile is 29 850 kg.The propulsion system parameters are displayed in Table 2. The aerodynamic shape parameters of warhead and trajectory parameters are shown in Table 3.The initial variable values obtained by the independent design provide the PSO for superior initial conditions.

Table 2 Designed parameters of propulsion system

Table 3 Designed aerodynamic shape parameters of warhead and trajectory parameters

6.2 Results analysis using uniform design

The following optimization work will be carried out on the basis of the independent design of propulsion,aerodynamic and trajectory parameters.When there are too many design variables and constraints,the calculation workload will be greatly increased, then reducing the optimization efficiency.Therefore,the selection of design variables and constraints is a necessary operation before further optimization.In this paper, the uniform design test technique is used to analyze and simplify the design variables and constraint functions.

There are total 29 design variables in the propulsion system, aerodynamic shape and trajectory for the three-stage ballistic missile.X1?X3respectively represent the engine diameter in each stage,X4?X6respectively represent the length of the cylindrical section in each stage,X7?X9respectively represent the average pressure of the combustion chamber in each stage,X10?X12respectively represent the working-time of each stage engine,X13?X15respectively represent the expansion ratio of the nozzle in each stage,X16?X18respectively represent the initial expansion semi-angle of the nozzle in each stage,X19?X21respectively represent the exit expansion semi-angle of the nozzle,X22?X23respectively represent the bluntness radius and length of warhead,X24?X29respectively represent the launch azimuth,vertical ascending time,turning constant coefficient, maximal absolute value of the AOA,the gradient of the pitch angle in Stage 2 and the gradient of the pitch angle in Stage 3.

There are total 32 constraint functions in the propulsion system,aerodynamic shape and trajectory,which are shown in Table 1.Through the uniform design technique,the constraint function ordering is shown in Table 4, the effect of design variables on objective function is shown in Table 5, and the effect of design variables on cumulative constraints is shown in Table 6.According to the results in these tables, the variables and constraints with significant influence are filtered by probability ranking to reduce the calculation burden of the PSO.

Table 4 Constraint function ordering

Table 5 Effect of the design variables on objective function

Table 4 shows the mean values,standard deviations and probabilities of the constraint functions, and the probabilitiesP(gi >0) are then sorted according to its numerical values. It can be seen that the average values of all constraint functions are less than zero, which satisfies the requirement of the constraint index. The probability ofgi >0 can be computed by the mean value and standard deviation.From Table 4,through sorting and screening the probabilities, the constraintsg2,g7,g12,g19,g20,g21,g26,g27,g28,g30,g31,g32have large probabilities of greater than zero, showing strong constraint capability. However, other constraints have small probability of greater than zero,which will be not considered in the PSO.

Table 5 shows the mean values,standard deviations and probabilities of the derivatives of the objective function to design variables, and the probabilitiesP(|?f/?Xi| > b)are then sorted according to its numerical values.The objective function is the take–off mass, not considering the warhead mass, therefore, trajectory design variables and warhead shape parameters do not have directly effect on the objective function. In order to better distinguish the probabilitiesP(|?f/?Xi| > b), parameterbis set as 14.From Table 5, derivative ?f/?Xipresents significant changes in the mean values and standard deviations. The design variables X1, X2, X3, X8, X9, X10, X16present remarkable influences on the objective function, and the probabilities are relatively large compared with other design variables.Therefore,these design variables are important ones and should be considered for further optimization.

Table 6 Effect of the design variables on cumulative constraints

Table 6 shows the mean values,standard deviations and probabilities of the derivatives of cumulative constraints to design variables.The probabilities P(|?C/?Xi| > c)are sorted according to its numerical values,in which parameter c is set as 1.4.It can be obtained that derivative ?C/?Xialso presents significant changes in the mean values and standard deviations. Different from Table 4 and Table 5,the probability changes in Table 6 are stable,not appearing tiny probabilities.Therefore,the design variables have significant effect on cumulative constraints, and all of them will be taken into account in the process of further optimization.

6.3 Results analysis using uniform particle swarm optimization

After initialization of the particle swarm conducted by the uniform design,the PSO is applied to optimize the design variables,whose flow chart is shown in Fig.6.

Relative parameters are set as follows.Iteration number itermaxis 300. The maximal value ωmaxand the minimal value ωminof inertia weight are respectively 0.9 and 0.1.Learning factors c1and c2are both 2.Constraint factor d is 1 and the number of particles n is 30. Under the above given simulation conditions,optimization values of variables through iteration of 300 generation are shown in Table 7 and the objective function curve during iteration is shown in Fig.7.The constraints satisfied during the optimization process are shown in Table 8. Meanwhile, the AOA, ascent normal overload, ascent dynamic pressure,velocity of the whole trajectory and the height of the whole trajectory are obtained to display results comparisons of pre–optimization and post-optimization from Fig. 8 to Fig.12.

Fig.6 Flow chart of PSO

Table 7 Optimization results of design variables

Fig.7 Variation curve of objective function in the iteration

Table 8 Optimization results of constraint functions

From Table 7, the post-optimization design variables are located between the upper and lower limit, satisfying the constraint boundaries.From Table 8,the selected constraint functions are all less than zero through the optimization,and the post-optimized design variables are within the limits of the constraints. In Fig. 7, the optimized take-off mass of the ballistic missile converges to 28 000 kg reduced by 1 850 kg compared with the pre-optimized value of 29 850 kg after 250 iterations.

Fig.8 Ascent AOA of pre-optimization and post-optimization

Fig. 9 Ascent normal overload of pre-optimization and postoptimization

Fig. 10 Ascent dynamic pressure of pre-optimization and postoptimization

Fig.11 Velocity of pre-optimization and post-optimization

Fig.12 Height of pre-optimization and post-optimization

The pure integrated design [21] does not take the independent design of each subsystem into consideration,which makes the take-off mass reduce too much with the value of 20 t.The parameters of each subsystem basically reach the suboptimal values through the independent design,providing superior initial values for further multidisciplinary integrated design and reducing 2 t for the take-off mass.

In addition, the uniform dispersion of the design variables can commendably present the intrinsic characteristic in the design space through the uniform test design,which makes the particles diverse and is conducive to the initial population for containing the global optimal information.When initializing the particles without using the uniform design technique,in some cases,i.e.,unevenness of space coverage,the optimization index (take-off mass) is larger than the value of 28 000 kg, which indicates the particles fall into the local optimum.Therefore,the uniform test design technique is a very practical initialization strategy to excellently overcome the prematurity of the PSO. The integrated design combined with the propulsion system,trajectory and aerodynamic shape reveals an excellent design efficiency by employing the optimization method of combination of uniform design and the PSO, and the global optimal solution of the take-off mass is acquired under the given payload.As the combination of the independent design and uniform design experiment provides a superior condition for the PSO, the variable optimization results present a better flight performance than baseline cases.

In very few cases, the optimization result index given by the uniform design is better than that given by the nonuniform design.In most cases, the optimization results of the non-uniform design are basically the same as the results of the uniform design.The differences are the iterations and optimization times.In addition,the uniform design density will affect optimization time. The higher the uniform design density, the longer the computation time.However,the optimization accuracy of design variables is basically consistent. Considering the calculation burden,total particles are set as 60, which can satisfy the coverage in the design space.

In terms of the algorithm,the characteristics of different optimization algorithms were illustrated in[28,32].For example,the PSO presents excellent global search ability and the GA presents excellent strong fault tolerance. In most cases,apart from differences in iterations and computation time,the results are basically consistent by using different intelligent optimization algorithms.

In Fig. 8 and Fig. 9, the maximum negative attack angle in the subsonic turning phrase reduces from 14?to 10?through the integrated optimization,and the normal overload value reduces from 0.6 to 0.4.It can be found that the PSO can effectively reduce the constraints and improve the aerodynamic characteristic of the ballistic missile. However, in the flight beginning from Stage 2, the AOA and normal overload through optimization are larger than the pre-optimization.The AOA value at the engine-cutoff increases from 3?to 6.8?,and the normal overload value increases from 0.49 to 0.88.The AOA and normal overload are significantly altered during the vacuum flight phrase to adjust the influence in the subsonic turning phrase.At the end of the first stage,the primary booster is separated,and the constraint range of the normal overload is expanded.Thus, the attitude angle of the missile can be changed to improve the missile’s range and other indicators.

In Fig.10,the dynamic pressure reaches the maximum value at the time of 40 s.The pre-optimized dynamic pressure is 0.12 MPa,exceeding the requirement of the maximum value 0.1 MPa.Through the optimization,the maximum dynamic pressure is 0.097 MPa,which meets the requirement of the index.As the result of decreasing the air density with increasing the flight height,the dynamic pressure decreases gradually and reduces to zero after the time of 90 s basically.

In Fig. 11 and Fig. 12, the whole flight time increases from 1 835 s to 1 920 s through the optimization. The flight height at the engine-cutoff increases from 320 km to 360 km, and the velocity at the engine-cutoff decreases from 6.9 km/s to 6.8 km/s. As the potential energy increases with the improvement of the fight height, the engine-cutoff velocity drops slightly.The height of the ballistic missile reaches the maximum value at the time of around 1 000 s from 1 250 km to 1 550 km. Similarly,the velocity at the peak height reduces from 5.6 km/s to 5.2 km/s as the result of increasing the potential energy.The range of the ballistic missile obviously improves from 8 000 km to 8 500 km. Based on the analysis presented above,the optimal design variables can be acquired under the condition of multi-constraints by employing the proposed integrated optimization technique through the combination of the uniform design and PSO,and the minimum take-off mass is achieved,which demonstrates that the proposed hybrid integrated design is effective and feasible.

7.Conclusions

In the case of the given launch point,target point and warhead mass, this paper studies the overall optimization design of propulsion, trajectory and aerodynamics for the long-range ballistic missile based on the index of the minimum take-off mass.Conclusions are listed as follows.

(i) The independent design of each subsystem is conducted for the propulsion,trajectory,and aerodynamics of the ballistic missile. In the design of the propulsion system,the Lagrange multiplier method is employed to compute the optimal propellant mass of each stage and take-off mass based on the minimum shutdown velocity derived by the MET.In the trajectory design,trajectory parameters are devised to obtain the ascent trajectory of the ballistic missile. In the aerodynamic design, aerodynamic coefficients and dynamic derivatives of the ballistic missile are computed by using the software of missile DATCOM developed by the US air force.Through the independent design of each subsystem,it provides references for the following integrated design.

(ii) The method of combination of the uniform design and PSO is proposed to settle the problems of the integrated design involving 15 variables and 18 constraints from multisystem.In purpose of avoiding premature problems of falling into the local optimal solution, the particle population is initialized by using the uniform design experiment to assure its complete coverage in the design space, followed by further variables optimization adopted by the PSO.

(iii)The take-off mass after further integrated design for the three-stage ballistic missile is 28 000 kg reduced by 1 850 kg compared with the pre-optimized value of 29 850 kg obtained by the independent design of each subsystem.The multidisciplinary integrated design for the long-range ballistic missile distinctly ameliorates ascent aerodynamic properties and augments height and range of the whole trajectory. From the simulation data, the optimized parameters can distinctly improve the flight performances in the integrated design.The proposed multidisciplinary integrated design technique combined with the independent and integrated design is impressive to cope with the overall parameters optimization and exhibits favorable feasibility.