# Multiconstrained Ascent Trajectory Optimization Using an Improved Particle Swarm Optimization Method

Authors: Mu Lin; Zhao-Huanyu Zhang; Hongyu Zhou; Yongtao Shui
Journal: International Journal of Aerospace Engineering
 (2021)

---

## 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.


---

## Body

## 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. Section2 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]:
(1)V̇=Pcosα−Xm+gv+ωv,γ̇=Y+PsinαmVcosσ+Vrcosγ+gγ+ωγ,ψ̇=Y+PsinαmVcosγsinσ+Vrcosγsinψtanϕ+gψ+ωψ,ṙ=Vsinγ,θ̇=Vcosγsinψ,rcosϕϕ̇=Vcosγcosψr,ṁ=−PIspg,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, Isp is the specific impulse of the propulsion, α is the AOA, σ is the bank angle, and g is the gravity acceleration.
(2)gv=grsinγ+gωsinγsinϕ+cosγcosϕcosψ,gγ=grVcosγ+gωVcosγsinϕ−sinγcosϕcosψ,gψ=−gωVcosγcosϕsinψ,ωv=ωe2rcosϕsinγcosϕ−cosγsinϕcosψ,ωγ=2ωecosϕsinψ+ωe2rcosϕVcosγcosϕ+sinγsinϕcosψ,ωψ=2ωecosϕcosψtanγ−sinϕ+ωe2rVcosγsinψsinϕcosϕ,where gr and gω represent components of the earth’s gravitational acceleration; gr 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]:
(3)XO=io,eo,ho,Ωo,ωoT,where io represents the orbital inclination, eo the eccentricity, ho the orbital angular momentum, Ωo the right ascension of ascending node (RAAN), and ωo the argument of perigee; eo and ho determine the shape of the orbit, and Ωo, ωo and io determine the position of the orbit in the inertial space.For an elliptical orbit, the following formula can be obtained:(4)ao=ho2/μE1−eo2,where ao 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:(5)ho=Vrcosγ.According to the law of conservation of energy,(6)V22−μEr=−μE2ao.The semimajor axis can be derived as follows:(7)ao=μEr2μE−rV2.Based on Equation (4), the eccentricity is calculated as follows:
(8)eo=1−ho2μEao.The orbit inclination is calculated as follows:(9)io=arccossinψcosϕ.Representeo as the eccentricity vector and ho as the angular momentum vector; then, eo andho can be represented in the inertial coordinate system and calculated by the following way:
(10)ho=r×v,eo=v×hoμE−rr,where r denotes the position vector of the vehicle and v the velocity vector.Definen1=0,0,1T and n2=n1×ho. Then, the RAAN and the argument of perigee can be calculated as follows:
(11)Ωo=arccosn2xn2,n2y≥0,2π−arccosn2xn2,n2y<0,ωo=arccosn2⋅eon2eo,eoz≥0,2π−arccosn2⋅eon2eo,eoz<0,where n2x denotes the first component of n2 and n2y the second; eoz is the third component of eo.The trajectory optimization problem can be represented as follows:(12)minJs.t.x0,Cu,fx⇒u,where x denotes the states of the system, for example, the trajectory states of the launch vehicle; x0 the initial states; u the optimal or feasible law of control; fx the differential equations of the system, for example, the equations of motion of the launch vehicle; Cu 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.


### 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 constraintsAt the end of the ascent phase, the trajectory states should meet the specified orbital elements in order to achieve accurate orbit insertion:(13)Cequ=eotf−eof∗aotf−aof∗Ωotf−Ωof∗iotf−iof∗ωotf−ωof∗θotf−θof∗.Equation (13) can also be represented as follows:
(14)hf∗=rf∗−ae=ho∗/μE1+eo∗cosθo∗−ae,Vf∗=2μE21+eo∗cosθo∗ho∗−μEao∗,γf∗=arccos1+eo∗cosθo∗eo∗2+2eo∗cosθo∗+1.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:
(15)Cequ=htf−hf∗γtf−γf∗Vtf−Vf∗ωotf−ωof∗.(b)
Path constraintsConstraints derived from the structural strength can be represented as follows:(16)Nm≤Nm,max,Nx≤Nx,max,Qy≤Qy,max,My≤My,max,where Nm represents the normal overload, Nx the axial overload, Qy the shear force, and My the bending moment; Nm,max, Nx,max, Qy,max, and My,max represent the allowed maximum values of Nm, Nx, Qy, and My, respectively.The normal and axial overloads can be calculated as follows:(17)Nm=Psinα+Ymg,Nx=P−Dbmg,where Db 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 intos units with (s+1) discrete points, as shown in Figure 1.Figure 1
Illustration of normal load, axial load, shear force, and bending moment.Consider that the inertia force equals to the external force derived from the structure; the longitudinal load of the sectioni (denoted as Nx,i) is obtained as follows:
(18)Nx,i=P−Db−P−Dbmm0−i=P1−m0−im+m0−im−1Db=P−Db1−m0−im,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(19)m0−i=m1−is.Therefore,(20)Nx,i=isP−Dbi=0,1,⋯,s.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:(21)α≤αmax,α̇≤α̇max,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:(22)minJopt=∫Nmdt.


## 2.1. Mathematical Models of the Launch Vehicle

The equations of motion of the launch vehicle are as follows [5]:
(1)V̇=Pcosα−Xm+gv+ωv,γ̇=Y+PsinαmVcosσ+Vrcosγ+gγ+ωγ,ψ̇=Y+PsinαmVcosγsinσ+Vrcosγsinψtanϕ+gψ+ωψ,ṙ=Vsinγ,θ̇=Vcosγsinψ,rcosϕϕ̇=Vcosγcosψr,ṁ=−PIspg,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, Isp is the specific impulse of the propulsion, α is the AOA, σ is the bank angle, and g is the gravity acceleration.
(2)gv=grsinγ+gωsinγsinϕ+cosγcosϕcosψ,gγ=grVcosγ+gωVcosγsinϕ−sinγcosϕcosψ,gψ=−gωVcosγcosϕsinψ,ωv=ωe2rcosϕsinγcosϕ−cosγsinϕcosψ,ωγ=2ωecosϕsinψ+ωe2rcosϕVcosγcosϕ+sinγsinϕcosψ,ωψ=2ωecosϕcosψtanγ−sinϕ+ωe2rVcosγsinψsinϕcosϕ,where gr and gω represent components of the earth’s gravitational acceleration; gr 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]:
(3)XO=io,eo,ho,Ωo,ωoT,where io represents the orbital inclination, eo the eccentricity, ho the orbital angular momentum, Ωo the right ascension of ascending node (RAAN), and ωo the argument of perigee; eo and ho determine the shape of the orbit, and Ωo, ωo and io determine the position of the orbit in the inertial space.For an elliptical orbit, the following formula can be obtained:(4)ao=ho2/μE1−eo2,where ao 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:(5)ho=Vrcosγ.According to the law of conservation of energy,(6)V22−μEr=−μE2ao.The semimajor axis can be derived as follows:(7)ao=μEr2μE−rV2.Based on Equation (4), the eccentricity is calculated as follows:
(8)eo=1−ho2μEao.The orbit inclination is calculated as follows:(9)io=arccossinψcosϕ.Representeo as the eccentricity vector and ho as the angular momentum vector; then, eo andho can be represented in the inertial coordinate system and calculated by the following way:
(10)ho=r×v,eo=v×hoμE−rr,where r denotes the position vector of the vehicle and v the velocity vector.Definen1=0,0,1T and n2=n1×ho. Then, the RAAN and the argument of perigee can be calculated as follows:
(11)Ωo=arccosn2xn2,n2y≥0,2π−arccosn2xn2,n2y<0,ωo=arccosn2⋅eon2eo,eoz≥0,2π−arccosn2⋅eon2eo,eoz<0,where n2x denotes the first component of n2 and n2y the second; eoz is the third component of eo.The trajectory optimization problem can be represented as follows:(12)minJs.t.x0,Cu,fx⇒u,where x denotes the states of the system, for example, the trajectory states of the launch vehicle; x0 the initial states; u the optimal or feasible law of control; fx the differential equations of the system, for example, the equations of motion of the launch vehicle; Cu 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.


## 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 constraintsAt the end of the ascent phase, the trajectory states should meet the specified orbital elements in order to achieve accurate orbit insertion:(13)Cequ=eotf−eof∗aotf−aof∗Ωotf−Ωof∗iotf−iof∗ωotf−ωof∗θotf−θof∗.Equation (13) can also be represented as follows:
(14)hf∗=rf∗−ae=ho∗/μE1+eo∗cosθo∗−ae,Vf∗=2μE21+eo∗cosθo∗ho∗−μEao∗,γf∗=arccos1+eo∗cosθo∗eo∗2+2eo∗cosθo∗+1.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:
(15)Cequ=htf−hf∗γtf−γf∗Vtf−Vf∗ωotf−ωof∗.(b)
Path constraintsConstraints derived from the structural strength can be represented as follows:(16)Nm≤Nm,max,Nx≤Nx,max,Qy≤Qy,max,My≤My,max,where Nm represents the normal overload, Nx the axial overload, Qy the shear force, and My the bending moment; Nm,max, Nx,max, Qy,max, and My,max represent the allowed maximum values of Nm, Nx, Qy, and My, respectively.The normal and axial overloads can be calculated as follows:(17)Nm=Psinα+Ymg,Nx=P−Dbmg,where Db 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 intos units with (s+1) discrete points, as shown in Figure 1.Figure 1
Illustration of normal load, axial load, shear force, and bending moment.Consider that the inertia force equals to the external force derived from the structure; the longitudinal load of the sectioni (denoted as Nx,i) is obtained as follows:
(18)Nx,i=P−Db−P−Dbmm0−i=P1−m0−im+m0−im−1Db=P−Db1−m0−im,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(19)m0−i=m1−is.Therefore,(20)Nx,i=isP−Dbi=0,1,⋯,s.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:(21)α≤αmax,α̇≤α̇max,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:(22)minJopt=∫Nmdt.


## 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:(23)φ=α+γ,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 Jopt (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:(24)φ=φ0+k1t,0≤t<t31,φ1+k2t−t31,t31≤t<t32,φ2+k3t−t32,t32≤t≤t3,where φ0 denotes the initial pitch angle of the third stage, φ1 the pitch angle at the time whent=t31, and φ2 the pitch angle at the time when t=t32; k1, k2, and k3 are unknown coefficients.t3 is the working time of the propulsion of the third stage; t31 and t32 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 k1, k2, and k3 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=t3) as the optimization parameters:
(25)u=φ1,φ2,φ3T.Then, unknown coefficients can be calculated as follows:(26)k1=φ3−φ2t31,k2=φ2−φ1t32−t31,k1=φ3−φ2t3−t32.The time history of the pitch angle is illustrated in Figure2.Figure 2
Schematic diagram of program angle segment design of vacuum segment.In order to find the optimalu that satisfies all the constraints, the iteration of u is detailed as follows:
(1)
Set the initial value ofu(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 statesx0, and the value of u, 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, reviseu using the optimization algorithm(6)
If the optimal solution is generated based on the convergence criterion of the optimization method, stop iteration and outputu, φt, and the trajectory; otherwise, go to step 2The optimization algorithm is proposed in Section4.


## 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:(27)vi,jG+1=wpsovi,jG+c1r1pi,j−xi,jG+c2r2pg,j−xi,jG,(28)xi,jG+1=xi,jG+vi,jG+1,where the optimal position that the ith particle has ever found is represented byPi=pi,1,pi,2,⋯,pi,DpsoT. The superscript “G” denotes the generation of the evolution, 0<wpso<1 denotes the inertia weight that indicates how much the particle is affected by its own velocity, and c1,c2∈0,4 denotes the learning factor; c1 denotes the influence of individual experience on the ith particle, and c2 the influence of group experience on the ith particle. r1,r2∈0,1 denotes uniformly distributed random numbers generated in each generation of evolution.Letζ1=c1r1+c2r2 and ζ2=c1r1pi,j+c2r2pg,j, then we have
(29)xi,jG+2+ζ1−wpso−1xi,jG+1+wpsoxi,jG=ζ2.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 (r1 and r2). 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:(30)xi,jG+2+ζ1−wpso−1xi,jG+1+wpsoxi,jG=ζ3,where
(31)ζ3=c1r1pi,j+r3ηi,jG+c2r2pg,j,ηi,jG=r4⋅RμR,σR2,where ηi,jG denotes the Gaussian perturbation term, which is added to the jth dimension of the ith particle’s position in the Gth generation; r3 and r4 are another two random numbers, and R is a random factor subject to the average value (denoted as μR) and the variance (denoted as σR2).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:(32)xi,jG=pi,j≠pg,j.The change in the velocity of the basic PSO and that of the improved PSO are as follows:(33)vi,jG+1=wpsovi,jG+c1r1pi,j−pi,j+c2r2pg,j−pi,j=wpsovi,jG+c2r2pg,j−pi,j,(34)vi,jG+1=wpsovi,jG+c1r1pi,j−pi,j+c2r2pg,j+r3ηi,jG−pi,j=wpsovi,jG+c2r2pg,j−pi,j+c2r2r3ηi,jG.Compare Equation (33) and Equation (34)), there is a deviation item:c2r2r3ηi,jG. Therefore, by proposing a Gaussian perturbation term in PSO, the particle has a larger searching space; so, the global searching capability is facilitated.
(35)xi,jG=pi,j=pg,j.In this case, Equations (33) and (34) are rewritten as follows:
(36)vi,jG+1=wpsovi,jG,vi,jG+1=wpsovi,jG+c2r2r3ηi,jG.In the basic PSO, withoutηi,jG, the searching space of the particle is quickly reduced (wpso<1) because effects of the random values (r1 and r2) do not exist. So, the particle can only search in the vicinity of xi,jG 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)
xi,jG≠pi,j, xi,jG≠pg,j, and pi,j=pg,jIn this case, Equations (33) and (34) are rewritten as follows:
(37)vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pi,j−xi,j=wpsovi,jG+c1r1+c2r2pi,j−xi,j,vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pi,j+r3ηi,jG−xi,j=wpsovi,jG+c1r1+c2r2pi,j−xi,j+c2r2r3ηi,jG.In PSO, theith particle moves to the global optimal position. But if pi,j is the local optimal position, the swarm will quickly gather there and finally obtains a suboptimal solution. But in the improved PSO, the ith particle obtains a chance of jumping out of the local optimal position.
(b)
xi,jG≠pi,j≠pg,j and pi,j≠pg,jIn this case, Equations (33) and (34) are rewritten as follows:
(38)vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pg,j−xi,j,vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pg,j+r3ηi,jG−xi,j.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:(39)G=GmaxorJpsoG−1/5∑i=04JpsoG−iJpsoG≤εJ,where JpsoG denotes the gbest in the Gth generation and εJ is the threshold.


## 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:(27)vi,jG+1=wpsovi,jG+c1r1pi,j−xi,jG+c2r2pg,j−xi,jG,(28)xi,jG+1=xi,jG+vi,jG+1,where the optimal position that the ith particle has ever found is represented byPi=pi,1,pi,2,⋯,pi,DpsoT. The superscript “G” denotes the generation of the evolution, 0<wpso<1 denotes the inertia weight that indicates how much the particle is affected by its own velocity, and c1,c2∈0,4 denotes the learning factor; c1 denotes the influence of individual experience on the ith particle, and c2 the influence of group experience on the ith particle. r1,r2∈0,1 denotes uniformly distributed random numbers generated in each generation of evolution.Letζ1=c1r1+c2r2 and ζ2=c1r1pi,j+c2r2pg,j, then we have
(29)xi,jG+2+ζ1−wpso−1xi,jG+1+wpsoxi,jG=ζ2.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 (r1 and r2). 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:(30)xi,jG+2+ζ1−wpso−1xi,jG+1+wpsoxi,jG=ζ3,where
(31)ζ3=c1r1pi,j+r3ηi,jG+c2r2pg,j,ηi,jG=r4⋅RμR,σR2,where ηi,jG denotes the Gaussian perturbation term, which is added to the jth dimension of the ith particle’s position in the Gth generation; r3 and r4 are another two random numbers, and R is a random factor subject to the average value (denoted as μR) and the variance (denoted as σR2).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:(32)xi,jG=pi,j≠pg,j.The change in the velocity of the basic PSO and that of the improved PSO are as follows:(33)vi,jG+1=wpsovi,jG+c1r1pi,j−pi,j+c2r2pg,j−pi,j=wpsovi,jG+c2r2pg,j−pi,j,(34)vi,jG+1=wpsovi,jG+c1r1pi,j−pi,j+c2r2pg,j+r3ηi,jG−pi,j=wpsovi,jG+c2r2pg,j−pi,j+c2r2r3ηi,jG.Compare Equation (33) and Equation (34)), there is a deviation item:c2r2r3ηi,jG. Therefore, by proposing a Gaussian perturbation term in PSO, the particle has a larger searching space; so, the global searching capability is facilitated.
(35)xi,jG=pi,j=pg,j.In this case, Equations (33) and (34) are rewritten as follows:
(36)vi,jG+1=wpsovi,jG,vi,jG+1=wpsovi,jG+c2r2r3ηi,jG.In the basic PSO, withoutηi,jG, the searching space of the particle is quickly reduced (wpso<1) because effects of the random values (r1 and r2) do not exist. So, the particle can only search in the vicinity of xi,jG 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)
xi,jG≠pi,j, xi,jG≠pg,j, and pi,j=pg,jIn this case, Equations (33) and (34) are rewritten as follows:
(37)vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pi,j−xi,j=wpsovi,jG+c1r1+c2r2pi,j−xi,j,vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pi,j+r3ηi,jG−xi,j=wpsovi,jG+c1r1+c2r2pi,j−xi,j+c2r2r3ηi,jG.In PSO, theith particle moves to the global optimal position. But if pi,j is the local optimal position, the swarm will quickly gather there and finally obtains a suboptimal solution. But in the improved PSO, the ith particle obtains a chance of jumping out of the local optimal position.
(b)
xi,jG≠pi,j≠pg,j and pi,j≠pg,jIn this case, Equations (33) and (34) are rewritten as follows:
(38)vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pg,j−xi,j,vi,jG+1=wpsovi,jG+c1r1pi,j−xi,j+c2r2pg,j+r3ηi,jG−xi,j.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:(39)G=GmaxorJpsoG−1/5∑i=04JpsoG−iJpsoG≤εJ,where JpsoG denotes the gbest in the Gth generation and εJ is the threshold.


## 5. Numerical Simulation

The initial trajectory states of the launch vehicle are listed in Table1. 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 t31=10 s and t32=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, Nm,max=10.0, Nx,max=3, Qy,max=400 N, and My,max=200 N·m. In addition, the terminal constraints are given in the form of target-orbital elements: eof∗=0.001532, aof∗=6770 km, Ωof∗=30°,iof∗=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.Table 1
The initial trajectory states.
VariableInitialTerminalUnitsAltitude474500kmVelocity66007669m/sFlight path angle7.60.1°Search for theu=φ1,φ2,φ3T 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.Figure 3
Time history of the velocity.Figure 4
Time history of the altitude.Figure 5
Time history of the flight path angle.Figure 6
The orbital semiaxis and eccentricity.Figure 7
Time histories of the normal load.Figure 8
Time histories of the AOA and the pitch angle.Figure 9
Time histories of the maximum axial load, shear force, and bending moment with improved PSO.Figure 10
The axial force, shear force, and bending moment at 45 s with improved PSO.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φ1=‐48.3°, φ2=‐48.4°, and φ3=‐52.5°; 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 Figures9 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.Figure 11
Time histories of the maximum axial load, shear force, and bending moment with basic PSO.Figure 12
The axial force, shear force, and bending moment at 45 s with basic PSO.To compare with basic PSO, the results in Table2.Table 2
Basic PSO compare with improved PSO.
AlgorithmsMax AOAMax rate of change of the AOAMax normal loadMax axial loadMax shear forceMax bending momentImproved PSO39.57°0.286°/s2.612.03391.6 N188.6 N·mBasic PSO46.28°0.395°/s2.772.05459.8 N224.6 N·mAccording 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.



---

*Source: 6647440-2021-02-25.xml*