Abstract
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.
1. 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–4], which usually appears in the indirect methods [5–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–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–17].
Based on the type of the objective function, the optimization problem can be classified by a dynamic or a static optimization problem [18–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–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 proportional-integral-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.
2. 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 is the velocity, is the flight path angle, is the heading angle, is the radial distance from the Earth center to the vehicle, is the longitude, is the geography latitude, is vehicle’s mass, is the thrust, is the specific impulse of the propulsion, is the AOA, is the bank angle, and is the gravity acceleration. where and represent components of the earth’s gravitational acceleration; lies in the direction of the geocentric vector, while is in the direction of the earth’s rotational angular velocity. represents the angular velocity of the earth’s rotation.
An orbit in the three-dimensional space can be represented by six orbital elements [16]: where represents the orbital inclination, the eccentricity, the orbital angular momentum, the right ascension of ascending node (RAAN), and the argument of perigee; and determine the shape of the orbit, and , and determine the position of the orbit in the inertial space.
For an elliptical orbit, the following formula can be obtained: where is the semimajor axis of the orbit and 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 as the eccentricity vector and as the angular momentum vector; then, and can be represented in the inertial coordinate system and calculated by the following way: where denotes the position vector of the vehicle and the velocity vector.
Define and . Then, the RAAN and the argument of perigee can be calculated as follows: where denotes the first component of and the second; is the third component of .
The trajectory optimization problem can be represented as follows: where denotes the states of the system, for example, the trajectory states of the launch vehicle; the initial states; the optimal or feasible law of control; the differential equations of the system, for example, the equations of motion of the launch vehicle; the constraints that should be satisfied when searching for , including equality and inequality constraints; and the performance index.
The formulation Equation (12) is detailed below.
2.2. 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:
Equation (13) can also be represented as follows:
The bank angle is usually neglected in the ascent trajectory optimization problem for a launch vehicle, so . 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: (b)Path constraints
Constraints derived from the structural strength can be represented as follows: where represents the normal overload, the axial overload, the shear force, and the bending moment; , , , and represent the allowed maximum values of , , , and , respectively.
The normal and axial overloads can be calculated as follows: where 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 units with () 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 (denoted as ) is obtained as follows: where the superscript () denotes the part between section 0 to section ().
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 is the allowed maximum AOA and 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:
3. 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. 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–24) and minimizes the (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 denotes the initial pitch angle of the third stage, the pitch angle at the time when, and the pitch angle at the time when ; , , and are unknown coefficients. is the working time of the propulsion of the third stage; and are intermediate variables.
Since the initial states of the third stage are given, however, the value of is known. In addition, the values of , , and are difficult to calculate because the numerical boundaries are generally not known in prior. Therefore, we chose , , and (where denotes the pitch angle at the time when ) 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 that satisfies all the constraints, the iteration of is detailed as follows: (1)Set the initial value of (2)Based on the equations of motion, the overall parameters, for example, the magnitude of the thrust and the mass of the vehicle, the initial states , and the value of , calculate the trajectory in the third stage of the vehicle via numerical integration(3)Obtain the satisfaction of the terminal constraints and the path constraints(4)Obtain the performance index(5)According to the satisfaction of the constraints and value of the performance index, revise using the optimization algorithm(6)If the optimal solution is generated based on the convergence criterion of the optimization method, stop iteration and output , , and the trajectory; otherwise, go to step 2
The optimization algorithm is proposed in Section 4.
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 th particle has ever found is represented by. The superscript “” denotes the generation of the evolution, denotes the inertia weight that indicates how much the particle is affected by its own velocity, and denotes the learning factor; denotes the influence of individual experience on the th particle, and the influence of group experience on the th particle. denotes uniformly distributed random numbers generated in each generation of evolution.
Let and , 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.
4.2. The Improved PSO
According to Equation (27), the evolution of PSO is directly influenced by two random numbers ( and ). 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 the searching capacity and the convergence accuracy, which mainly include evolution formula improvement, parameter selection, variation strategy, and machine learning [30–33].
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 denotes the Gaussian perturbation term, which is added to the th dimension of the th particle’s position in the th generation; and are another two random numbers, and is a random factor subject to the average value (denoted as ) and the variance (denoted as ).
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:. 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 , the searching space of the particle is quickly reduced () because effects of the random values ( and ) do not exist. So, the particle can only search in the vicinity of 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. (a), , and
In this case, Equations (33) and (34) are rewritten as follows:
In PSO, the th particle moves to the global optimal position. But if is the local optimal position, the swarm will quickly gather there and finally obtains a suboptimal solution. But in the improved PSO, the th particle obtains a chance of jumping out of the local optimal position. (b) and
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 denotes the in the th generation and is the threshold.
5. 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 s and 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: °, °/s, , , N, and N·m. In addition, the terminal constraints are given in the form of target-orbital elements: , km, °,°, °, and °. 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 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 optimal solution derived from the improved PSO is °, °, and °; the corresponding pitch angle can be seen in Figure 8.
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 moment when the flight time of the third stage equals to 45 s; for example, the distributions of the axial load, the shear force, and the bending moment are, respectively, shown in Figures 10 and 12.


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.
6. 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.
Acknowledgments
This paper is supported by the China Postdoctoral Science Foundation (2019M661290) and Postdoctoral Science Foundation of Heilongjiang Province of China (LBH-Z19060).