Multiconstrained Ascent Trajectory Optimization Using an Improved Particle Swarm Optimization Method

This paper researches the ascent trajectory optimization problem in view of multiple constraints that effect on the launch vehicle. First, a series of common constraints that effect on the ascent trajectory are formulated for the trajectory optimization problem. Then, in order to reduce the computational burden on the optimal solution, the restrictions on the angular momentum and the eccentricity of the target orbit are converted into constraints on the terminal altitude, velocity, and flight path angle. In this way, the requirement on accurate orbit insertion can be easily realized by solving a three-parameter optimization problem. Next, an improved particle swarm optimization algorithm is developed based on the Gaussian perturbation method to generate the optimal trajectory. Finally, the algorithm is verified by numerical simulation.


Introduction
As one of the most multifunctional vehicles, the launch vehicle has been extensively applied to carry out spacecrafts since 1960s. It can effectively support the requirements on satellite launching, for example, meteorological monitoring, deep space exploration, and manned space flight. Nowadays, technologies about launch vehicle reflect the highest level of science and technology of a country.
In order to send the upper stage payload into the predetermined orbit accurately, the launch vehicle should take off in accordance with the specified launch direction at the specified time and finally reach the specified altitude, velocity, and flight path angle at the time of shutdown. In addition, a series of path constraints should be taken into consideration during the ascent phase, so the requirements on flight security and the performance on the control system can be satisfied. Finally, in order to maximize the performance of the ascent trajectory in a certain launching mission, the specified objective function should be incorporated into the design of the trajectory. As a consequence, an optimization problem is formulated for the ascent trajectory of the launch vehicle.
Trajectory optimization of the launch vehicle is in general considered a typical nonlinear optimal control problem [1]. Initially, analytical methods are developed to solve ascent trajectory optimization problems based on the optimal control theory [2][3][4], which usually appears in the indirect methods [5][6][7]. In this kind of method, the optimization problem is first converted by formulating a Hamilton function and then solved by solving a two-point boundary value problem using the maximum principle [8]. But in many multiconstrained optimization problems that are of great nonlinearity, analytical solutions are not easy to obtain, limiting the application of the maximum principle. Numerical methods suit better for multiconstrained problems and are extensively applied to the design of ascent trajectory. In numerical methods [9], the original problem is converted to a nonlinear system with discretized constraints, states, and optimization parameters. This kind of conversion can be realized by many developed methods, for example, the shooting method, the pseudospectral method, the convex method, and heuristic algorithms [10][11][12][13]. Then, the nonlinear planning problem is solved by a group of classical optimization methods, including the sequential quadratic programming, the simulated annealing, particle swarm optimization (PSO) method, and the genetic algorithm [14][15][16][17].
Based on the type of the objective function, the optimization problem can be classified by a dynamic or a static optimization problem [18][19][20]. In general, a trajectory optimization problem is considered a dynamic problem because the states, the constraints, and the objective function are considered to be dynamical. However, by converting the original problem into a static problem, the calculation of the optimal ascent trajectory can be facilitated because the solution is determined by a series of independent characteristic parameters; then, the nonlinearity of the problem is largely relieved and can be tackled by a developed optimization algorithm.
Generally, heuristic algorithms have fewer restrictions on the objective function and can be applied to more optimization problem in comparison with direct and indirect algorithms. In addition, without any requirements on the gradient information, heuristic algorithm could search the searching space better and easily find the global optimal solution. Thus, heuristic algorithms have been extensively applied to different kinds of optimization problems and have been obtaining rapid developments from researchers. However, the computational efficiency of heuristic algorithms is generally affected by two problems: premature and local optimum. Since proposed by Eberhart and Kennedy in 1995 [21], many improvement methods have been proposed to address the shortcomings of PSO, for example, adaptive PSO, local searching PSO, quantum PSO, and immune PSO [22][23][24][25].
A number of academic contributions have been published since the PSO algorithm has been well studied. Becerra [26] proposed the PSOPT method to achieve sparse nonlinear programming and automatic differentiation, and it incorporates automatic scaling and mesh refinement facilities. To avoid the calculations needed in the common analytical approaches, a new method is proposed using a PSO method [27] for solving an optimal control problem. Chaotic PSO [28] is used to solve the classification problem. An optimization method based on the PSO and LQR technique [29] was proposed to tune the controller gains of digital proportionalintegral-derivative (PID) parameters for a DC motor. This paper researches the optimal ascent trajectory for the launch vehicle. A series of constraints are formulated for the trajectory optimization problem, including various orbital elements, overload, shear force, and bending moment. Constraints, especially equity constraints, can increase the difficulties on solving an optimization problem; so, in this paper, the terminal constraints are first converted into forms that can be easily addressed. In this way, the multiconstrained, highly nonlinear optimization problem is converted to a parameter optimization problem that is easier to solve. Next, an improved PSO method is developed by proposing a Gaussian perturbation term. Theoretical analysis indicates that the revised PSO has better searching capacity in comparison with the basic PSO and numerical simulation results demonstrate the efficiency of the algorithm.
The rest of the paper is organized as follows. Section 2 details the mathematical models of the ascent trajectory optimization problem. Section 3 formulates the approach that a feasible trajectory is designed. Section 4 proposes an improved PSO optimization method to solve the problem. Simulation is carried out in Section 5, and conclusions are drawn in Section 6.

Formulation of the Ascent Trajectory
Optimization Problem 2.1. Mathematical Models of the Launch Vehicle. The equations of motion of the launch vehicle are as follows [5]: where V is the velocity, γ is the flight path angle, ψ is the heading angle, r is the radial distance from the Earth center to the vehicle, θ is the longitude, ϕ is the geography latitude, m is vehicle's mass, P is the thrust, I sp is the specific impulse of the propulsion, α is the AOA, σ is the bank angle, and g is the gravity acceleration.
g v = g r sin γ + g ω sin γ sin ϕ + cos γ cos ϕ cos ψ ð Þ , ω v = ω 2 e r cos ϕ sin γ cos ϕ − cos γ sin ϕ cos ψ ð Þ , ω γ = 2ω e cos ϕ sin ψ + ω 2 e r cos ϕ V cos γ cos ϕ + sin γ sin ϕ cos ψ ð Þ , ω ψ = 2ω e cos ϕ cos ψ tan γ − sin ϕ ð Þ + ω 2 e r V cos γ sin ψ sin ϕ cos ϕ, where g r and g ω represent components of the earth's gravitational acceleration; g r lies in the direction of the geocentric vector, while g ω is in the direction of the earth's rotational angular velocity. ω e represents the angular velocity of the earth's rotation. An orbit in the three-dimensional space can be represented by six orbital elements [16]: For an elliptical orbit, the following formula can be obtained: where a o is the semimajor axis of the orbit and μ E is the gravitational parameter of Earth. The angular momentum can be calculated with the flight states as follows: According to the law of conservation of energy, The semimajor axis can be derived as follows: Based on Equation (4), the eccentricity is calculated as follows: The orbit inclination is calculated as follows: Represent e o as the eccentricity vector and h o as the angular momentum vector; then, e o andh o can be represented in the inertial coordinate system and calculated by the following way: where r denotes the position vector of the vehicle and v the velocity vector.
Define n 1 = ½0, 0, 1 T and n 2 = n 1 × h o . Then, the RAAN and the argument of perigee can be calculated as follows: where n 2x denotes the first component of n 2 and n 2y the second; e oz is the third component of e o . The trajectory optimization problem can be represented as follows: where x denotes the states of the system, for example, the trajectory states of the launch vehicle; x 0 the initial states; u the optimal or feasible law of control; f ðxÞ the differential equations of the system, for example, the equations of motion of the launch vehicle; C u the constraints that should be satisfied when searching for u, including equality and inequality constraints; and J the performance index. The formulation Equation (12) is detailed below.

Trajectory Constraints.
In the ascending phase of a launch vehicle, a series of constraints, which can be generally classified by terminal constraints and path constraints, should be carefully considered. In this paper, the terminal constraints are proposed to satisfy the requirement on accurate target-orbit insertion; at the same time, the path constraints are proposed to guarantee the flight security of the launch vehicle. The trajectory constraints are formulated in the following.

(a) Terminal constraints
At the end of the ascent phase, the trajectory states should meet the specified orbital elements in order to achieve accurate orbit insertion:  (13) can also be represented as follows: The bank angle is usually neglected in the ascent trajectory optimization problem for a launch vehicle, so σ = 0. When the launch position is given, the orbital inclination can be satisfied by determining the launching azimuth, and the RAAN can be satisfied by setting the launching time. Therefore, the number of the terminal constrains is reduced to four: Constraints derived from the structural strength can be represented as follows: where N m represents the normal overload, N x the axial overload, Q y the shear force, and M y the bending moment; N m,max , N x,max , Q y,max , and M y,max represent the allowed maximum values of N m , N x , Q y , and M y , respectively. The normal and axial overloads can be calculated as follows: where D b denotes the axial aerodynamics force. In order to find the distribution of the axial load along the axial direction of the launch vehicle, it is necessary to establish a discrete model of load calculation. From the top to the bottom, the vehicle is divided into several units along the axial direction; therefore, a series of discrete points, which are connected successively by elastic units without mass, are obtained to accomplish load analysis.
Consider the vehicle is divided into s units with (s + 1) discrete points, as shown in Figure 1.
Consider that the inertia force equals to the external force derived from the structure; the longitudinal load of the section i (denoted as N x,i ) is obtained as follows: where the superscript (0 − i) denotes the part between section 0 to section i (i = 0, 1, ⋯, s).
If the mass of the vehicle is uniformly distributed along the axial direction, then Therefore, In order to facilitate the performance of the control system, the angle of attack (AOA) and its rate of change should be constrained as follows: where α max is the allowed maximum AOA and _ α max is the allowed maximum change rate of AOA.
Obviously, there are four terminal constraints and six path constraints in the optimization problem.
Finally, in order to minimize the requirement on the lateral maneuver acceleration, the performance index is selected as follows:

Design of the Optimal Ascent Trajectory
According to the equations of motion of the launch vehicle, the trajectory is mainly determined by the AOA and the bank angle if we take the magnitude of the thrust as constant and the aerodynamic forces as the function of the AOA and the velocity. In general, the value of the bank angle is very small (equals to zero mostly) when we design the optimal ascent trajectory of the launch vehicle, so the trajectory is determined only by the AOA. In addition, the pitch angle is generally used to design the trajectory of the launch vehicle.

International Journal of Aerospace Engineering
The pitch angle is calculated as follows: where φ denoted the pitch angle. Thus, the trajectory optimization problem requires us to find the time history of φ that satisfies Equations (18)(19)(20)(21)(22)(23)(24) and minimizes the J opt (see Equation (22)). This paper takes the third stage of the launch vehicle for example and studies the trajectory optimization algorithm. During the third stage in the ascent phase, the aerodynamic forces have such a small impact on the trajectory compared to the thrust and the gravity that can be neglected in this optimization problem.
Here, the pitch angle during the third stage is designed as a three-segment piecewise linear function as follows: where φ 0 denotes the initial pitch angle of the third stage, φ 1 the pitch angle at the time whent = t 31 , and φ 2 the pitch angle at the time when t = t 32 ; k 1 , k 2 , and k 3 are unknown coefficients. t 3 is the working time of the propulsion of the third stage; t 31 and t 32 are intermediate variables.
Since the initial states of the third stage are given, however, the value of φ 0 is known. In addition, the values of k 1 , k 2 , and k 3 are difficult to calculate because the numerical boundaries are generally not known in prior. Therefore, we chose φ 1 , φ 2 , andφ 3 (where φ 3 denotes the pitch angle at the time when t = t 3 ) as the optimization parameters: Then, unknown coefficients can be calculated as follows: The time history of the pitch angle is illustrated in Figure 2.
In order to find the optimal u that satisfies all the constraints, the iteration of u is detailed as follows: (1) Set the initial value of u The optimization algorithm is proposed in Section 4.

The Improved Particle Swarm Optimization Algorithm
4.1. The Basic PSO. This paper develops an improved PSO algorithm to calculate the optimal ascent trajectory. In PSO, a particle revises its position and velocity vectors by tracking its individual optimal solution as well as the global optimal solution of the swarm. This kind of updating gives the global searching capacity of PSO in the early stage of evolution and the local convergence capacity at the end of evolution. The evolution of the swarm is as follows: where the optimal position that the i th particle has ever found is The superscript "G" denotes the generation of the evolution, 0 < w pso < 1 denotes the inertia weight that indicates how much the particle is affected by its own velocity, and c 1 , c 2 ∈ ½0, 4 denotes the learning factor; c 1 denotes the influence of individual experience on the i th particle, and c 2 the influence of group experience on the i th particle. r 1 , r 2 ∈ ½0, 1 denotes uniformly distributed random numbers generated in each generation of evolution. Let ζ 1 = c 1 r 1 + c 2 r 2 and ζ 2 = c 1 r 1 p i,j + c 2 r 2 p g,j , then we have If the movement of a particle is regarded as a continuous process, Equation (29) can be treated as a classical inhomogeneous second-order differential equation without velocity term. Equation (29) also indicates that the motion of a particle does not need to be described using the velocity, so the process of evolution can be simplified to improve the searching efficiency. (27), the evolution of PSO is directly influenced by two random numbers (r 1 and r 2 ). This kind of randomness can increase the probability of finding a better position but can also affect the convergence and the optimality of the solution. Many scholars have taken extensive efforts on the PSO method to balance   International Journal of Aerospace Engineering the searching capacity and the convergence accuracy, which mainly include evolution formula improvement, parameter selection, variation strategy, and machine learning [30][31][32][33].

The Improved PSO. According to Equation
In this paper, the basic PSO is modified by introducing the Gaussian perturbation into the evolution of the particle's position. The improved PSO is as follows: where where η G i,j denotes the Gaussian perturbation term, which is added to the j th dimension of the i th particle's position in the G th generation; r 3 and r 4 are another two random numbers, and R is a random factor subject to the average value (denoted as μ R ) and the variance (denoted as σ 2 R ). In order to demonstrate the improvements of PSO on the computation efficiency, the basic PSO and the proposed improved PSO are compared in view of four common cases as follows: The change in the velocity of the basic PSO and that of the improved PSO are as follows: Compare Equation (33) and Equation (34)), there is a deviation item:c 2 r 2 r 3 η G i,j . Therefore, by proposing a Gaussian perturbation term in PSO, the particle has a larger searching space; so, the global searching capability is facilitated.
In this case, Equations (33) and (34) are rewritten as follows: In the basic PSO, without η G i,j , the searching space of the particle is quickly reduced (w pso < 1) because effects of the random values (r 1 and r 2 ) do not exist. So, the particle can only search in the vicinity of x G i,j and can be easily trapped in the local optimum. By using a Gaussian perturbation term, local optimum can be better avoided and thus, larger opportunity of finding the best position is obtained.
In this case, Equations (33) and (34) are rewritten as follows: In PSO, the i th particle moves to the global optimal position. But if p i,j is the local optimal position, the swarm will quickly gather there and finally obtains a suboptimal solution. But in the improved PSO, the i th particle obtains a chance of jumping out of the local optimal position.
In this case, Equations (33) and (34) are rewritten as follows: It can be seen that the Gaussian perturbation term affords the particle a larger searching space, which facilitates global searching.
Finally, the PSO algorithm terminates when the following condition is satisfied: where J ðGÞ pso denotes the g best in the G th generation and ε J is the threshold.

Numerical Simulation
The initial trajectory states of the launch vehicle are listed in Table 1. The engine thrust is 26.1 kN, the specified impulse is 240 s, the engine working time is 48.0 s, and the initial mass of the vehicle is 1200 kg. In order to calculate the shear force and the bending moment, set the length of the vehicle as 2.45 m. In addition, set t 31 = 10 s and t 32 = 40 s; so the time duration of each subphase in the third phase (see Figure 2) is 10 s, 30 s, and 8 s, respectively; set the initial pitch angle as -43.4°.
Then, the path constraints of the launch vehicle are set as follows: α max = 40°, _ α max = 5°/s, N m,max = 10:0, N x,max = 3, Q y,max = 400 N, and M y,max = 200 N·m. In addition, the terminal constraints are given in the form of target-orbital elements: e * of = 0:001532, a * of = 6770 km, Ω * of = 30°,i * of = 60°, ω * of = 50°, and θ * of = 40°. The launching position is located at (E100°, N40°); so, the initial heading angle equals to 40.75°according to Equation (9). After conversion, the terminal constraints are represented by using the altitude, the velocity, and the flight path angle which are also listed in Table 1.
Search for the u = ½φ 1 , φ 2 , φ 3 T by using the improved PSO, the results of the optimization problem are illustrated in Figures 3-10. In order to compare with the basic PSO algorithm, we solve the problem using the same initial conditions.
The simulation results indicate that all the path constraints and terminal constraints are well satisfied: in view of the absolute value, the maximum AOA is 39.57°, the maximum rate of change of the AOA is 0.286°/s, the maximum normal load is 2.61, the maximum axial load is 2.03, the maximum shear force is 391.6 N, and the maximum bending moment is 188.6 N·m.
The axial load, the shear force, and the bending moment vary along the body of the launch vehicle, so only the maximum values are given in the simulation results, as illustrated in Figures 9 and 11. In addition, take the values at the To compare with basic PSO, the results in Table 2.
According to the numerical results of the basic PSO, in view of the absolute value, the maximum AOA is 46.28°, the maximum rate of change of the AOA is 0.395°/s, the maximum normal load is 2.77, the maximum axial load is 2.05, the maximum shear force is 459.8 N, and the maximum bending moment is 224.6 N·m. The simulation results indicate that basic PSO algorithm fails to achieve the path constraints and terminal constraints.

Conclusions
This paper researches the optimal ascent trajectory for the launch vehicle. A series of constraints are formulated for the trajectory optimization problem, which dramatically increase the nonlinearity of the optimization problem and makes the solution more difficult. In order to facilitate the computation and the convergence of the optimization method, the multiconstrained, nonlinear optimization problem is converted into a parameter optimization problem that is easier to solve. Next, an improved PSO method is developed by proposing a Gaussian perturbation term. By comparing the improved PSO and the basic PSO in consideration of four common cases, the improvement in the searching capacity is theoretically proved.
Finally, numerical simulation demonstrates the efficiency of the algorithm: in view of over ten constraints, the optimal, feasible ascent trajectory is obtained with the performance index minimized.

Data Availability
The data used to support the findings of this study were related to trade secrets. Requests for data will be considered by the corresponding author in the future if necessary.

Conflicts of Interest
There is no conflict of interest regarding the publication of this paper.