# Endoatmospheric Ascent Optimal Guidance with Analytical Nonlinear Trajectory Prediction

Authors: Shilei Zhao; Wanchun Chen; Liang Yang
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

In this paper, an endoatmospheric ascent optimal guidance law with terminal constraint is proposed, which is under the framework of predictor-corrector algorithm. Firstly, a precise analytical nonlinear trajectory prediction with arbitrary Angle of Attack (AOA) profile is derived. This derivation process is divided into two steps. The first step is to derive the analytical trajectory with zero AOA using a regular perturbation method. The other step is to employ pseudospectral collocation scheme and regular perturbation method to solve the increment equation so as to derive the analytical solution with arbitrary AOA profile. The increment equation is formulated by Taylor expansion around the trajectory with zero AOA which remains the second order increment terms. Therefore, the resulting analytical solutions are the nonlinear functions of high order terms of arbitrary AOA values discretized in Chebyshev-Gauss-Legendre points, which has high accuracy. Secondly, an iterative correction scheme using analytical gradient is proposed to solve the endoatmospheric ascent optimal guidance problem, in which the dynamical constraint is enforced by the resulting analytical solutions. It only takes a fraction of a second to get the guidance command. Nominal simulations, Monte Carlo simulations, and optimality verification are carried out to test the performance of the proposed guidance law. The results show that it not only performs well in providing the optimal guidance command, but also has great applicability, high guidance accuracy and computational efficiency. Moreover, it has great robustness even in large dispersions and uncertainties.


---

## Body

## 1. Introduction

For rocket-powered launch vehicles, the guidance accuracy of ascent phase has a significant effect on the flight performance in subsequent missions [1]. In this phase, it involves strong nonlinearity and coupling dynamics resulting from the complicated forces such as thrust, aerodynamic forces, and gravity. That brings some troublesome but intriguing challenges for the ascent guidance design if the terminal trajectory constraint should be enforced strictly. Plenty of efforts have been made to this research field.Optimal control is widely used for ascent guidance and trajectory tracking [2–4]. The majority methods can be divided into two categories: indirect method and direct method [5]. The indirect method gets the optimal solution by solving two-point boundary-value problem (TPBVP). In [6], Lu et al. proposed a closed-loop endoatmospheric ascent guidance law by using the classical finite difference method for solving TPBVP. Further research and development on improving this method are referred to [7–9]. In [10], Lu et al. presented a numerical trajectory reconstruction algorithm based on a finite element method, which satisfies the real-time requirement of generating the guidance commands. In [11, 12], Dukeman et al. proposed an endoatmospheric ascent guidance for rocket-powered launch vehicles, in which multiple shooting method is employed to solving TPBVP with a high computational efficiency. In addition, many optimization algorithms such as convex optimization [13] and swarming algorithms [14–16] are also used to solve TPBVP resulting from ascent guidance problem. But the indirect method needs to derive the complicated necessary conditions including adjoint equations, transversality conditions, and maximum principle. And it is sensitive to the initial guess, which leads to it being unsuitable for onboard guidance. The direct method transforms the optimal control problem into a nonlinear programming problem (NLP), and then well-developed numerical algorithms can be used to solve it in a discrete manner. In [17, 18], Ghose et al. adopted a differential evolution method to obtain the optimal ascent phase trajectory for the hypersonic vehicle. In recent years, benefiting from the high accuracy in both the primal and dual solutions and the equivalence between the direct and indirect forms, the pseudospectral method has been widely used in trajectory optimization and guidance design for the ascent phase [19–27]. However, it requires the third-party software to solve the complex NLP problem, which has huge codes and occupies large computational sources. Therefore, it is not suitable for onboard application. In order to alleviate the computing burden of online guidance, several other methods have been applied. The predictor-corrector method has great potential for ascent phase guidance. In [28], Li et al. proposed a novel segmented and weighted adaptive predictor-corrector guidance method for the ascent phase of hypersonic vehicle. In [29], Yang et al. proposed a linear Gauss pseudospectral model predictive control for solving the nonlinear optimal control problem without solving NLP, which has high computational efficiency. It has been successfully applied into entry guidance for the hypersonic vehicle [30–32].However, all of the mentioned methods heavily rely on numerical integration prediction, which leads to the heavy computation and time-consuming. On the contrary, if the trajectory prediction can be provided in an analytical manner, the computational efficiency can be significantly improved. For the gliding trajectory of hypersonic vehicle, the analytical solution has been widely studied and applied to the entry guidance [33–38]. But there is no study on precise analytical solution to ascent trajectory. The existing ascent guidance laws in an analytic manner ignore the aerodynamic force, which results in large prediction error as for the endoatmospheric trajectory [39, 40]. In endoatmospheric ascent phase, the variations in velocity, altitude, and flight-path angle are much great, which leads to the great difference in a flight environment. It also involves the controlled thrust and aerodynamics force which results that its dynamics has strong state coupling and nonlinearity. Therefore, it is difficult to find the analytical solution using traditional method.This paper aims at developing an efficient optimal ascent guidance law which considers multiple constraints and quadratic performance index. This method is under the framework of predictor-corrector algorithm. Firstly, an analytical trajectory prediction with high accuracy is proposed using a regular perturbation method. The solving process consists of two steps. In the first step, the analytical trajectory with zero Angle of Attack (AOA) is derived. In the derivation process, an assumptive constant axial force is introduced to decompose the dynamics, and therefore, the resulting subdynamics is uncoupled and can be solved by integrating analytically. In the second step, Taylor expansion around the analytical trajectory with zero AOA is used to formulate the increment dynamic equation which fully considers the second increment terms. And then, pseudospectral collocation scheme and regular perturbation method are employed to derive the analytical solution of the terminal state, which is a nonlinear function of high order terms of arbitrary AOA values discretized in Chebyshev-Gauss-Legendre points. Secondly, by introducing a quadric performance index, an endoatmospheric ascent optimal guidance problem is formulated, in which the dynamical constraint is enforced by the resulting analytical solutions. It should be noted that it is impossible to derive the analytical solution to such problem because the high order terms are involved. Then, an iterative solving scheme using analytical gradient is proposed to solve it. This method is attractive from the point of high efficiency, high accuracy with few points, and character that the optimal guidance command for nonlinear guidance problem is provided. The most difference from the previous methods is that the analytical prediction is involved to remove the trajectory integration, which significantly improves the computational efficiency. Furthermore, the high order increment terms are fully considered to bring the high accuracy for analytical prediction. Finally, the accuracy verification for analytical solution, several nominal simulations, comparison with the optimal solution, and Monte Carlo simulations have been carried out to evaluate the performance of the proposed method. The results show that the analytical solutions are highly consistent with the numerical solutions. The proposed method not only performs well in providing the optimal guidance command but also has high computational efficiency and guidance accuracy. Moreover, it has great robustness even in large dispersions and uncertainties.The remainder of this paper is organized as follows: Section2 establishes the formulation of ascent phase guidance. Section 3 derives the analytical solution to ascent trajectory. Section 4 provides the optimal AOA to eliminate the terminal errors. And then, the implementation of ascent guidance law using optimal analytical correction is presented in Section 5. Finally, Section 6 provides some simulations to evaluate the performance of the proposed guidance law.


## 2. Problem Formulation

### 2.1. Ascent Dynamics

During the ascent phase for a point-mass vehicle model, the equations of motion for the vehicle over the nonrotating spherical Earth in the longitudinal plane are established as(1)V̇=Pcosα−Dm−gsinγ,γ̇=Psinα+LmV+Vr−gVcosγ,ḣ=Vsinγ,where m, V, γ, h, and α are the mass, velocity, flight-path angle, altitude, and AOA, respectively; r=Re+h is the distance from the center of Earth to the vehicle, where Re is the average radius of the Earth, and g=μ/r2 is the gravitational acceleration; P is the thrust of the vehicle; and L,D are the aerodynamic drag and lift expressed as
(2)L=12ρV2SCl,D=12ρV2SCd,where ρ is the atmospheric density and S is the reference area. The terms, Cl and Cd, are the lift and drag coefficients which are dependent on Mach number and AOA.


### 2.2. Terminal Constraints and Performance Index

In order to guarantee the successful handover with the subsequent phase, some terminal state, mainly flight-path angle and altitude, must be enforced to the desired ones in the ascent phase. For simplification, the terminal constraints are represented in a general form as(3)ψxtf=γf−γfc=0,hf−hfc=0,where γf and hf are the ascent phase terminal flight-path angle and altitude and γfc and hfc are the desired ones, respectively.For the ascent phase guidance, the control variable is AOA, which determines the direction of the thrust. And the autopilot is used to track the AOA profile to achieve the flight mission. In order to provide the stable input, the weighted squared sum of AOA is selected as the performance index.(4)J=12∫t0tfα2tf−tndt.It is noted that the denominator,tf−tn, is used to shape the AOA profile and makes the terminal AOA zero if n>0 is selected.


### 2.3. Vehicle Description

For rocket-powered launch vehicle, its thrust can be formulated as(5)P=Ispqmg0−SePa,where Isp is specific impulse and qm is the mass flow which can be treated as a constant for the solid launch vehicle; Se is the area of nozzle; and Pa is the atmospheric pressure at the flight altitude. Noted that, although the thrust is related to altitude, the change resulting from the altitude is very small. Additionally, the lift and drag coefficients can be fitted to be linear function and parabolic function of AOA due to the axisymmetric aerodynamic configuration, which are formulated as
(6)Cl=Clαα,Cd=Cd0+Cdαα2,where Clα, Cd0, and Cdα are all related to Mach number, which can be fitted as the polynomial of Mach number, i.e.,
(7)Clα=∑i=0pCl,iαMai,Cd0=∑i=0pCd0,iMai,Cdα=∑i=0pCd,iαMai.In the solving of the analytical solutions, these coefficients are regarded as the constants,C¯lα, C¯d0, and C¯dα, which are the average values during the flight. Because the aerodynamic coefficients do not vary greatly with the Mach when the vehicle is flying at supersonic speed, the error caused by this simplified hypothesis is small. The second stage of a two-stage launch vehicle model, the parameters of which listed in Table 1, is used for all simulations.Table 1
Parameter of the vehicle model.
ParameterValueParameterValueCl,iα[3.774, 0.451, -0.075,4.98e−03, −1.19e−04]m0kg46680Cd0,i[0.846, -0.152,1.58e−03, −5.67e−04]qmkg/s479.9Cd,iα[4.963, 0.598, -0.097,7.12e−03, −1.93e−04]IspkN264.7S (m2)4.3tfs60


## 2.1. Ascent Dynamics

During the ascent phase for a point-mass vehicle model, the equations of motion for the vehicle over the nonrotating spherical Earth in the longitudinal plane are established as(1)V̇=Pcosα−Dm−gsinγ,γ̇=Psinα+LmV+Vr−gVcosγ,ḣ=Vsinγ,where m, V, γ, h, and α are the mass, velocity, flight-path angle, altitude, and AOA, respectively; r=Re+h is the distance from the center of Earth to the vehicle, where Re is the average radius of the Earth, and g=μ/r2 is the gravitational acceleration; P is the thrust of the vehicle; and L,D are the aerodynamic drag and lift expressed as
(2)L=12ρV2SCl,D=12ρV2SCd,where ρ is the atmospheric density and S is the reference area. The terms, Cl and Cd, are the lift and drag coefficients which are dependent on Mach number and AOA.


## 2.2. Terminal Constraints and Performance Index

In order to guarantee the successful handover with the subsequent phase, some terminal state, mainly flight-path angle and altitude, must be enforced to the desired ones in the ascent phase. For simplification, the terminal constraints are represented in a general form as(3)ψxtf=γf−γfc=0,hf−hfc=0,where γf and hf are the ascent phase terminal flight-path angle and altitude and γfc and hfc are the desired ones, respectively.For the ascent phase guidance, the control variable is AOA, which determines the direction of the thrust. And the autopilot is used to track the AOA profile to achieve the flight mission. In order to provide the stable input, the weighted squared sum of AOA is selected as the performance index.(4)J=12∫t0tfα2tf−tndt.It is noted that the denominator,tf−tn, is used to shape the AOA profile and makes the terminal AOA zero if n>0 is selected.


## 2.3. Vehicle Description

For rocket-powered launch vehicle, its thrust can be formulated as(5)P=Ispqmg0−SePa,where Isp is specific impulse and qm is the mass flow which can be treated as a constant for the solid launch vehicle; Se is the area of nozzle; and Pa is the atmospheric pressure at the flight altitude. Noted that, although the thrust is related to altitude, the change resulting from the altitude is very small. Additionally, the lift and drag coefficients can be fitted to be linear function and parabolic function of AOA due to the axisymmetric aerodynamic configuration, which are formulated as
(6)Cl=Clαα,Cd=Cd0+Cdαα2,where Clα, Cd0, and Cdα are all related to Mach number, which can be fitted as the polynomial of Mach number, i.e.,
(7)Clα=∑i=0pCl,iαMai,Cd0=∑i=0pCd0,iMai,Cdα=∑i=0pCd,iαMai.In the solving of the analytical solutions, these coefficients are regarded as the constants,C¯lα, C¯d0, and C¯dα, which are the average values during the flight. Because the aerodynamic coefficients do not vary greatly with the Mach when the vehicle is flying at supersonic speed, the error caused by this simplified hypothesis is small. The second stage of a two-stage launch vehicle model, the parameters of which listed in Table 1, is used for all simulations.Table 1
Parameter of the vehicle model.
ParameterValueParameterValueCl,iα[3.774, 0.451, -0.075,4.98e−03, −1.19e−04]m0kg46680Cd0,i[0.846, -0.152,1.58e−03, −5.67e−04]qmkg/s479.9Cd,iα[4.963, 0.598, -0.097,7.12e−03, −1.93e−04]IspkN264.7S (m2)4.3tfs60


## 3. Analytical Ascent-Trajectory Prediction

In this section, the analytical ascent trajectory prediction is developed. Firstly, a simplified condition, in which the AOA is set to zero, is considered. Under the simplification, the analytical solution to ascent trajectory is only related to the flight time for the given initial values and vehicle parameters. Next, based on the presented analytical solution with zero AOA, the analytical solution to trajectory increment corresponding to nonzero AOA is derived. Then, the analytical solution with arbitrary AOA profile, which is the sum of analytical solution with zero AOA and trajectory increment corresponding to nonzero AOA, is provided. In the derivation of these analytical solutions, the regular perturbation method is adopted to decompose the coupled dynamics into a Taylor series which can be solved sequentially.


### 3.1. Dynamics for Analytical Trajectory

In order to simplify the solutions, the index atmospheric density model,ρ=ρ0exp−βh, where ρ0=1.225kg/m3 is the sea-level atmospheric density and β=1.389×10−4m−1, is adopted for the solving of the analytical solutions.In order to avoid confusion, the variables corresponding to zero AOA are expressed with subscript “b.” For the case of zero AOA, it is obvious that P¯sinα+L=0 and D=ρbVb2SC0/2; then, the dynamics is simplified to
(8)V̇b=1mP¯−ρbVb2SC¯d02−gbsinγb,γ̇b=Vbrb−gbVbcosγb,ḣb=Vbsinγb.For simplifying the derivation, a variable related toγb is introduced.
(9)θb=ln1+sinγb1−sinγb.Replacingγb with θb as the state variable, the equation of motion presented in (8) can be rewritten as
(10)V̇b=1mP¯−ρbVb2SC¯d02−gbsinγb,θ̇b=2Vbrb−gbVb,ḣb=Vbsinγb.The increment of the velocity, flight-path angle, and altitude with nonzero AOA relative to that with zero AOA are defined asΔV, Δγ, and Δh, which are used to derive the analytical solutions with arbitrary AOA profile. It is obvious that the differential equations of those trajectory increments are the difference of (1) and (8). For the sake of simplification, the following Taylor’s expansions are adopted.
(11)exp−βΔh≈1−βΔh,sinα≈α,1V≈1Vb−ΔVVb2,cosα≈1−α22,sinγ≈sinγb+cosγbΔγ−sinγbΔγ22,cosγ≈cosγb−sinγbΔγ−cosγbΔγ22.Then, the dynamics equations of increments are established as(12)ΔV̇=−12mP¯+ρbVb2SC¯dαα2+2ρbVbSCdΔV−βρbVb2SCdΔh−gicosγbΔγ+12gisinγbΔγ2,Δγ̇=1mP¯Vb+12ρbVbSC¯lα1−βΔhα+1r0+giVb2cosγbΔV−Vbr0−giVbsinγb+0.5cosγbΔγΔγ,Δḣ=sinγb+cosγbΔγΔV+VbcosγbΔγ−sinγbΔγ22.Becauser≫Δh, the altitude increment on the value of 1/r is ignored. It should be noted that the second order terms related to the AOA and state increment are remained to achieve high accuracy.


### 3.2. Analytical Ascent Trajectory with Zero AOA

In (10), the right side of the differential equation of Vb can be divided into the following two parts:
(13)V̇b=Favem+Fεm,where Fave is the major axial force, which is the sum of thrust, average atmospheric drag, and average axial component of gravity, and Fε is the perturbation axial force:
(14)Fave=P¯−ρaveVave2SC¯d02−mavegavesinγave,Fε=ρaveVave2SC¯d02+mavegavesinγave−ρbVb2SC¯d02−mgbsinγb,where mave,γave, Vave, ρave, and gave are the average mass, flight-path angle, velocity, atmospheric density, and gravitational acceleration, respectively, which can be previously determined by the average of its initial and terminal values. In the boost phase, the thrust is usually much greater than the change of aerodynamic drag and the axial component of gravity. Thus, the perturbation axial force is small enough and can be regarded as the correction terms. On the other hand, the trigonometric term sinγb can be written as follows:
(15)sinγb=sinγ0+c1Δθb+sinγb−sinγ0−c1Δθb,Δθb=θb−θ0,where c1 is the constant coefficient, which can be obtained by the initial and terminal values of γb and θb:
(16)c1=sinγbf−sinγ0θbf−θ0.In the boost phase, the flight-path angle varies slightly; therefore, the linear termsinγ0+c1Δθb is much larger than the rest terms. Therefore, according to the regular perturbation method presented in Appendix A, (10) can be modified as
(17)V̇bθ̇bḣb=Favem2Vbrb−gbVbVbsinγ0+c1Δθb+εkFεm0Vbsinγb−sinγ0−c1Δθb,where ε is a small parameter for the subsequent regular perturbation order and k is a constant, which is set to be equal to ε. According to (A.5), the zeroth-order dynamics is
(18)V̇b0=Favem,θ̇b0=2Vb0r0−giVb0,ḣb0=Vb0sinγ0+c1Δθb0.It is noted that, in order to decouple the differential equations ofθb0 and hb0, the terms related to r are expanded at initial point since the change of r is much smaller than its initial value, r0. The first-order dynamics is
(19)V̇b1=1mP¯−Fave−12ρ0exp−βhb0Vb02SC¯d0−gavesinγ0+c1Δθb0,θ̇b1=−2Vb0r0−2giVb0Δhb0r0+2Vb0r0+giVb0Vb1Vb0,ḣb1=Vb0sinγb0−sinγ0−c1Δθb0+Vb1sinγb0+Vb0cos2γb0θb12,where γb0 is the zeroth-order solution to flight-path angle, which can be derived from θb0.
(20)γb0=sin−1expθb0−1expθb0+1.It should be noted that the terms related tohb1 in the differential equations of θb1 are ignored due to hb1≪Δhb0. Additionally, according to the boundary conditions, the initial values of zeroth-order and first-order states are
(21)Vb00=V0,θb00=θ0,hb00=h0,Vb01=0,θb01=0,hb01=0.


#### 3.2.1. Zeroth-Order Solution

It is obvious that the three differential equations in (18) are not coupled, so that the analytical solution to Vb0, θb0, and hb0 can be derived by integrating these three differential equations in sequence.Integrating the first equation of (18) from t0 to t, the analytical solution of Vb0 at time t can be easily obtained as
(22)Vb0=V0−Veln1−t−t0T,where T=m0/qm and Ve=Fave/qm. Obviously, the differential equations of other variables are all related to Vb0; thus, the analytical solutions of these variables are the functions of Vb0. For the sake of simplifying the expression, a dimensionless variable, V¯=Vb0/Ve, is employed, which is regarded as an independent variable in the subsequent derivation. From (22), the expression of V¯ is
(23)V¯=V¯0−ln1−t−t0T,where V¯0=V0/Ve is initial value of V¯. Substituting (22) into the first equation of (18), the differential equation of V¯ can be expressed as a simple function.
(24)V¯̇=FavemVe=expV¯−V¯0T.Dividing the second and third equations of (18) by (24) yields
(25)dθb0dV¯=2K1V¯−K2V¯exp−V¯,(26)dhb0dV¯=K1r0sinγ0V¯+c1V¯Δθb0exp−V¯,where K1=TVeexpV¯0/r0and K2=giTexpV¯0/Ve. Integrating (25) from V¯0 to V¯, θb0 can be obtained as
(27)θb0V¯=θ0+2K1fV¯V¯;1,1−2K2ξV¯V¯;1,1,where fV¯V¯;n,mn≥0,m≥0 is
(28)fV¯V¯;n,m=V¯n+1−V¯0n+1n+1,m=0,n!mn+1∑k=0nmkk!V¯0kexp−mV¯0−V¯kexp−mV¯,m>0.ξV¯V¯;n,mn>0,m>0 is the integral term expressed as
(29)ξV¯V¯;n,m=∫V¯0V¯exp−mxV¯xV¯ndxV¯.Although the integral term is quite simple, it has no analytical solution. In order to derive its approximate analytical solution, Chebyshev interpolation polynomial is used to approximate the integrand. Let(30)exp−mxV¯xV¯n≈∑i=0Nexp−mxV¯ixV¯inLix=∑i=0NpV¯in,mxi,where xV¯i is the value of V¯ at the i-th interpolation point xi; Lixi=0,⋯,N is the interpolation basis function; and x=xV¯−V¯2/V¯1 is the dimensionless variable with −1,1, where V¯1=V¯−V¯0/2 and V¯2=V¯+V¯0/2. Then, ξV¯ can be approximated as
(31)ξV¯V¯;n,m=12V¯−V¯0∑i=0NWixexp−mxV¯ixV¯in,where Wix is the integral weight coefficient. The derivation process is presented in Appendix B in detail. The error analysis is carried out by comparing the solution obtained by numerical integration and the 6th-degree interpolation polynomials in Figure 1, from which we can see that the interpolation polynomials have high accuracy.Figure 1
Comparison ofξV¯V¯ with various n and m.In order to derive the analytical solution tohb0, an integral related to Δθb0 is introduced.
(32)IθxV¯;n,m=∫V¯0xV¯Δθb0xV¯nexp−mxV¯dxV¯,n≥0,m>0.Using the formula of partial integration, the analytical solution toIθxV¯;n,m is expressed as
(33)IθxV¯;n,m=−n!mn+1∑k=0nmkk!Δθb0xV¯xV¯kexp−mxV¯−θrxV¯;k,m,where θrxV¯;n,m is the integral shown below:
(34)θrxV¯;n,m=∫θ0θb0xV¯nexp−mxV¯dθb0;n≥0,m>0.It is obvious that the analytical solution toθr is similar to that of θb0. Substituting (25) into (34), θr can be expressed as
(35)θrxV¯;n,m=2∫V¯0xV¯K1xV¯n+1exp−m+1xV¯−K2xV¯n−1exp−m+1xV¯dxV¯=2K1fV¯xV¯;1,m+1−2K2ξV¯xV¯;1,m+1,n=0,2K1fV¯xV¯;n+1,m+1−2K2fV¯xV¯;n−1,m+1,n>0.Then, integrating (26) from V¯0 to V¯, the analytical solution to hb0 is
(36)hb0V¯=h0+K1r0sinγ0fV¯V¯;1,1+c1IθV¯;1,1,where
(37)IθV¯;1,1=2K1fV¯V¯;1,2+fV¯V¯;2,2−1+V¯exp−V¯fV¯V¯;1,1−2K2fV¯V¯;0,2+ξV¯V¯;1,2−1+V¯exp−V¯ξV¯V¯;1,1.Finally, the zeroth-order solutions with zero AOA are derived, which are the fundamental for deriving the first-order correction.


#### 3.2.2. First-Order Correction

In this subsection, let us focus on the first-order dynamics. Regarding the dimensionless velocity,V¯, as the dependent variable, the first-order dynamics can be formulated as
(38)dVb1dV¯=VeKP−1−KaC¯d0V¯2exp−β¯Δh¯b0−K2sinγ0+c1Δθb0exp−V¯,dθb1dV¯=−2K1V¯−2K2V¯exp−V¯Δh¯b0+2K1+K2V¯2exp−V¯V¯b1,dhb1dV¯=K1r0exp−V¯V¯sinγb0−c1Δθb0−sinγ0+V¯b1sinγb0+V¯cos2γb0θb12,where Δh¯b0=hb0−h0/r0 and V¯b1=Vb1/Ve are the dimensionless variables and β¯=βr0,KP=P¯/Fave, and Ka=ρiSVe2/2Fave are the constant coefficients.The analytical solution toVb1 can be derived by integrating the first equation of (38). But the term, which is the integration of V¯2exp−β¯Δh¯b0, cannot be solved analytically because Δh¯b0 is not the linear function of V¯. Similar to (30), let
(39)fρxV¯=xV¯2exp−β¯Δh¯b0≈∑i=0NfρxV¯iLix,where fρxV¯i=xV¯i2exp−β¯K1sinγ0fV¯xV¯i;1,1+c1IθxV¯i;1,1. As shown in Figure 2, the approximation polynomials are accurate enough in comparison with the original functions fρxV¯. Then, the analytical solution to Vb1 is
(40)Vb1V¯=VeKP−1V¯−V¯0−12V¯−V¯0KaC¯d0∑i=0NfρxV¯iWix−K2sinγ0exp−V¯0−exp−V¯−K2c1IθV¯;0,1,where according to (33), the expression of IθV¯;0,1 is
(41)IθV¯;0,1=2K1fV¯V¯;1,2−exp−V¯fV¯V¯;1,1−2K2ξV¯V¯;1,2−exp−V¯ξV¯V¯;1,1.Figure 2
Comparison betweenfρV¯ and the corresponding approximate polynomial.By integrating the second equation of (38), the analytical solution to θb1 can be obtained as
(42)θb1V¯=−2K1h¯r00+h¯r10−Δh¯b01+V¯exp−V¯+2K1V¯r1−V¯b1exp−V¯+4K2∫V0V¯Δh¯b0exp−xV¯xV¯dxV¯+2K2∫V0V¯V¯b1exp−xV¯xV¯2dxV¯,where h¯r00, h¯r10, and V¯r1 are all the functions of V¯, which are expressed as
(43)h¯r00V¯=K1sinγ0fV¯V¯;1,2+c1IθV¯;1,2,h¯r10V¯=K1sinγ0fV¯V¯;2,2+c1IθV¯;2,2,V¯r1V¯=KP−1fV¯V¯;0,1−12V¯−V¯0KaC¯d0∑i=0NfρxV¯iexp−xV¯iWix−K2sinγ0fV¯V¯;0,2−K2c1IθV¯;0,2,where the expression of IθV¯;0,2, IθV¯;1,2, and IθV¯;1,2 can be obtained from (33).As for the two integration terms in (42), there are no analytical solutions. Chebyshev interpolation polynomials are also adopted to approximate the integrands. Let
(44)exp−mxV¯xV¯n=∑i=0NcV¯in,mxV¯i,where the coefficient cV¯in,m is
(45)cV¯in,m=∑k=iNpV¯kn,mCnk−V¯2/V¯1n−kV¯1k.Then, the analytical solutions to the two integration terms are(46)∫V0V¯Δh¯b0exp−xV¯xV¯dxV¯=∑i=0NcV¯i1,1IhV¯;i,∫V0V¯V¯b1exp−xV¯xV¯2dxV¯=∑i=0NcV¯i2,1IV¯V¯;i,where
(47)IhV¯;n=∫V¯0V¯Δh¯b0xV¯nexp−xV¯dxV¯,IV¯V¯;n=∫V¯0V¯V¯b1xV¯nexp−xV¯dxV¯.By the formula of partial integration, the analytical solution toIhV¯;n and IV¯V¯;n can be derived as
(48)IhV¯;n=−n!∑k=0n1k!Δh¯b0V¯kexp−V¯−ϑhV¯;k,IV¯V¯;n=−n!∑k=0n1k!V¯b1V¯kexp−V¯−ϑV¯V¯;k,where the expressions of ϑhV¯;n and ϑV¯V¯;n are
(49)ϑhV¯;n=K1sinγ0fV¯V¯;n+1,2+c1IθV¯;n+1,2,ϑV¯V¯;n=KP−1fV¯V¯;n,1−12V¯−V¯0KaCd0∑i=0NfρxV¯ixV¯inexp−xV¯iWiV¯−K2sinγ0fV¯V¯;n,2−K2c1IθV¯;n,2.Then, the analytical solution toθb1 is
(50)θb1V¯=−2K1h¯r00+h¯r10−Δh¯b01+V¯exp−V¯+4K2∑i=0NcV¯i1,1IhV¯;i+2K1V¯r1−V¯b1exp−V¯+2K2∑i=0NcV¯i2,1IV¯V¯;i.As for the analytical solution tohb1, it is the integral of dhb0/dV¯ from V¯0 to V¯. It is noted that all the variables in dhb0/dV¯, such as γb0, θb0, V¯b1, and θb1, are expressed as the function of V¯. Therefore, we can use Chebyshev interpolation polynomial to approximate the integrand. Let
(51)fhxV¯=exp−xV¯xV¯sinγb0−c1Δθb0−sinγ0+V¯b1sinγb0+xV¯cos2γb0θb12≈∑i=0NfhxV¯iLix.Then, the solution tohb1 can be obtained as
(52)hb1V¯=K1r0∑i=0NfhxV¯iWiV¯.Finally, the first-order solutions are derived. According to (A.6), the analytical solution with zero AOA is the sum of the zero-order and first-order solutions.


### 3.3. Analytic Ascent Trajectory with Arbitrary AOA Profile

In the last subsection, the analytical solutions to ascent trajectory with zero AOA are provided. In order to obtain the full analytical solutions to ascent trajectory with arbitrary AOA profile, the analytical solutions to trajectory increment corresponding to no-zero AOA are presented firstly in this subsection. Then, the full analytical solutions are the sum of analytical trajectory with zero AOA and trajectory increment corresponding to nonzero AOA.


#### 3.3.1. Regular Perturbation Model with AOA Increment

Because the independent variable for the analytical solution with zero AOA is the dimensionless velocity,V¯, the differential equations of trajectory increments should be formulated as
(53)dΔVdV¯=VefV1α2+fV2ΔVVe+fV3Δγ+fV4Δhr0+fV5Δhα2r0+fV6Δγ2,dΔγdV¯=fγ1α+fγ2ΔVVe+fγ3Δγ+fγ4Δhαr0+fγ5Δγ2,dΔh¯dV¯=r0fh1ΔVVe+fh2Δγ+fh3ΔVΔγVe+fh4Δγ2,,where fV1∼fV6, fγ1∼fγ5, and fh1∼fh4 are all the functions of V¯, which are expressed as
(54)fV1V¯=−KP2+KaC¯dαV¯b2exp−βΔhb,fV2V¯=−2KaC¯d0V¯bexp−βΔhb,fV3V¯=−K2cosγbexp−V¯,fV4V¯=βr0KaCd0V¯b2exp−βΔhb,fV5V¯=βr0KaC¯dαV¯b2exp−βΔhb,fV6V¯=K2sinγbexp−V¯2,fγ1V¯=KPV¯b+KaC¯lαV¯bexp−βΔhb,fγ2V¯=K1+K2V¯b2cosγbexp−V¯,fγ3V¯=−K1V¯b−K2V¯bsinγbexp−V¯,fγ4V¯=−βr0KaC¯lαV¯bexp−βΔhb,fγ5V¯=−K1V¯b−K2/V¯bcosγbexp−V¯2,fh1V¯=K1sinγbexp−V¯,fh2V¯=K1V¯bcosγbexp−V¯,fh3V¯=K1cosγbexp−V¯,fh4V¯=−K1V¯bsinγbexp−V¯2..Generally, the AOA profile is expressed as the function of the flight time. The AOA profile with the dimensionless velocity can be formulated by using the relationship between the flight timet and V¯, which can be easily expressed as
(55)t=t0+T1−expV¯0−V¯.In the differential equation ofΔγ, it is easy to find that the term, fγ1α, is much larger than the sum of other terms, fγ1α≫fγ2ΔV/Ve+fγ3Δγ+fγ4Δhα/r0+fγ5Δγ2, which is the dominant term. Similarly, the term, fV1α2+fV3Δγ, is the dominant term of the differential equation of ΔV. As for the differential equation of Δh, due to Vb≫ΔV and Δγ2≤1, it is obvious that the term, fh1ΔV/Ve+fh2Δγ, is the dominant term. Therefore, the dynamics can be modified as
(56)dΔγdV¯dΔVdV¯dΔhdV¯=fγ1αVefV1α2+VefV3Δγr0Vefh1ΔV+r0fh2Δγ+εk1Vefγ2ΔV+fγ3Δγ+1r0fγ4Δhα+fγ5Δγ2fV2ΔV+Ver0fV4Δh+Ver0fV5Δhα2+VefV6Δγ2r0Vefh3ΔVΔγ+r0fh4Δγ2.According to the regular perturbation method presented in AppendixB, the zeroth-order and first-order dynamics are
(57)dΔγ0dV¯=fγ1V¯α,dΔV0dV¯=VefV1V¯α2+fV3V¯Δγ0,dΔh0dV¯=r0fh1V¯ΔV0Ve+fh2V¯Δγ0,,(58)dΔγ1dV¯=fγ2V¯ΔV0Ve+fγ3V¯Δγ0+fγ4V¯Δh0αr0+fγ5V¯Δγ02,dΔV1dV¯=VefV2V¯ΔV0Ve+fV4V¯Δh0r0+fV5V¯Δh0α2r0+fV6V¯Δγ02,dΔh1dV¯=r0fh1V¯ΔV1Ve+fh2V¯Δγ1+fh3V¯ΔV0Δγ0Ve+fh3V¯Δγ02..The initial values of zeroth-order and first-order states are all zeros.


#### 3.3.2. Zeroth-Order Solution to Trajectory Increment

It is obvious that the three differential equations of zero-order states are not coupled. Therefore, by integrating these differential equations orderly, the analytical solutions toΔγ0, ΔV0, and Δh0 can be obtained as
(59)Δγ0V¯=∫V¯0V¯fγ1xV¯αxV¯dxV¯,(60)ΔV0V¯=Ve∫V¯0V¯fV1xV¯α2xV¯dxV¯+Ve∫V¯0V¯fV3xV¯Δγ0xV¯dxV¯,(61)Δh0V¯=r0Ve∫V¯0V¯fh1xV¯ΔV0xV¯dxV¯+r0∫V¯0V¯fh2xV¯Δγ0xV¯dxV¯.Although the integral terms in the above formulas are quite complex, the integrands are related only to the independent variable,V¯. Here, Chebyshev polynomials are also used to approximate the integrands. Let
(62)fγ1xV¯αxV¯≈∑i=0Nfγ1xV¯iαxV¯iLix,(63)fV1xV¯α2xV¯≈∑i=0NfV1xV¯iα2xV¯iLix,(64)fV3xV¯Δγ0xV¯≈∑i=0NfV3xV¯iΔγ0xV¯iLix,(65)fh1xV¯ΔV¯0xV¯≈∑i=0Nfh1xV¯iΔV¯0xV¯iLix,(66)fh2xV¯Δγ0xV¯≈∑i=0Nfh2xV¯iΔγ0xV¯iLix.Substituting (62) into (59), the analytical solution to the zeroth-order flight-path angle increment can be obtained as
(67)Δγ0V¯=12V¯−V¯0∑i=0Nfγ1xV¯iαxV¯iWix=γα10V¯α,where α is the vector consisting of the value of AOA at interpolating points, i.e.,
(68)α=αxV¯0,αxV¯1,⋯,αxV¯NT,γα10V¯ is the coefficient matrix expressed as
(69)γα10V¯=12V¯−V¯0W0V¯fγ1xV¯0,⋯,WNV¯fγ1xV¯N=WTV¯Fγ1,where
(70)WV¯=12V¯−V¯0W0x,⋯,WNxT.Fγ1 is the diagonal matrix whose diagonal elements are the value of fγ1 at interpolating points, i.e.,
(71)Fγ1=fγ1xV¯0⋯0⋮⋱⋮0⋯fγ1xV¯N.Similarly,Fγ2∼Fγ3, FV1∼FV4, and Fh1∼Fh2 mentioned below are all the diagonal matrices whose diagonal elements are the value of corresponding function at interpolating points.Substituting (63) and (64) into (60), the analytical solution to the zeroth-order velocity increment can be obtained as
(72)ΔV0V¯=12V¯−V¯0Ve∑i=0NfV1xV¯iα2xV¯i+fV3xV¯iΔγ0xV¯iWix=12V¯−V¯0Ve∑i=0NfV1xV¯iα2xV¯i+fV3xV¯iγα10xV¯iαWix=Vα10V¯α+Vα20V¯αsqu,where αsqu is the vector consisting of the square of AOA at interpolating points, i.e.,
(73)αsqu=α2xV¯0,α2xV¯1,⋯,α2xV¯NT,Vα10V¯ and Vα20V¯ are the coefficient matrices expressed as
(74)Vα10V¯=VeWTV¯FV3AcbFγ1,Vα20V¯=VeWTV¯FV1,where Acb is
(75)Acb=WxV¯0,⋯,WxV¯NT.Substituting (65) and (66) into (61), the analytical solution to the zeroth-order altitude increment can be obtained as
(76)Δh0V¯=12V¯−V¯0r0∑i=0Nfh1xV¯iΔV¯0xV¯i+fh2xV¯iΔγ0xV¯iWix=12V¯−V¯0r0∑i=0Nfh1xV¯iVα10V¯α+Vα20V¯αsquVeWix+12V¯−V¯0r0∑i=0Nfh2xV¯iWTV¯Fγ1αWix=hα10V¯α+hα20V¯αsqu,where
(77)hα10V¯=r0WTV¯Fh1AcbFV3+Fh2AcbFγ1,hα20V¯=r0WV¯Fh1AcbFV1..


#### 3.3.3. First-Order Correction to Trajectory Increment

As for the first-order corrections, their analytical solutions can be obtained by integrating the three differential equations of (58) in sequence. Similar to the zeroth-order solution, Chebyshev polynomials are also used to approximate the integrands. Therefore, the analytical solution to Δγ1 is
(78)Δγ1V¯=12V¯−V¯0∑i=0NWixfγ2xV¯iΔV0xV¯iVe+fγ3xV¯iΔγ0xV¯i+12V¯−V¯0∑i=0NWixfγ4xV¯iΔh0xV¯iαxV¯ir0fγ5xV¯iΔγ0xV¯i2.Similar to (76), the first term of (78) can be rewritten as
(79)Δγ11V¯=γα11V¯α+γα2V¯αsqu,where
(80)γα11V¯=WTV¯Fγ2AcbFV3+Fγ3AcbFγ1,γα2V¯=WTV¯Fγ2AcbFV1..Substituting (69) and (76) into the second term of (78), it can be expressed as
(81)Δγ21V¯=12V¯−V¯01r0∑i=0NWixfγ4xV¯ihα10xV¯iα+hα20xV¯iαsquαxV¯i+12V¯−V¯0∑i=0NWixfγ5xV¯iγα10xV¯iα2=αTγα3V¯α+αsquTγα4V¯α,where γα3V¯ and γα4V¯ are the coefficient matric expressed as
(82)γα3V¯=12V¯−V¯01r0∑i=0NWiV¯fγ4xV¯ihα10xV¯iTIi+12V¯−V¯0∑i=0NWiV¯fγ5xV¯iγα10xV¯iTγα10xV¯i,γα4V¯=12V¯−V¯01r0∑i=0NWiV¯fγ4xV¯ihα20xV¯iTIi,where Ii is the i-th row of the identity matrix with the dimension N+1. Then, combining (79) with (81), we get
(83)Δγ1V¯=γα11V¯α+γα2V¯αsqu+αTγα3V¯α+αsquTγα4V¯α.By using Chebyshev polynomials to approximate the integrands, the analytical solution toΔV1 can be obtained as
(84)ΔV1V¯=12V¯−V¯0Ve∑i=0NWixfV2xV¯iΔV0xV¯iVe+fV4xV¯iΔh0xV¯ir0+12V¯−V¯0Ve∑i=0NWixfV5xV¯iΔh0xV¯iα2xV¯ir0+fV6xV¯iΔγ0xV¯i2.Similar to (76), the first term of (84) can be rewritten as
(85)ΔV11V¯=Vα10V¯α+Vα20V¯αsqu,where
(86)Vα11V¯=VeWTV¯FV2AcbFV3+FV4AcbFh1AcbFV3+Fh2AcbFγ1,Vα21V¯=VeWTV¯FV2+FV4AcbFh1AcbFV1.Similarly, substituting (69) and (76) into the second term of (84), it can be expressed as
(87)ΔV21V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯ihα10xV¯iα+hα20xV¯iαsquα2xV¯i+12V¯−V¯0Ve∑i=0NWixfV6xV¯iγα10xV¯iα2=αTVα3V¯α+αsquTVα4V¯α+αsquTVα5V¯αsqu,where
(88)Vα3V¯=12V¯−V¯0Ve∑i=0NWixfV6xV¯iγα10xV¯iTγα10xV¯i,Vα4V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯iITihα10xV¯i,Vα5V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯iITihα20xV¯i.Then, combining (85) with (87), the analytical solution to ΔV1 is
(89)ΔV1V¯=Vα10V¯α+Vα20V¯αsqu+αTVα3V¯α+αsquTVα4V¯α+αsquTVα5V¯αsqu.Similar to (78) and (84), the analytical solution to Δh1 can be approximated as
(90)Δh1V¯=12V¯−V¯0r0∑i=0NWixfh1xV¯iΔV1xV¯iVe+r0∑i=0NWiV¯fh2xV¯iΔγ1xV¯i+12V¯−V¯0r0∑i=0NWixfh3xV¯iΔV0xV¯iΔγ0xV¯iVe+fh4xV¯iΔγ0xV¯i2.Substituting (89) into (90), the first term of (90) can be rewritten as
(91)Δh11V¯=12V¯−V¯0r0Ve∑i=0NWixfh1xV¯iVα11xV¯iα+Vα21xV¯iαsqu+αTVα3xV¯iα+αsquTVα4xV¯iα+αsquTVα5xV¯iαsqu=hα1,11V¯α+hα2,11V¯αsqu+αThα3,1V¯α+αsquThα4,1V¯α+αsquThα5V¯α,where
(92)hα1,11V¯=r0WTV¯Fh1AcbFV2AcbFV3+FV4AcbFh1AcbFV3+Fh2AcbFγ1,hα2,11V¯=r0WTV¯Fh1AcbFV2+FV4AcbFh1AcbFV1,hα3,1V¯=12V¯−V¯0r0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV6xV¯jγα10xV¯jTγα10xV¯j,hα4,1V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV5xV¯jITjhα10xV¯j,hα5V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV5xV¯jITjhα20xV¯j.Similar to (91), the second term of (90) can be rewritten as
(93)Δh21V¯=hα1,21V¯α+hα2,21V¯αsqu+αThα3,2V¯α+αsquThα4,2V¯α,where
(94)hα1,21V¯=r0WTV¯Fh2AcbFγ2AcbFV3+Fγ3AcbFγ1,hα2,21V¯=r0WTV¯Fh2AcbFγ2AcbFV1,hα3,2V¯=12V¯−V¯0r0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ5xV¯jγα10xV¯jTγα10xV¯j+12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ4xV¯jhα10xV¯jTIj,hα4,2V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ4xV¯jhα20xV¯jTIj.Substituting (69) and (72) into the third term of (90), it can be expressed as
(95)Δh31V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα10xV¯iα+Vα20xV¯iαsquγα10xV¯iα+12V¯−V¯0r0∑i=0NWixfh4xV¯iγα10xV¯iα2=αThα3,3V¯α+αsquThα4,3V¯α,where
(96)hα3,3V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα10xV¯iTγα10xV¯i+12V¯−V¯0r0∑i=0NWixfh4xV¯iγα10xV¯iTγα10xV¯i,hα4,3V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα20xV¯iTγα10xV¯i.Therefore, combining (91), (93), and (95), the analytical solution to Δh1 is expressed as
(97)Δh1V¯=hα11V¯α+hα21V¯αsqu+αThα3V¯α+αsquThα4V¯α+αsquThα5V¯α,where
(98)hα11V¯=hα1,11V¯+hα1,21V¯,hα21V¯=hα2,11V¯+hα2,21V¯,hα3V¯=hα3,1V¯+hα3,2V¯+hα3,3V¯,hα4V¯=hα4,1V¯+hα4,2V¯+hα4,3V¯..


#### 3.3.4. Full Analytical Solution to Trajectory with Arbitrary AOA

When the zero-order and first-order solutions are derived, the full analytical solutions to trajectory increment can be formulated as(99)ΔxV¯=Xα1V¯α+Xα2V¯αsqu+ATXα3V¯α+AsquTXα4V¯α+AsquTXα5V¯αsqu,where A and Asqu are the partitioned diagonal matrices expressed as (100), respectively.
(100)A=diagα,α,α,Asqu=diagαsqu,αsqu,αsqu.The coefficients are(101)XαiV¯=VαiTV¯,γαiTV¯,hαiTV¯T,i=1,2,3,4,5,and each coefficient is determined by the sum of the zeroth-order and first-order terms.Then, the analytical solution to ascent trajectory can be expressed as(102)x=xbV¯+Xα1V¯α+Xα2V¯αsqu+ATXα3V¯α+AsquTXα4V¯α+AsquTXα5V¯αsqu,where xbV¯ is the analytical solution to the trajectory with zero AOA, which is expressed as
(103)xbV¯=Vb0+Vb1,γb0+γb1,hb0+hb1T.It should be noted that the analytical solution includes the high order of AOA which leads to higher accuracy for states in all region than that only first order of AOA is considered.


## 3.1. Dynamics for Analytical Trajectory

In order to simplify the solutions, the index atmospheric density model,ρ=ρ0exp−βh, where ρ0=1.225kg/m3 is the sea-level atmospheric density and β=1.389×10−4m−1, is adopted for the solving of the analytical solutions.In order to avoid confusion, the variables corresponding to zero AOA are expressed with subscript “b.” For the case of zero AOA, it is obvious that P¯sinα+L=0 and D=ρbVb2SC0/2; then, the dynamics is simplified to
(8)V̇b=1mP¯−ρbVb2SC¯d02−gbsinγb,γ̇b=Vbrb−gbVbcosγb,ḣb=Vbsinγb.For simplifying the derivation, a variable related toγb is introduced.
(9)θb=ln1+sinγb1−sinγb.Replacingγb with θb as the state variable, the equation of motion presented in (8) can be rewritten as
(10)V̇b=1mP¯−ρbVb2SC¯d02−gbsinγb,θ̇b=2Vbrb−gbVb,ḣb=Vbsinγb.The increment of the velocity, flight-path angle, and altitude with nonzero AOA relative to that with zero AOA are defined asΔV, Δγ, and Δh, which are used to derive the analytical solutions with arbitrary AOA profile. It is obvious that the differential equations of those trajectory increments are the difference of (1) and (8). For the sake of simplification, the following Taylor’s expansions are adopted.
(11)exp−βΔh≈1−βΔh,sinα≈α,1V≈1Vb−ΔVVb2,cosα≈1−α22,sinγ≈sinγb+cosγbΔγ−sinγbΔγ22,cosγ≈cosγb−sinγbΔγ−cosγbΔγ22.Then, the dynamics equations of increments are established as(12)ΔV̇=−12mP¯+ρbVb2SC¯dαα2+2ρbVbSCdΔV−βρbVb2SCdΔh−gicosγbΔγ+12gisinγbΔγ2,Δγ̇=1mP¯Vb+12ρbVbSC¯lα1−βΔhα+1r0+giVb2cosγbΔV−Vbr0−giVbsinγb+0.5cosγbΔγΔγ,Δḣ=sinγb+cosγbΔγΔV+VbcosγbΔγ−sinγbΔγ22.Becauser≫Δh, the altitude increment on the value of 1/r is ignored. It should be noted that the second order terms related to the AOA and state increment are remained to achieve high accuracy.


## 3.2. Analytical Ascent Trajectory with Zero AOA

In (10), the right side of the differential equation of Vb can be divided into the following two parts:
(13)V̇b=Favem+Fεm,where Fave is the major axial force, which is the sum of thrust, average atmospheric drag, and average axial component of gravity, and Fε is the perturbation axial force:
(14)Fave=P¯−ρaveVave2SC¯d02−mavegavesinγave,Fε=ρaveVave2SC¯d02+mavegavesinγave−ρbVb2SC¯d02−mgbsinγb,where mave,γave, Vave, ρave, and gave are the average mass, flight-path angle, velocity, atmospheric density, and gravitational acceleration, respectively, which can be previously determined by the average of its initial and terminal values. In the boost phase, the thrust is usually much greater than the change of aerodynamic drag and the axial component of gravity. Thus, the perturbation axial force is small enough and can be regarded as the correction terms. On the other hand, the trigonometric term sinγb can be written as follows:
(15)sinγb=sinγ0+c1Δθb+sinγb−sinγ0−c1Δθb,Δθb=θb−θ0,where c1 is the constant coefficient, which can be obtained by the initial and terminal values of γb and θb:
(16)c1=sinγbf−sinγ0θbf−θ0.In the boost phase, the flight-path angle varies slightly; therefore, the linear termsinγ0+c1Δθb is much larger than the rest terms. Therefore, according to the regular perturbation method presented in Appendix A, (10) can be modified as
(17)V̇bθ̇bḣb=Favem2Vbrb−gbVbVbsinγ0+c1Δθb+εkFεm0Vbsinγb−sinγ0−c1Δθb,where ε is a small parameter for the subsequent regular perturbation order and k is a constant, which is set to be equal to ε. According to (A.5), the zeroth-order dynamics is
(18)V̇b0=Favem,θ̇b0=2Vb0r0−giVb0,ḣb0=Vb0sinγ0+c1Δθb0.It is noted that, in order to decouple the differential equations ofθb0 and hb0, the terms related to r are expanded at initial point since the change of r is much smaller than its initial value, r0. The first-order dynamics is
(19)V̇b1=1mP¯−Fave−12ρ0exp−βhb0Vb02SC¯d0−gavesinγ0+c1Δθb0,θ̇b1=−2Vb0r0−2giVb0Δhb0r0+2Vb0r0+giVb0Vb1Vb0,ḣb1=Vb0sinγb0−sinγ0−c1Δθb0+Vb1sinγb0+Vb0cos2γb0θb12,where γb0 is the zeroth-order solution to flight-path angle, which can be derived from θb0.
(20)γb0=sin−1expθb0−1expθb0+1.It should be noted that the terms related tohb1 in the differential equations of θb1 are ignored due to hb1≪Δhb0. Additionally, according to the boundary conditions, the initial values of zeroth-order and first-order states are
(21)Vb00=V0,θb00=θ0,hb00=h0,Vb01=0,θb01=0,hb01=0.


### 3.2.1. Zeroth-Order Solution

It is obvious that the three differential equations in (18) are not coupled, so that the analytical solution to Vb0, θb0, and hb0 can be derived by integrating these three differential equations in sequence.Integrating the first equation of (18) from t0 to t, the analytical solution of Vb0 at time t can be easily obtained as
(22)Vb0=V0−Veln1−t−t0T,where T=m0/qm and Ve=Fave/qm. Obviously, the differential equations of other variables are all related to Vb0; thus, the analytical solutions of these variables are the functions of Vb0. For the sake of simplifying the expression, a dimensionless variable, V¯=Vb0/Ve, is employed, which is regarded as an independent variable in the subsequent derivation. From (22), the expression of V¯ is
(23)V¯=V¯0−ln1−t−t0T,where V¯0=V0/Ve is initial value of V¯. Substituting (22) into the first equation of (18), the differential equation of V¯ can be expressed as a simple function.
(24)V¯̇=FavemVe=expV¯−V¯0T.Dividing the second and third equations of (18) by (24) yields
(25)dθb0dV¯=2K1V¯−K2V¯exp−V¯,(26)dhb0dV¯=K1r0sinγ0V¯+c1V¯Δθb0exp−V¯,where K1=TVeexpV¯0/r0and K2=giTexpV¯0/Ve. Integrating (25) from V¯0 to V¯, θb0 can be obtained as
(27)θb0V¯=θ0+2K1fV¯V¯;1,1−2K2ξV¯V¯;1,1,where fV¯V¯;n,mn≥0,m≥0 is
(28)fV¯V¯;n,m=V¯n+1−V¯0n+1n+1,m=0,n!mn+1∑k=0nmkk!V¯0kexp−mV¯0−V¯kexp−mV¯,m>0.ξV¯V¯;n,mn>0,m>0 is the integral term expressed as
(29)ξV¯V¯;n,m=∫V¯0V¯exp−mxV¯xV¯ndxV¯.Although the integral term is quite simple, it has no analytical solution. In order to derive its approximate analytical solution, Chebyshev interpolation polynomial is used to approximate the integrand. Let(30)exp−mxV¯xV¯n≈∑i=0Nexp−mxV¯ixV¯inLix=∑i=0NpV¯in,mxi,where xV¯i is the value of V¯ at the i-th interpolation point xi; Lixi=0,⋯,N is the interpolation basis function; and x=xV¯−V¯2/V¯1 is the dimensionless variable with −1,1, where V¯1=V¯−V¯0/2 and V¯2=V¯+V¯0/2. Then, ξV¯ can be approximated as
(31)ξV¯V¯;n,m=12V¯−V¯0∑i=0NWixexp−mxV¯ixV¯in,where Wix is the integral weight coefficient. The derivation process is presented in Appendix B in detail. The error analysis is carried out by comparing the solution obtained by numerical integration and the 6th-degree interpolation polynomials in Figure 1, from which we can see that the interpolation polynomials have high accuracy.Figure 1
Comparison ofξV¯V¯ with various n and m.In order to derive the analytical solution tohb0, an integral related to Δθb0 is introduced.
(32)IθxV¯;n,m=∫V¯0xV¯Δθb0xV¯nexp−mxV¯dxV¯,n≥0,m>0.Using the formula of partial integration, the analytical solution toIθxV¯;n,m is expressed as
(33)IθxV¯;n,m=−n!mn+1∑k=0nmkk!Δθb0xV¯xV¯kexp−mxV¯−θrxV¯;k,m,where θrxV¯;n,m is the integral shown below:
(34)θrxV¯;n,m=∫θ0θb0xV¯nexp−mxV¯dθb0;n≥0,m>0.It is obvious that the analytical solution toθr is similar to that of θb0. Substituting (25) into (34), θr can be expressed as
(35)θrxV¯;n,m=2∫V¯0xV¯K1xV¯n+1exp−m+1xV¯−K2xV¯n−1exp−m+1xV¯dxV¯=2K1fV¯xV¯;1,m+1−2K2ξV¯xV¯;1,m+1,n=0,2K1fV¯xV¯;n+1,m+1−2K2fV¯xV¯;n−1,m+1,n>0.Then, integrating (26) from V¯0 to V¯, the analytical solution to hb0 is
(36)hb0V¯=h0+K1r0sinγ0fV¯V¯;1,1+c1IθV¯;1,1,where
(37)IθV¯;1,1=2K1fV¯V¯;1,2+fV¯V¯;2,2−1+V¯exp−V¯fV¯V¯;1,1−2K2fV¯V¯;0,2+ξV¯V¯;1,2−1+V¯exp−V¯ξV¯V¯;1,1.Finally, the zeroth-order solutions with zero AOA are derived, which are the fundamental for deriving the first-order correction.


### 3.2.2. First-Order Correction

In this subsection, let us focus on the first-order dynamics. Regarding the dimensionless velocity,V¯, as the dependent variable, the first-order dynamics can be formulated as
(38)dVb1dV¯=VeKP−1−KaC¯d0V¯2exp−β¯Δh¯b0−K2sinγ0+c1Δθb0exp−V¯,dθb1dV¯=−2K1V¯−2K2V¯exp−V¯Δh¯b0+2K1+K2V¯2exp−V¯V¯b1,dhb1dV¯=K1r0exp−V¯V¯sinγb0−c1Δθb0−sinγ0+V¯b1sinγb0+V¯cos2γb0θb12,where Δh¯b0=hb0−h0/r0 and V¯b1=Vb1/Ve are the dimensionless variables and β¯=βr0,KP=P¯/Fave, and Ka=ρiSVe2/2Fave are the constant coefficients.The analytical solution toVb1 can be derived by integrating the first equation of (38). But the term, which is the integration of V¯2exp−β¯Δh¯b0, cannot be solved analytically because Δh¯b0 is not the linear function of V¯. Similar to (30), let
(39)fρxV¯=xV¯2exp−β¯Δh¯b0≈∑i=0NfρxV¯iLix,where fρxV¯i=xV¯i2exp−β¯K1sinγ0fV¯xV¯i;1,1+c1IθxV¯i;1,1. As shown in Figure 2, the approximation polynomials are accurate enough in comparison with the original functions fρxV¯. Then, the analytical solution to Vb1 is
(40)Vb1V¯=VeKP−1V¯−V¯0−12V¯−V¯0KaC¯d0∑i=0NfρxV¯iWix−K2sinγ0exp−V¯0−exp−V¯−K2c1IθV¯;0,1,where according to (33), the expression of IθV¯;0,1 is
(41)IθV¯;0,1=2K1fV¯V¯;1,2−exp−V¯fV¯V¯;1,1−2K2ξV¯V¯;1,2−exp−V¯ξV¯V¯;1,1.Figure 2
Comparison betweenfρV¯ and the corresponding approximate polynomial.By integrating the second equation of (38), the analytical solution to θb1 can be obtained as
(42)θb1V¯=−2K1h¯r00+h¯r10−Δh¯b01+V¯exp−V¯+2K1V¯r1−V¯b1exp−V¯+4K2∫V0V¯Δh¯b0exp−xV¯xV¯dxV¯+2K2∫V0V¯V¯b1exp−xV¯xV¯2dxV¯,where h¯r00, h¯r10, and V¯r1 are all the functions of V¯, which are expressed as
(43)h¯r00V¯=K1sinγ0fV¯V¯;1,2+c1IθV¯;1,2,h¯r10V¯=K1sinγ0fV¯V¯;2,2+c1IθV¯;2,2,V¯r1V¯=KP−1fV¯V¯;0,1−12V¯−V¯0KaC¯d0∑i=0NfρxV¯iexp−xV¯iWix−K2sinγ0fV¯V¯;0,2−K2c1IθV¯;0,2,where the expression of IθV¯;0,2, IθV¯;1,2, and IθV¯;1,2 can be obtained from (33).As for the two integration terms in (42), there are no analytical solutions. Chebyshev interpolation polynomials are also adopted to approximate the integrands. Let
(44)exp−mxV¯xV¯n=∑i=0NcV¯in,mxV¯i,where the coefficient cV¯in,m is
(45)cV¯in,m=∑k=iNpV¯kn,mCnk−V¯2/V¯1n−kV¯1k.Then, the analytical solutions to the two integration terms are(46)∫V0V¯Δh¯b0exp−xV¯xV¯dxV¯=∑i=0NcV¯i1,1IhV¯;i,∫V0V¯V¯b1exp−xV¯xV¯2dxV¯=∑i=0NcV¯i2,1IV¯V¯;i,where
(47)IhV¯;n=∫V¯0V¯Δh¯b0xV¯nexp−xV¯dxV¯,IV¯V¯;n=∫V¯0V¯V¯b1xV¯nexp−xV¯dxV¯.By the formula of partial integration, the analytical solution toIhV¯;n and IV¯V¯;n can be derived as
(48)IhV¯;n=−n!∑k=0n1k!Δh¯b0V¯kexp−V¯−ϑhV¯;k,IV¯V¯;n=−n!∑k=0n1k!V¯b1V¯kexp−V¯−ϑV¯V¯;k,where the expressions of ϑhV¯;n and ϑV¯V¯;n are
(49)ϑhV¯;n=K1sinγ0fV¯V¯;n+1,2+c1IθV¯;n+1,2,ϑV¯V¯;n=KP−1fV¯V¯;n,1−12V¯−V¯0KaCd0∑i=0NfρxV¯ixV¯inexp−xV¯iWiV¯−K2sinγ0fV¯V¯;n,2−K2c1IθV¯;n,2.Then, the analytical solution toθb1 is
(50)θb1V¯=−2K1h¯r00+h¯r10−Δh¯b01+V¯exp−V¯+4K2∑i=0NcV¯i1,1IhV¯;i+2K1V¯r1−V¯b1exp−V¯+2K2∑i=0NcV¯i2,1IV¯V¯;i.As for the analytical solution tohb1, it is the integral of dhb0/dV¯ from V¯0 to V¯. It is noted that all the variables in dhb0/dV¯, such as γb0, θb0, V¯b1, and θb1, are expressed as the function of V¯. Therefore, we can use Chebyshev interpolation polynomial to approximate the integrand. Let
(51)fhxV¯=exp−xV¯xV¯sinγb0−c1Δθb0−sinγ0+V¯b1sinγb0+xV¯cos2γb0θb12≈∑i=0NfhxV¯iLix.Then, the solution tohb1 can be obtained as
(52)hb1V¯=K1r0∑i=0NfhxV¯iWiV¯.Finally, the first-order solutions are derived. According to (A.6), the analytical solution with zero AOA is the sum of the zero-order and first-order solutions.


## 3.2.1. Zeroth-Order Solution

It is obvious that the three differential equations in (18) are not coupled, so that the analytical solution to Vb0, θb0, and hb0 can be derived by integrating these three differential equations in sequence.Integrating the first equation of (18) from t0 to t, the analytical solution of Vb0 at time t can be easily obtained as
(22)Vb0=V0−Veln1−t−t0T,where T=m0/qm and Ve=Fave/qm. Obviously, the differential equations of other variables are all related to Vb0; thus, the analytical solutions of these variables are the functions of Vb0. For the sake of simplifying the expression, a dimensionless variable, V¯=Vb0/Ve, is employed, which is regarded as an independent variable in the subsequent derivation. From (22), the expression of V¯ is
(23)V¯=V¯0−ln1−t−t0T,where V¯0=V0/Ve is initial value of V¯. Substituting (22) into the first equation of (18), the differential equation of V¯ can be expressed as a simple function.
(24)V¯̇=FavemVe=expV¯−V¯0T.Dividing the second and third equations of (18) by (24) yields
(25)dθb0dV¯=2K1V¯−K2V¯exp−V¯,(26)dhb0dV¯=K1r0sinγ0V¯+c1V¯Δθb0exp−V¯,where K1=TVeexpV¯0/r0and K2=giTexpV¯0/Ve. Integrating (25) from V¯0 to V¯, θb0 can be obtained as
(27)θb0V¯=θ0+2K1fV¯V¯;1,1−2K2ξV¯V¯;1,1,where fV¯V¯;n,mn≥0,m≥0 is
(28)fV¯V¯;n,m=V¯n+1−V¯0n+1n+1,m=0,n!mn+1∑k=0nmkk!V¯0kexp−mV¯0−V¯kexp−mV¯,m>0.ξV¯V¯;n,mn>0,m>0 is the integral term expressed as
(29)ξV¯V¯;n,m=∫V¯0V¯exp−mxV¯xV¯ndxV¯.Although the integral term is quite simple, it has no analytical solution. In order to derive its approximate analytical solution, Chebyshev interpolation polynomial is used to approximate the integrand. Let(30)exp−mxV¯xV¯n≈∑i=0Nexp−mxV¯ixV¯inLix=∑i=0NpV¯in,mxi,where xV¯i is the value of V¯ at the i-th interpolation point xi; Lixi=0,⋯,N is the interpolation basis function; and x=xV¯−V¯2/V¯1 is the dimensionless variable with −1,1, where V¯1=V¯−V¯0/2 and V¯2=V¯+V¯0/2. Then, ξV¯ can be approximated as
(31)ξV¯V¯;n,m=12V¯−V¯0∑i=0NWixexp−mxV¯ixV¯in,where Wix is the integral weight coefficient. The derivation process is presented in Appendix B in detail. The error analysis is carried out by comparing the solution obtained by numerical integration and the 6th-degree interpolation polynomials in Figure 1, from which we can see that the interpolation polynomials have high accuracy.Figure 1
Comparison ofξV¯V¯ with various n and m.In order to derive the analytical solution tohb0, an integral related to Δθb0 is introduced.
(32)IθxV¯;n,m=∫V¯0xV¯Δθb0xV¯nexp−mxV¯dxV¯,n≥0,m>0.Using the formula of partial integration, the analytical solution toIθxV¯;n,m is expressed as
(33)IθxV¯;n,m=−n!mn+1∑k=0nmkk!Δθb0xV¯xV¯kexp−mxV¯−θrxV¯;k,m,where θrxV¯;n,m is the integral shown below:
(34)θrxV¯;n,m=∫θ0θb0xV¯nexp−mxV¯dθb0;n≥0,m>0.It is obvious that the analytical solution toθr is similar to that of θb0. Substituting (25) into (34), θr can be expressed as
(35)θrxV¯;n,m=2∫V¯0xV¯K1xV¯n+1exp−m+1xV¯−K2xV¯n−1exp−m+1xV¯dxV¯=2K1fV¯xV¯;1,m+1−2K2ξV¯xV¯;1,m+1,n=0,2K1fV¯xV¯;n+1,m+1−2K2fV¯xV¯;n−1,m+1,n>0.Then, integrating (26) from V¯0 to V¯, the analytical solution to hb0 is
(36)hb0V¯=h0+K1r0sinγ0fV¯V¯;1,1+c1IθV¯;1,1,where
(37)IθV¯;1,1=2K1fV¯V¯;1,2+fV¯V¯;2,2−1+V¯exp−V¯fV¯V¯;1,1−2K2fV¯V¯;0,2+ξV¯V¯;1,2−1+V¯exp−V¯ξV¯V¯;1,1.Finally, the zeroth-order solutions with zero AOA are derived, which are the fundamental for deriving the first-order correction.


## 3.2.2. First-Order Correction

In this subsection, let us focus on the first-order dynamics. Regarding the dimensionless velocity,V¯, as the dependent variable, the first-order dynamics can be formulated as
(38)dVb1dV¯=VeKP−1−KaC¯d0V¯2exp−β¯Δh¯b0−K2sinγ0+c1Δθb0exp−V¯,dθb1dV¯=−2K1V¯−2K2V¯exp−V¯Δh¯b0+2K1+K2V¯2exp−V¯V¯b1,dhb1dV¯=K1r0exp−V¯V¯sinγb0−c1Δθb0−sinγ0+V¯b1sinγb0+V¯cos2γb0θb12,where Δh¯b0=hb0−h0/r0 and V¯b1=Vb1/Ve are the dimensionless variables and β¯=βr0,KP=P¯/Fave, and Ka=ρiSVe2/2Fave are the constant coefficients.The analytical solution toVb1 can be derived by integrating the first equation of (38). But the term, which is the integration of V¯2exp−β¯Δh¯b0, cannot be solved analytically because Δh¯b0 is not the linear function of V¯. Similar to (30), let
(39)fρxV¯=xV¯2exp−β¯Δh¯b0≈∑i=0NfρxV¯iLix,where fρxV¯i=xV¯i2exp−β¯K1sinγ0fV¯xV¯i;1,1+c1IθxV¯i;1,1. As shown in Figure 2, the approximation polynomials are accurate enough in comparison with the original functions fρxV¯. Then, the analytical solution to Vb1 is
(40)Vb1V¯=VeKP−1V¯−V¯0−12V¯−V¯0KaC¯d0∑i=0NfρxV¯iWix−K2sinγ0exp−V¯0−exp−V¯−K2c1IθV¯;0,1,where according to (33), the expression of IθV¯;0,1 is
(41)IθV¯;0,1=2K1fV¯V¯;1,2−exp−V¯fV¯V¯;1,1−2K2ξV¯V¯;1,2−exp−V¯ξV¯V¯;1,1.Figure 2
Comparison betweenfρV¯ and the corresponding approximate polynomial.By integrating the second equation of (38), the analytical solution to θb1 can be obtained as
(42)θb1V¯=−2K1h¯r00+h¯r10−Δh¯b01+V¯exp−V¯+2K1V¯r1−V¯b1exp−V¯+4K2∫V0V¯Δh¯b0exp−xV¯xV¯dxV¯+2K2∫V0V¯V¯b1exp−xV¯xV¯2dxV¯,where h¯r00, h¯r10, and V¯r1 are all the functions of V¯, which are expressed as
(43)h¯r00V¯=K1sinγ0fV¯V¯;1,2+c1IθV¯;1,2,h¯r10V¯=K1sinγ0fV¯V¯;2,2+c1IθV¯;2,2,V¯r1V¯=KP−1fV¯V¯;0,1−12V¯−V¯0KaC¯d0∑i=0NfρxV¯iexp−xV¯iWix−K2sinγ0fV¯V¯;0,2−K2c1IθV¯;0,2,where the expression of IθV¯;0,2, IθV¯;1,2, and IθV¯;1,2 can be obtained from (33).As for the two integration terms in (42), there are no analytical solutions. Chebyshev interpolation polynomials are also adopted to approximate the integrands. Let
(44)exp−mxV¯xV¯n=∑i=0NcV¯in,mxV¯i,where the coefficient cV¯in,m is
(45)cV¯in,m=∑k=iNpV¯kn,mCnk−V¯2/V¯1n−kV¯1k.Then, the analytical solutions to the two integration terms are(46)∫V0V¯Δh¯b0exp−xV¯xV¯dxV¯=∑i=0NcV¯i1,1IhV¯;i,∫V0V¯V¯b1exp−xV¯xV¯2dxV¯=∑i=0NcV¯i2,1IV¯V¯;i,where
(47)IhV¯;n=∫V¯0V¯Δh¯b0xV¯nexp−xV¯dxV¯,IV¯V¯;n=∫V¯0V¯V¯b1xV¯nexp−xV¯dxV¯.By the formula of partial integration, the analytical solution toIhV¯;n and IV¯V¯;n can be derived as
(48)IhV¯;n=−n!∑k=0n1k!Δh¯b0V¯kexp−V¯−ϑhV¯;k,IV¯V¯;n=−n!∑k=0n1k!V¯b1V¯kexp−V¯−ϑV¯V¯;k,where the expressions of ϑhV¯;n and ϑV¯V¯;n are
(49)ϑhV¯;n=K1sinγ0fV¯V¯;n+1,2+c1IθV¯;n+1,2,ϑV¯V¯;n=KP−1fV¯V¯;n,1−12V¯−V¯0KaCd0∑i=0NfρxV¯ixV¯inexp−xV¯iWiV¯−K2sinγ0fV¯V¯;n,2−K2c1IθV¯;n,2.Then, the analytical solution toθb1 is
(50)θb1V¯=−2K1h¯r00+h¯r10−Δh¯b01+V¯exp−V¯+4K2∑i=0NcV¯i1,1IhV¯;i+2K1V¯r1−V¯b1exp−V¯+2K2∑i=0NcV¯i2,1IV¯V¯;i.As for the analytical solution tohb1, it is the integral of dhb0/dV¯ from V¯0 to V¯. It is noted that all the variables in dhb0/dV¯, such as γb0, θb0, V¯b1, and θb1, are expressed as the function of V¯. Therefore, we can use Chebyshev interpolation polynomial to approximate the integrand. Let
(51)fhxV¯=exp−xV¯xV¯sinγb0−c1Δθb0−sinγ0+V¯b1sinγb0+xV¯cos2γb0θb12≈∑i=0NfhxV¯iLix.Then, the solution tohb1 can be obtained as
(52)hb1V¯=K1r0∑i=0NfhxV¯iWiV¯.Finally, the first-order solutions are derived. According to (A.6), the analytical solution with zero AOA is the sum of the zero-order and first-order solutions.


## 3.3. Analytic Ascent Trajectory with Arbitrary AOA Profile

In the last subsection, the analytical solutions to ascent trajectory with zero AOA are provided. In order to obtain the full analytical solutions to ascent trajectory with arbitrary AOA profile, the analytical solutions to trajectory increment corresponding to no-zero AOA are presented firstly in this subsection. Then, the full analytical solutions are the sum of analytical trajectory with zero AOA and trajectory increment corresponding to nonzero AOA.


### 3.3.1. Regular Perturbation Model with AOA Increment

Because the independent variable for the analytical solution with zero AOA is the dimensionless velocity,V¯, the differential equations of trajectory increments should be formulated as
(53)dΔVdV¯=VefV1α2+fV2ΔVVe+fV3Δγ+fV4Δhr0+fV5Δhα2r0+fV6Δγ2,dΔγdV¯=fγ1α+fγ2ΔVVe+fγ3Δγ+fγ4Δhαr0+fγ5Δγ2,dΔh¯dV¯=r0fh1ΔVVe+fh2Δγ+fh3ΔVΔγVe+fh4Δγ2,,where fV1∼fV6, fγ1∼fγ5, and fh1∼fh4 are all the functions of V¯, which are expressed as
(54)fV1V¯=−KP2+KaC¯dαV¯b2exp−βΔhb,fV2V¯=−2KaC¯d0V¯bexp−βΔhb,fV3V¯=−K2cosγbexp−V¯,fV4V¯=βr0KaCd0V¯b2exp−βΔhb,fV5V¯=βr0KaC¯dαV¯b2exp−βΔhb,fV6V¯=K2sinγbexp−V¯2,fγ1V¯=KPV¯b+KaC¯lαV¯bexp−βΔhb,fγ2V¯=K1+K2V¯b2cosγbexp−V¯,fγ3V¯=−K1V¯b−K2V¯bsinγbexp−V¯,fγ4V¯=−βr0KaC¯lαV¯bexp−βΔhb,fγ5V¯=−K1V¯b−K2/V¯bcosγbexp−V¯2,fh1V¯=K1sinγbexp−V¯,fh2V¯=K1V¯bcosγbexp−V¯,fh3V¯=K1cosγbexp−V¯,fh4V¯=−K1V¯bsinγbexp−V¯2..Generally, the AOA profile is expressed as the function of the flight time. The AOA profile with the dimensionless velocity can be formulated by using the relationship between the flight timet and V¯, which can be easily expressed as
(55)t=t0+T1−expV¯0−V¯.In the differential equation ofΔγ, it is easy to find that the term, fγ1α, is much larger than the sum of other terms, fγ1α≫fγ2ΔV/Ve+fγ3Δγ+fγ4Δhα/r0+fγ5Δγ2, which is the dominant term. Similarly, the term, fV1α2+fV3Δγ, is the dominant term of the differential equation of ΔV. As for the differential equation of Δh, due to Vb≫ΔV and Δγ2≤1, it is obvious that the term, fh1ΔV/Ve+fh2Δγ, is the dominant term. Therefore, the dynamics can be modified as
(56)dΔγdV¯dΔVdV¯dΔhdV¯=fγ1αVefV1α2+VefV3Δγr0Vefh1ΔV+r0fh2Δγ+εk1Vefγ2ΔV+fγ3Δγ+1r0fγ4Δhα+fγ5Δγ2fV2ΔV+Ver0fV4Δh+Ver0fV5Δhα2+VefV6Δγ2r0Vefh3ΔVΔγ+r0fh4Δγ2.According to the regular perturbation method presented in AppendixB, the zeroth-order and first-order dynamics are
(57)dΔγ0dV¯=fγ1V¯α,dΔV0dV¯=VefV1V¯α2+fV3V¯Δγ0,dΔh0dV¯=r0fh1V¯ΔV0Ve+fh2V¯Δγ0,,(58)dΔγ1dV¯=fγ2V¯ΔV0Ve+fγ3V¯Δγ0+fγ4V¯Δh0αr0+fγ5V¯Δγ02,dΔV1dV¯=VefV2V¯ΔV0Ve+fV4V¯Δh0r0+fV5V¯Δh0α2r0+fV6V¯Δγ02,dΔh1dV¯=r0fh1V¯ΔV1Ve+fh2V¯Δγ1+fh3V¯ΔV0Δγ0Ve+fh3V¯Δγ02..The initial values of zeroth-order and first-order states are all zeros.


### 3.3.2. Zeroth-Order Solution to Trajectory Increment

It is obvious that the three differential equations of zero-order states are not coupled. Therefore, by integrating these differential equations orderly, the analytical solutions toΔγ0, ΔV0, and Δh0 can be obtained as
(59)Δγ0V¯=∫V¯0V¯fγ1xV¯αxV¯dxV¯,(60)ΔV0V¯=Ve∫V¯0V¯fV1xV¯α2xV¯dxV¯+Ve∫V¯0V¯fV3xV¯Δγ0xV¯dxV¯,(61)Δh0V¯=r0Ve∫V¯0V¯fh1xV¯ΔV0xV¯dxV¯+r0∫V¯0V¯fh2xV¯Δγ0xV¯dxV¯.Although the integral terms in the above formulas are quite complex, the integrands are related only to the independent variable,V¯. Here, Chebyshev polynomials are also used to approximate the integrands. Let
(62)fγ1xV¯αxV¯≈∑i=0Nfγ1xV¯iαxV¯iLix,(63)fV1xV¯α2xV¯≈∑i=0NfV1xV¯iα2xV¯iLix,(64)fV3xV¯Δγ0xV¯≈∑i=0NfV3xV¯iΔγ0xV¯iLix,(65)fh1xV¯ΔV¯0xV¯≈∑i=0Nfh1xV¯iΔV¯0xV¯iLix,(66)fh2xV¯Δγ0xV¯≈∑i=0Nfh2xV¯iΔγ0xV¯iLix.Substituting (62) into (59), the analytical solution to the zeroth-order flight-path angle increment can be obtained as
(67)Δγ0V¯=12V¯−V¯0∑i=0Nfγ1xV¯iαxV¯iWix=γα10V¯α,where α is the vector consisting of the value of AOA at interpolating points, i.e.,
(68)α=αxV¯0,αxV¯1,⋯,αxV¯NT,γα10V¯ is the coefficient matrix expressed as
(69)γα10V¯=12V¯−V¯0W0V¯fγ1xV¯0,⋯,WNV¯fγ1xV¯N=WTV¯Fγ1,where
(70)WV¯=12V¯−V¯0W0x,⋯,WNxT.Fγ1 is the diagonal matrix whose diagonal elements are the value of fγ1 at interpolating points, i.e.,
(71)Fγ1=fγ1xV¯0⋯0⋮⋱⋮0⋯fγ1xV¯N.Similarly,Fγ2∼Fγ3, FV1∼FV4, and Fh1∼Fh2 mentioned below are all the diagonal matrices whose diagonal elements are the value of corresponding function at interpolating points.Substituting (63) and (64) into (60), the analytical solution to the zeroth-order velocity increment can be obtained as
(72)ΔV0V¯=12V¯−V¯0Ve∑i=0NfV1xV¯iα2xV¯i+fV3xV¯iΔγ0xV¯iWix=12V¯−V¯0Ve∑i=0NfV1xV¯iα2xV¯i+fV3xV¯iγα10xV¯iαWix=Vα10V¯α+Vα20V¯αsqu,where αsqu is the vector consisting of the square of AOA at interpolating points, i.e.,
(73)αsqu=α2xV¯0,α2xV¯1,⋯,α2xV¯NT,Vα10V¯ and Vα20V¯ are the coefficient matrices expressed as
(74)Vα10V¯=VeWTV¯FV3AcbFγ1,Vα20V¯=VeWTV¯FV1,where Acb is
(75)Acb=WxV¯0,⋯,WxV¯NT.Substituting (65) and (66) into (61), the analytical solution to the zeroth-order altitude increment can be obtained as
(76)Δh0V¯=12V¯−V¯0r0∑i=0Nfh1xV¯iΔV¯0xV¯i+fh2xV¯iΔγ0xV¯iWix=12V¯−V¯0r0∑i=0Nfh1xV¯iVα10V¯α+Vα20V¯αsquVeWix+12V¯−V¯0r0∑i=0Nfh2xV¯iWTV¯Fγ1αWix=hα10V¯α+hα20V¯αsqu,where
(77)hα10V¯=r0WTV¯Fh1AcbFV3+Fh2AcbFγ1,hα20V¯=r0WV¯Fh1AcbFV1..


### 3.3.3. First-Order Correction to Trajectory Increment

As for the first-order corrections, their analytical solutions can be obtained by integrating the three differential equations of (58) in sequence. Similar to the zeroth-order solution, Chebyshev polynomials are also used to approximate the integrands. Therefore, the analytical solution to Δγ1 is
(78)Δγ1V¯=12V¯−V¯0∑i=0NWixfγ2xV¯iΔV0xV¯iVe+fγ3xV¯iΔγ0xV¯i+12V¯−V¯0∑i=0NWixfγ4xV¯iΔh0xV¯iαxV¯ir0fγ5xV¯iΔγ0xV¯i2.Similar to (76), the first term of (78) can be rewritten as
(79)Δγ11V¯=γα11V¯α+γα2V¯αsqu,where
(80)γα11V¯=WTV¯Fγ2AcbFV3+Fγ3AcbFγ1,γα2V¯=WTV¯Fγ2AcbFV1..Substituting (69) and (76) into the second term of (78), it can be expressed as
(81)Δγ21V¯=12V¯−V¯01r0∑i=0NWixfγ4xV¯ihα10xV¯iα+hα20xV¯iαsquαxV¯i+12V¯−V¯0∑i=0NWixfγ5xV¯iγα10xV¯iα2=αTγα3V¯α+αsquTγα4V¯α,where γα3V¯ and γα4V¯ are the coefficient matric expressed as
(82)γα3V¯=12V¯−V¯01r0∑i=0NWiV¯fγ4xV¯ihα10xV¯iTIi+12V¯−V¯0∑i=0NWiV¯fγ5xV¯iγα10xV¯iTγα10xV¯i,γα4V¯=12V¯−V¯01r0∑i=0NWiV¯fγ4xV¯ihα20xV¯iTIi,where Ii is the i-th row of the identity matrix with the dimension N+1. Then, combining (79) with (81), we get
(83)Δγ1V¯=γα11V¯α+γα2V¯αsqu+αTγα3V¯α+αsquTγα4V¯α.By using Chebyshev polynomials to approximate the integrands, the analytical solution toΔV1 can be obtained as
(84)ΔV1V¯=12V¯−V¯0Ve∑i=0NWixfV2xV¯iΔV0xV¯iVe+fV4xV¯iΔh0xV¯ir0+12V¯−V¯0Ve∑i=0NWixfV5xV¯iΔh0xV¯iα2xV¯ir0+fV6xV¯iΔγ0xV¯i2.Similar to (76), the first term of (84) can be rewritten as
(85)ΔV11V¯=Vα10V¯α+Vα20V¯αsqu,where
(86)Vα11V¯=VeWTV¯FV2AcbFV3+FV4AcbFh1AcbFV3+Fh2AcbFγ1,Vα21V¯=VeWTV¯FV2+FV4AcbFh1AcbFV1.Similarly, substituting (69) and (76) into the second term of (84), it can be expressed as
(87)ΔV21V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯ihα10xV¯iα+hα20xV¯iαsquα2xV¯i+12V¯−V¯0Ve∑i=0NWixfV6xV¯iγα10xV¯iα2=αTVα3V¯α+αsquTVα4V¯α+αsquTVα5V¯αsqu,where
(88)Vα3V¯=12V¯−V¯0Ve∑i=0NWixfV6xV¯iγα10xV¯iTγα10xV¯i,Vα4V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯iITihα10xV¯i,Vα5V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯iITihα20xV¯i.Then, combining (85) with (87), the analytical solution to ΔV1 is
(89)ΔV1V¯=Vα10V¯α+Vα20V¯αsqu+αTVα3V¯α+αsquTVα4V¯α+αsquTVα5V¯αsqu.Similar to (78) and (84), the analytical solution to Δh1 can be approximated as
(90)Δh1V¯=12V¯−V¯0r0∑i=0NWixfh1xV¯iΔV1xV¯iVe+r0∑i=0NWiV¯fh2xV¯iΔγ1xV¯i+12V¯−V¯0r0∑i=0NWixfh3xV¯iΔV0xV¯iΔγ0xV¯iVe+fh4xV¯iΔγ0xV¯i2.Substituting (89) into (90), the first term of (90) can be rewritten as
(91)Δh11V¯=12V¯−V¯0r0Ve∑i=0NWixfh1xV¯iVα11xV¯iα+Vα21xV¯iαsqu+αTVα3xV¯iα+αsquTVα4xV¯iα+αsquTVα5xV¯iαsqu=hα1,11V¯α+hα2,11V¯αsqu+αThα3,1V¯α+αsquThα4,1V¯α+αsquThα5V¯α,where
(92)hα1,11V¯=r0WTV¯Fh1AcbFV2AcbFV3+FV4AcbFh1AcbFV3+Fh2AcbFγ1,hα2,11V¯=r0WTV¯Fh1AcbFV2+FV4AcbFh1AcbFV1,hα3,1V¯=12V¯−V¯0r0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV6xV¯jγα10xV¯jTγα10xV¯j,hα4,1V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV5xV¯jITjhα10xV¯j,hα5V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV5xV¯jITjhα20xV¯j.Similar to (91), the second term of (90) can be rewritten as
(93)Δh21V¯=hα1,21V¯α+hα2,21V¯αsqu+αThα3,2V¯α+αsquThα4,2V¯α,where
(94)hα1,21V¯=r0WTV¯Fh2AcbFγ2AcbFV3+Fγ3AcbFγ1,hα2,21V¯=r0WTV¯Fh2AcbFγ2AcbFV1,hα3,2V¯=12V¯−V¯0r0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ5xV¯jγα10xV¯jTγα10xV¯j+12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ4xV¯jhα10xV¯jTIj,hα4,2V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ4xV¯jhα20xV¯jTIj.Substituting (69) and (72) into the third term of (90), it can be expressed as
(95)Δh31V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα10xV¯iα+Vα20xV¯iαsquγα10xV¯iα+12V¯−V¯0r0∑i=0NWixfh4xV¯iγα10xV¯iα2=αThα3,3V¯α+αsquThα4,3V¯α,where
(96)hα3,3V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα10xV¯iTγα10xV¯i+12V¯−V¯0r0∑i=0NWixfh4xV¯iγα10xV¯iTγα10xV¯i,hα4,3V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα20xV¯iTγα10xV¯i.Therefore, combining (91), (93), and (95), the analytical solution to Δh1 is expressed as
(97)Δh1V¯=hα11V¯α+hα21V¯αsqu+αThα3V¯α+αsquThα4V¯α+αsquThα5V¯α,where
(98)hα11V¯=hα1,11V¯+hα1,21V¯,hα21V¯=hα2,11V¯+hα2,21V¯,hα3V¯=hα3,1V¯+hα3,2V¯+hα3,3V¯,hα4V¯=hα4,1V¯+hα4,2V¯+hα4,3V¯..


### 3.3.4. Full Analytical Solution to Trajectory with Arbitrary AOA

When the zero-order and first-order solutions are derived, the full analytical solutions to trajectory increment can be formulated as(99)ΔxV¯=Xα1V¯α+Xα2V¯αsqu+ATXα3V¯α+AsquTXα4V¯α+AsquTXα5V¯αsqu,where A and Asqu are the partitioned diagonal matrices expressed as (100), respectively.
(100)A=diagα,α,α,Asqu=diagαsqu,αsqu,αsqu.The coefficients are(101)XαiV¯=VαiTV¯,γαiTV¯,hαiTV¯T,i=1,2,3,4,5,and each coefficient is determined by the sum of the zeroth-order and first-order terms.Then, the analytical solution to ascent trajectory can be expressed as(102)x=xbV¯+Xα1V¯α+Xα2V¯αsqu+ATXα3V¯α+AsquTXα4V¯α+AsquTXα5V¯αsqu,where xbV¯ is the analytical solution to the trajectory with zero AOA, which is expressed as
(103)xbV¯=Vb0+Vb1,γb0+γb1,hb0+hb1T.It should be noted that the analytical solution includes the high order of AOA which leads to higher accuracy for states in all region than that only first order of AOA is considered.


## 3.3.1. Regular Perturbation Model with AOA Increment

Because the independent variable for the analytical solution with zero AOA is the dimensionless velocity,V¯, the differential equations of trajectory increments should be formulated as
(53)dΔVdV¯=VefV1α2+fV2ΔVVe+fV3Δγ+fV4Δhr0+fV5Δhα2r0+fV6Δγ2,dΔγdV¯=fγ1α+fγ2ΔVVe+fγ3Δγ+fγ4Δhαr0+fγ5Δγ2,dΔh¯dV¯=r0fh1ΔVVe+fh2Δγ+fh3ΔVΔγVe+fh4Δγ2,,where fV1∼fV6, fγ1∼fγ5, and fh1∼fh4 are all the functions of V¯, which are expressed as
(54)fV1V¯=−KP2+KaC¯dαV¯b2exp−βΔhb,fV2V¯=−2KaC¯d0V¯bexp−βΔhb,fV3V¯=−K2cosγbexp−V¯,fV4V¯=βr0KaCd0V¯b2exp−βΔhb,fV5V¯=βr0KaC¯dαV¯b2exp−βΔhb,fV6V¯=K2sinγbexp−V¯2,fγ1V¯=KPV¯b+KaC¯lαV¯bexp−βΔhb,fγ2V¯=K1+K2V¯b2cosγbexp−V¯,fγ3V¯=−K1V¯b−K2V¯bsinγbexp−V¯,fγ4V¯=−βr0KaC¯lαV¯bexp−βΔhb,fγ5V¯=−K1V¯b−K2/V¯bcosγbexp−V¯2,fh1V¯=K1sinγbexp−V¯,fh2V¯=K1V¯bcosγbexp−V¯,fh3V¯=K1cosγbexp−V¯,fh4V¯=−K1V¯bsinγbexp−V¯2..Generally, the AOA profile is expressed as the function of the flight time. The AOA profile with the dimensionless velocity can be formulated by using the relationship between the flight timet and V¯, which can be easily expressed as
(55)t=t0+T1−expV¯0−V¯.In the differential equation ofΔγ, it is easy to find that the term, fγ1α, is much larger than the sum of other terms, fγ1α≫fγ2ΔV/Ve+fγ3Δγ+fγ4Δhα/r0+fγ5Δγ2, which is the dominant term. Similarly, the term, fV1α2+fV3Δγ, is the dominant term of the differential equation of ΔV. As for the differential equation of Δh, due to Vb≫ΔV and Δγ2≤1, it is obvious that the term, fh1ΔV/Ve+fh2Δγ, is the dominant term. Therefore, the dynamics can be modified as
(56)dΔγdV¯dΔVdV¯dΔhdV¯=fγ1αVefV1α2+VefV3Δγr0Vefh1ΔV+r0fh2Δγ+εk1Vefγ2ΔV+fγ3Δγ+1r0fγ4Δhα+fγ5Δγ2fV2ΔV+Ver0fV4Δh+Ver0fV5Δhα2+VefV6Δγ2r0Vefh3ΔVΔγ+r0fh4Δγ2.According to the regular perturbation method presented in AppendixB, the zeroth-order and first-order dynamics are
(57)dΔγ0dV¯=fγ1V¯α,dΔV0dV¯=VefV1V¯α2+fV3V¯Δγ0,dΔh0dV¯=r0fh1V¯ΔV0Ve+fh2V¯Δγ0,,(58)dΔγ1dV¯=fγ2V¯ΔV0Ve+fγ3V¯Δγ0+fγ4V¯Δh0αr0+fγ5V¯Δγ02,dΔV1dV¯=VefV2V¯ΔV0Ve+fV4V¯Δh0r0+fV5V¯Δh0α2r0+fV6V¯Δγ02,dΔh1dV¯=r0fh1V¯ΔV1Ve+fh2V¯Δγ1+fh3V¯ΔV0Δγ0Ve+fh3V¯Δγ02..The initial values of zeroth-order and first-order states are all zeros.


## 3.3.2. Zeroth-Order Solution to Trajectory Increment

It is obvious that the three differential equations of zero-order states are not coupled. Therefore, by integrating these differential equations orderly, the analytical solutions toΔγ0, ΔV0, and Δh0 can be obtained as
(59)Δγ0V¯=∫V¯0V¯fγ1xV¯αxV¯dxV¯,(60)ΔV0V¯=Ve∫V¯0V¯fV1xV¯α2xV¯dxV¯+Ve∫V¯0V¯fV3xV¯Δγ0xV¯dxV¯,(61)Δh0V¯=r0Ve∫V¯0V¯fh1xV¯ΔV0xV¯dxV¯+r0∫V¯0V¯fh2xV¯Δγ0xV¯dxV¯.Although the integral terms in the above formulas are quite complex, the integrands are related only to the independent variable,V¯. Here, Chebyshev polynomials are also used to approximate the integrands. Let
(62)fγ1xV¯αxV¯≈∑i=0Nfγ1xV¯iαxV¯iLix,(63)fV1xV¯α2xV¯≈∑i=0NfV1xV¯iα2xV¯iLix,(64)fV3xV¯Δγ0xV¯≈∑i=0NfV3xV¯iΔγ0xV¯iLix,(65)fh1xV¯ΔV¯0xV¯≈∑i=0Nfh1xV¯iΔV¯0xV¯iLix,(66)fh2xV¯Δγ0xV¯≈∑i=0Nfh2xV¯iΔγ0xV¯iLix.Substituting (62) into (59), the analytical solution to the zeroth-order flight-path angle increment can be obtained as
(67)Δγ0V¯=12V¯−V¯0∑i=0Nfγ1xV¯iαxV¯iWix=γα10V¯α,where α is the vector consisting of the value of AOA at interpolating points, i.e.,
(68)α=αxV¯0,αxV¯1,⋯,αxV¯NT,γα10V¯ is the coefficient matrix expressed as
(69)γα10V¯=12V¯−V¯0W0V¯fγ1xV¯0,⋯,WNV¯fγ1xV¯N=WTV¯Fγ1,where
(70)WV¯=12V¯−V¯0W0x,⋯,WNxT.Fγ1 is the diagonal matrix whose diagonal elements are the value of fγ1 at interpolating points, i.e.,
(71)Fγ1=fγ1xV¯0⋯0⋮⋱⋮0⋯fγ1xV¯N.Similarly,Fγ2∼Fγ3, FV1∼FV4, and Fh1∼Fh2 mentioned below are all the diagonal matrices whose diagonal elements are the value of corresponding function at interpolating points.Substituting (63) and (64) into (60), the analytical solution to the zeroth-order velocity increment can be obtained as
(72)ΔV0V¯=12V¯−V¯0Ve∑i=0NfV1xV¯iα2xV¯i+fV3xV¯iΔγ0xV¯iWix=12V¯−V¯0Ve∑i=0NfV1xV¯iα2xV¯i+fV3xV¯iγα10xV¯iαWix=Vα10V¯α+Vα20V¯αsqu,where αsqu is the vector consisting of the square of AOA at interpolating points, i.e.,
(73)αsqu=α2xV¯0,α2xV¯1,⋯,α2xV¯NT,Vα10V¯ and Vα20V¯ are the coefficient matrices expressed as
(74)Vα10V¯=VeWTV¯FV3AcbFγ1,Vα20V¯=VeWTV¯FV1,where Acb is
(75)Acb=WxV¯0,⋯,WxV¯NT.Substituting (65) and (66) into (61), the analytical solution to the zeroth-order altitude increment can be obtained as
(76)Δh0V¯=12V¯−V¯0r0∑i=0Nfh1xV¯iΔV¯0xV¯i+fh2xV¯iΔγ0xV¯iWix=12V¯−V¯0r0∑i=0Nfh1xV¯iVα10V¯α+Vα20V¯αsquVeWix+12V¯−V¯0r0∑i=0Nfh2xV¯iWTV¯Fγ1αWix=hα10V¯α+hα20V¯αsqu,where
(77)hα10V¯=r0WTV¯Fh1AcbFV3+Fh2AcbFγ1,hα20V¯=r0WV¯Fh1AcbFV1..


## 3.3.3. First-Order Correction to Trajectory Increment

As for the first-order corrections, their analytical solutions can be obtained by integrating the three differential equations of (58) in sequence. Similar to the zeroth-order solution, Chebyshev polynomials are also used to approximate the integrands. Therefore, the analytical solution to Δγ1 is
(78)Δγ1V¯=12V¯−V¯0∑i=0NWixfγ2xV¯iΔV0xV¯iVe+fγ3xV¯iΔγ0xV¯i+12V¯−V¯0∑i=0NWixfγ4xV¯iΔh0xV¯iαxV¯ir0fγ5xV¯iΔγ0xV¯i2.Similar to (76), the first term of (78) can be rewritten as
(79)Δγ11V¯=γα11V¯α+γα2V¯αsqu,where
(80)γα11V¯=WTV¯Fγ2AcbFV3+Fγ3AcbFγ1,γα2V¯=WTV¯Fγ2AcbFV1..Substituting (69) and (76) into the second term of (78), it can be expressed as
(81)Δγ21V¯=12V¯−V¯01r0∑i=0NWixfγ4xV¯ihα10xV¯iα+hα20xV¯iαsquαxV¯i+12V¯−V¯0∑i=0NWixfγ5xV¯iγα10xV¯iα2=αTγα3V¯α+αsquTγα4V¯α,where γα3V¯ and γα4V¯ are the coefficient matric expressed as
(82)γα3V¯=12V¯−V¯01r0∑i=0NWiV¯fγ4xV¯ihα10xV¯iTIi+12V¯−V¯0∑i=0NWiV¯fγ5xV¯iγα10xV¯iTγα10xV¯i,γα4V¯=12V¯−V¯01r0∑i=0NWiV¯fγ4xV¯ihα20xV¯iTIi,where Ii is the i-th row of the identity matrix with the dimension N+1. Then, combining (79) with (81), we get
(83)Δγ1V¯=γα11V¯α+γα2V¯αsqu+αTγα3V¯α+αsquTγα4V¯α.By using Chebyshev polynomials to approximate the integrands, the analytical solution toΔV1 can be obtained as
(84)ΔV1V¯=12V¯−V¯0Ve∑i=0NWixfV2xV¯iΔV0xV¯iVe+fV4xV¯iΔh0xV¯ir0+12V¯−V¯0Ve∑i=0NWixfV5xV¯iΔh0xV¯iα2xV¯ir0+fV6xV¯iΔγ0xV¯i2.Similar to (76), the first term of (84) can be rewritten as
(85)ΔV11V¯=Vα10V¯α+Vα20V¯αsqu,where
(86)Vα11V¯=VeWTV¯FV2AcbFV3+FV4AcbFh1AcbFV3+Fh2AcbFγ1,Vα21V¯=VeWTV¯FV2+FV4AcbFh1AcbFV1.Similarly, substituting (69) and (76) into the second term of (84), it can be expressed as
(87)ΔV21V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯ihα10xV¯iα+hα20xV¯iαsquα2xV¯i+12V¯−V¯0Ve∑i=0NWixfV6xV¯iγα10xV¯iα2=αTVα3V¯α+αsquTVα4V¯α+αsquTVα5V¯αsqu,where
(88)Vα3V¯=12V¯−V¯0Ve∑i=0NWixfV6xV¯iγα10xV¯iTγα10xV¯i,Vα4V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯iITihα10xV¯i,Vα5V¯=12V¯−V¯0Ver0∑i=0NWixfV5xV¯iITihα20xV¯i.Then, combining (85) with (87), the analytical solution to ΔV1 is
(89)ΔV1V¯=Vα10V¯α+Vα20V¯αsqu+αTVα3V¯α+αsquTVα4V¯α+αsquTVα5V¯αsqu.Similar to (78) and (84), the analytical solution to Δh1 can be approximated as
(90)Δh1V¯=12V¯−V¯0r0∑i=0NWixfh1xV¯iΔV1xV¯iVe+r0∑i=0NWiV¯fh2xV¯iΔγ1xV¯i+12V¯−V¯0r0∑i=0NWixfh3xV¯iΔV0xV¯iΔγ0xV¯iVe+fh4xV¯iΔγ0xV¯i2.Substituting (89) into (90), the first term of (90) can be rewritten as
(91)Δh11V¯=12V¯−V¯0r0Ve∑i=0NWixfh1xV¯iVα11xV¯iα+Vα21xV¯iαsqu+αTVα3xV¯iα+αsquTVα4xV¯iα+αsquTVα5xV¯iαsqu=hα1,11V¯α+hα2,11V¯αsqu+αThα3,1V¯α+αsquThα4,1V¯α+αsquThα5V¯α,where
(92)hα1,11V¯=r0WTV¯Fh1AcbFV2AcbFV3+FV4AcbFh1AcbFV3+Fh2AcbFγ1,hα2,11V¯=r0WTV¯Fh1AcbFV2+FV4AcbFh1AcbFV1,hα3,1V¯=12V¯−V¯0r0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV6xV¯jγα10xV¯jTγα10xV¯j,hα4,1V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV5xV¯jITjhα10xV¯j,hα5V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh1xV¯ifV5xV¯jITjhα20xV¯j.Similar to (91), the second term of (90) can be rewritten as
(93)Δh21V¯=hα1,21V¯α+hα2,21V¯αsqu+αThα3,2V¯α+αsquThα4,2V¯α,where
(94)hα1,21V¯=r0WTV¯Fh2AcbFγ2AcbFV3+Fγ3AcbFγ1,hα2,21V¯=r0WTV¯Fh2AcbFγ2AcbFV1,hα3,2V¯=12V¯−V¯0r0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ5xV¯jγα10xV¯jTγα10xV¯j+12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ4xV¯jhα10xV¯jTIj,hα4,2V¯=12V¯−V¯0∑j=0N∑i=0NWixWjxV¯ifh2xV¯ifγ4xV¯jhα20xV¯jTIj.Substituting (69) and (72) into the third term of (90), it can be expressed as
(95)Δh31V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα10xV¯iα+Vα20xV¯iαsquγα10xV¯iα+12V¯−V¯0r0∑i=0NWixfh4xV¯iγα10xV¯iα2=αThα3,3V¯α+αsquThα4,3V¯α,where
(96)hα3,3V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα10xV¯iTγα10xV¯i+12V¯−V¯0r0∑i=0NWixfh4xV¯iγα10xV¯iTγα10xV¯i,hα4,3V¯=12V¯−V¯0r0Ve∑i=0NWixfh3xV¯iVα20xV¯iTγα10xV¯i.Therefore, combining (91), (93), and (95), the analytical solution to Δh1 is expressed as
(97)Δh1V¯=hα11V¯α+hα21V¯αsqu+αThα3V¯α+αsquThα4V¯α+αsquThα5V¯α,where
(98)hα11V¯=hα1,11V¯+hα1,21V¯,hα21V¯=hα2,11V¯+hα2,21V¯,hα3V¯=hα3,1V¯+hα3,2V¯+hα3,3V¯,hα4V¯=hα4,1V¯+hα4,2V¯+hα4,3V¯..


## 3.3.4. Full Analytical Solution to Trajectory with Arbitrary AOA

When the zero-order and first-order solutions are derived, the full analytical solutions to trajectory increment can be formulated as(99)ΔxV¯=Xα1V¯α+Xα2V¯αsqu+ATXα3V¯α+AsquTXα4V¯α+AsquTXα5V¯αsqu,where A and Asqu are the partitioned diagonal matrices expressed as (100), respectively.
(100)A=diagα,α,α,Asqu=diagαsqu,αsqu,αsqu.The coefficients are(101)XαiV¯=VαiTV¯,γαiTV¯,hαiTV¯T,i=1,2,3,4,5,and each coefficient is determined by the sum of the zeroth-order and first-order terms.Then, the analytical solution to ascent trajectory can be expressed as(102)x=xbV¯+Xα1V¯α+Xα2V¯αsqu+ATXα3V¯α+AsquTXα4V¯α+AsquTXα5V¯αsqu,where xbV¯ is the analytical solution to the trajectory with zero AOA, which is expressed as
(103)xbV¯=Vb0+Vb1,γb0+γb1,hb0+hb1T.It should be noted that the analytical solution includes the high order of AOA which leads to higher accuracy for states in all region than that only first order of AOA is considered.


## 4. Optimal AOA Correction with Terminal Constraint

In the last section, the analytical solutions to ascent trajectory with arbitrary AOA profile are derived. The optimal AOA correction can be developed to remove terminal errors. Now, let us consider the nonlinear optimal control problem with the performance index expressed as (4) which subjects to dynamical constraints presented in (1) and the hard terminal constraints (3).The dynamic constraints presented in (1) can be transformed into the following algebraic constraints using the analytic solutions previously derived:
(104)xtf=xbV¯f+Xα1V¯fα+Xα2V¯fαsqu+ATXα3V¯fα+AsquTXα4V¯fα+AsquTXα5V¯fαsqu,where V¯f is the value of V¯ at the terminal time. Then, substituting (104) into (3), the terminal constraint functions can be converted to the functions about the AOA vector, α:
(105)Φα=ψxtf=0.Furthermore, using Chebyshev polynomial approximation, the performance index expressed in (4) can also be converted into
(106)J=12TexpV¯0−n+1∫V¯0V¯fα2xV¯exp−xV¯exp−xV¯−exp−V¯fndxV¯=14V¯f−V¯0TexpV¯0−n+1∑i=0NWiexp−xV¯iα2xV¯iexp−xV¯i−exp−V¯fn=12αTQα,where Q∈RN×N is a positive definite diagonal matrix whose elements along the diagonal are
(107)Qi=12V¯f−V¯0TexpV¯0−n+1Wiexp−xV¯iexp−xV¯i−exp−V¯fn,i=0,1,⋯,N.Finally, the nonlinear optimal control problem is transformed into a nonlinear constrained optimization problem, of which the performance index and constraints are (106) and (104), respectively. And the optimal variables are the AOA at the interpolation points, αxV¯ii=0,1,⋯,N. The nonlinear constrained optimization problem can further be transformed into a series of quadratic programming problems if linearization is applied for the terminal constraints. The optimal variables are obtained by iteratively solving this series of quadratic programming problems. The derivation process is presented as follows.Expanding (105) in Taylor series around the trajectory with the nominal AOA curve, αp. Then, a set of linear terminal constraints (108) using the deviations from the nominal AOA as the independent variables are formulated by neglecting the higher order terms:
(108)ψxpf+∂ψ∂xfxpf∂δxf∂ααpδα=0,where αp is normal AOA vector whose elements are the value of nominal AOA profile at interpolating points, i.e.,
(109)αp=αpxV¯0,αpxV¯1,⋯,αpxV¯NT.δα=α−αp is the deviation of the real AOA from the nominal AOA; ∂ψ/∂xf⋅ is the partial derivative of the terminal constraint function ψ⋅ with respect to the terminal state xf; and ∂δxf/∂α⋅ is the partial derivative of the terminal state with respect to the vector of AOA at the interpolation points, α, which can be expressed as
(110)∂δxf∂αα=Xα1+2Xα2α~+ATXα3+X~α3+AsquTXα4+2ATX~α4α~+2AsquTXα5+X~α5α~,where α~=diagα and X~αi=Vαi,γαi,hαiT,i=3,4,5. xpf is the terminal value of trajectory with the nominal AOA profile:
(111)xpf=xbV¯f+Xα1V¯fαp+Xα2V¯fαp,squ+ApTXα3V¯fαp+Ap,squTXα4V¯fαp+Ap,squTXα5V¯fαp,squ.As for the performance index (106), using the δα as the optimal variables, it can also be rewritten as
(112)J=12αpTQαp+αpTQδα+12δαTQδα.Considering the terminal constraints, the augmented performance index is(113)J¯=J+υTψxpf+∂ψ∂xfxpf∂δxf∂ααpδα,where υ∈RS is the Lagrange multiplier vector. According to numerical optimization theory, the corresponding KKT conditions are
(114)∂J¯∂δα=Qαp+δα+∂ψ∂xfxpf∂δxf∂ααpTυ=0,∂J¯∂υ=ψxpf+ψ̇xpfδẋfαpδα=0.It is easy to find that the KKT conditions are a set of linear algebraic equations which can be expressed as(115)Sz=K,where z=δαT,υTT is the requested variables; The elements of S and K are defined as
(116)S=Q∂ψ∂xfxpf∂δxf∂ααpT∂ψ∂xfxpf∂δxf∂ααp0S×S,K=−Qαp;ψxpf.Note thatK includes the deviation of terminal constraints. It is easy to solve the above linear algebraic equations to get the improvement of AOA, δα. The updating AOA at the interpolation points is
(117)α=αp+δα.Then, using the previous calculation process, we iteratively calculateδα until ψxtf<ε, where ε is the prescribed tolerance. It is noted that the coefficients in (111) are dependent with the AOA profile; therefore, only one time is needed to calculate those coefficients at each guidance period. Additionally, the value of AOA at the arbitrary time t can be approximated by Chebyshev interpolating polynomial:
(118)αt=∑i=0NαxV¯iLixt,where xt is the value of x at the time t, which can be obtained by
(119)xt=V¯−V¯2V¯1=−2V¯f−V¯0ln1−t−t0T−1.


## 5. Implementation of Ascent Guidance Law

In this section, an ascent guidance law with terminal constraints is proposed, in which the current AOA command is generated by iteratively solving the nonlinear optimal control problem using the analytical trajectory prediction. The procedure for implementing the proposed guidance law is included in the flow chart in Figure3 and summarized in detail as follows.
(1)
In the initialization, set initial parameters and select zero AOA as the initial nominal AOAαp(2)
Calculate the analytical solutions toxbf and the coefficients Xα1~Xα5 using the current states(3)
Apply the method of linearization around the nominal AOAαp to formulate the matrices S and K, and calculate control deviation δα by solving linear algebraic equations(4)
Generate the new nominal AOA, and repeat step 3 untilψxtf is within the prescribed tolerance, then go to Step 5(5)
Update the current guidance commandαt0 by (118). Then, go to next step(6)
Applyαt0 to the real vehicle and then return Step 2. Also, the AOA obtained at Step 4 is used as the initial nominal AOA for the next guidance stepFigure 3
Flow chart for implementing the proposed guidance law.Different from the traditional guidance laws, the proposed guidance does not need to conduct any offline planning for specific mission. Therefore, this guidance can handle various urgent tasks. Besides, it can significantly reduce the computational load compared to the method based on the numerical solutions because the latter needs to conduct thousands of numerical trajectory integrations for iterations in practice.


## 6. Simulation and Analysis

In this section, a comparison with the numerical integration is provided to verified the accuracy and efficiency of the analytical trajectory prediction. Then, several nominal cases and Monte Carlo simulations with dispersions and uncertainties are carried out to evaluate the performance and robustness of the proposed method. In which, a typical launch vehicle described in previous section is used. All programs are implemented on a personal computer with 2.6GHz processor and 2020a MATLAB environment.


### 6.1. Accuracy Verification for Analytical Solution with Arbitrary AOA

The accuracy of proposed guidance law is highly dependent with the analytical trajectory prediction, which provides the nominal trajectory for Taylor expansion. Therefore, an accuracy verification for analytical solution is provided in this subsection. For all simulations, the initial conditions are set toV0=875m/s, γ0=20deg, and h0=32km. And three types of AOA curve listed in Table 2 are considered.Table 2
Angle of Attack curves.
CaseAOA curve1αt=α06.25t/tf−0.62−12αt=α0t/tf−123αt=α01−2t/tfFirstly, the cases of three types of AOA curve withα0=10°,5°,0°,−5°,−10° are simulated. Figure 4 shows the time histories of velocity, flight-path angle, and altitude for the different AOA curves. Obviously, the errors between the analytical solutions and numerical solutions are very small. In order to provide the quantitative description on the errors, the terminal states for analytical solutions and numerical solutions of α0=−10° are presented in Table 3. Table 3 shows that the relative errors are all less than 1%, which verifies the accuracy of analytical solutions. It is noted that the accuracy will be further improved as the total flight time decreases.Figure 4
Ascent state histories with various AOA curves.
(a)
AOA curve 1(b)
AOA curve 2(c)
AOA curve 3Table 3
Comparison of the analytical solutions and numerical solutions.
CaseMethodsVfm/sγfdeghfmComputing time (s)1Numerical solutions3225.576.438348217.100.0230Analytical solutions3227.146.449248279.288.03E−04Relative errors0.067%0.079%0.38%—2Numerical solutions3272.49-5.008934662.200.0267Analytical solutions3288.19-4.994834929.097.16E−04Relative errors0.48%-0.28%0.77%—3Numerical solutions3250.870.582238319.760.0234Analytical solutions3259.150.452838297.287.86E−04Relative errors0.35%0.67%-0.36%—Besides, the computing time of the analytical and numerical solutions obtained by fourth-order Runge-Kutta integral with the step length of 0.1 s is presented in Table3. The results indicate that the computing time of the analytical solutions is less than 1 ms, which is about 3% of that of the numerical solutions. Obviously, the computational load of the trajectory planning method based on the numerical solutions is heavy because thousands of numerical trajectory integrations in practice are needed. However, if the analytical solutions are applied to trajectory planning, only one calculation is needed so that the consumed time will be greatly reduced. Conclusively, the analytical trajectory prediction has high accuracy and efficiency.


### 6.2. Nominal Case with Different Terminal Constraints

In this subsection, optimality and applicability of the proposed guidance law are evaluated by several different terminal constraints. It is noted that the initial conditions are the same as the last subsection and the coefficientn is selected to be 0.5.


#### 6.2.1. Optimality Verification

In this subsection, the optimality of the proposed guidance law is verified by a comparison with the results generated by the world famous nonlinear optimal control solver software, GPOPS-II [41, 42]. There are four cases of terminal constraint which are listed in Table 4 are considered. It should be noted that all the results in this subsection are generated by open-loop optimal guidance.Table 4
Terminal constraints for different cases.
ConstraintsCase 1Case 2Case 3Case 4γfcdeg0055hfckm40454540Figures5 and 6 depict the altitude vs. flight-path angle histories and AOA profiles obtained by the two methods, respectively. It is obvious that the all the curves are the same, which implies that the proposed method is capable of providing the optimal trajectory. The terminal errors of both methods for different cases are presented in Table 5. It can be seen from the table that although the terminal errors of the proposed method are slightly larger than that of GPOPS-II, it is still quite small, just less than 0.2% except case 4. In case 4, the terminal relative errors are a little large that about 3%. This is because the AOA is larger than the other three cases, which results in the relatively larger errors for the analytical integration. The errors absolutely met the requirement of the online guidance. Moreover, the accuracy will be improved as the time-to-go decreases.Figure 5
Altitude vs. flight-path angle histories.Figure 6
AOA profiles.Table 5
Terminal relative errors and computation times for different cases.
CaseMethodFlight-path angle errorAltitude errorComputing time (s)1GPOPS-II0.0025%-0.0053%0.171Proposed method0.063%-0.042%2.2E-032GPOPS-II-0.0036%0.0042%0.119Proposed method-0.069%0.11%2.5E-033GPOPS-II0.0044%-0.0032%0.280Proposed method-0.078%0.15%2.8E-034GPOPS-II0.0034%-0.0038%0.6772Proposed method2.42%-2.97%3.0E-03Table5 also lists the computing time of both methods. Though the accuracy of the proposed method is slightly lower than that of GPOPS-II, its computing time is less than 3 ms, which is much less than that of GPOPS-II.


#### 6.2.2. Applicability of Different Cases

In this subsection, the applicability of the proposed guidance law is tested by different initial conditions and terminal constraints, which are listed in Table6. All the cases can be divided into two groups. The first one is cases 1-4, in which the initial flight-path angle and altitude are same and the terminal constraints are different. And the second one is cases 5-8, in which the initial flight-path angle and altitude are different but the terminal constraints are same. Noted that the initial velocity is the same and set to 875 m/s2 in all these cases. In these simulations, the real thrust of vehicle is regarded as the expression of (5), and the aerodynamic coefficients, Clα, Cd0, and Cdα, are approximated by the polynomials of Mach number as shown in (7). Additionally, International Standard Atmosphere (ISA) mode, in which the density and pressure are the piecewise functions of altitude, is adopted.Table 6
Initial conditions and terminal constraints for different cases.
Caseγ0degh0kmγfcdeghfckm1203204022032045320325404203254551537040625370407152704082527040The flight-path angle and altitude histories are depicted in Figure7. It is obvious that the proposed guidance is able to guide the vehicle to the desired terminal states with high accuracy. The terminal errors are presented in Table 7. The maximum error for all cases is 2.04E−04 deg for fight-path angle and 8.8906e−04 m for altitude. It is very small and meets the requirements on the subsequent successful handover. Note that the flight-path angle curves for cases 1 and 2 are monotonous, but the ones for cases 3 and 4 decrease at the beginning and then increase at the end. The reason is that the terminal altitude constrains for cases 3 and 4 are relatively low so that the flight-path angles need to decrease at the beginning so as to meet the altitude constraints. That is consistent with the requirement of the energy management. As for cases 5 and 6, the altitude increases at the beginning and then decreases at the end. This is because the initial altitudes for these two case are large; the vehicle needs to climb and then descend to meet the terminal altitude constraints.Figure 7
Flight-path angle and altitude histories for different cases.
(a)
Flight-path angle histories(b)
Altitude historiesTable 7
Terminal errors with close-loop guidance.
CaseFlight-path angle error (deg)Altitude error (m)CaseFlight-path angle error (deg)Altitude error (m)1−1.15E−041.71E−035−1.7930e−051.9401e−042−2.04E−04−4.94E−0464.2814e−058.8906e−043−8.29E−05−3.28E−0471.0245e−041.5770e−0441.04E−046.61E−0582.7622e−052.6255e−05The AOA profiles obtained by proposed guidance are presented in Figure8. Obviously, the terminal AOA for all cases converge to zero, which result from the weighted squared sum performance index whose denominator is flight time to go. Additionally, all AOA profiles vary smoothly and its maximal value happens at the beginning. Moreover, as seen from Figure 9, let us focus on the consuming time for the proposed method at each guidance period. It is apparent that most of them distribute from 1 ms to 3 ms, and its maximum is only 5 ms. Therefore, the proposed guidance law has high computational efficiency and guidance accuracy.Figure 8
AOA profiles for different cases.Figure 9
Consuming time at each guidance period.


### 6.3. Monte Carlo Simulations

In this subsection, 1000 times Monte Carlo simulations are performed to further verify the robustness of the guidance law. Dispersions and uncertainties including aerodynamics error, atmospheric uncertainties, thrust uncertainties, and wind interference are fully considered. In these cases, the percentage deviations of the lift and drag coefficients typically vary with Mach number and AOA. And the percentage deviations of the thrust typically vary with altitude. For simplicity, the linear dispersion model is used in the simulations.(120)δP=δP0+kδPh91000−h91000,δCl=δCl0+kδClMaMa−78+kδClα18πα−π18,δCd=δCd0+kδCdMaMa−78+kδCdα18πα−π18.Dispersion parameters used in simulations are listed in Table8. VWind represents the wind speed, and δρ represents the percentage deviation of the atmospheric density. The first four cases listed in Table 6 are simulated in this subsection.Table 8
The dispersion parameters of the Monte Carlo simulations.
Variables3σDistributionVariables3σDistributionδP015%GaussiankδPh15%GaussianδCl015%GaussianδCd015%GaussiankδClMa30%GaussiankδCdMa30%GaussiankδClα45%GaussiankδCdα45%GaussianVwind300 m/sGaussianδρ30%GaussianThe dispersion of terminal states are presented in Figure10. The terminal altitude errors are mostly near 0.01 m, and the maximum is no more than 0.02 m. The terminal flight-path angle errors are all within 6×10-3 deg. The statistical distribution is further demonstrated in Table 9. The means and standard deviations of the terminal constraints are all within a small range, which reveals that the proposed phase guidance has great robustness and performance in guidance accuracy.Figure 10
Terminal deviations of the Monte Carlo simulations.
(a)
Case 1(b)
Case 2(c)
Case 3(d)
Case 4Table 9
The statistics on the terminal deviations of the Monte Carlo simulations.
CaseMeanδγfdegMeanδhfmStdδγfdegStdδhfm14.90E−041.28E−039.36E−058.86E−042−7.68E−04−2.47E−049.58E−053.86E−033−2.91E−051.02E−031.73E−044.62E−0342.61E−057.10E−041.62E−043.4E−03


## 6.1. Accuracy Verification for Analytical Solution with Arbitrary AOA

The accuracy of proposed guidance law is highly dependent with the analytical trajectory prediction, which provides the nominal trajectory for Taylor expansion. Therefore, an accuracy verification for analytical solution is provided in this subsection. For all simulations, the initial conditions are set toV0=875m/s, γ0=20deg, and h0=32km. And three types of AOA curve listed in Table 2 are considered.Table 2
Angle of Attack curves.
CaseAOA curve1αt=α06.25t/tf−0.62−12αt=α0t/tf−123αt=α01−2t/tfFirstly, the cases of three types of AOA curve withα0=10°,5°,0°,−5°,−10° are simulated. Figure 4 shows the time histories of velocity, flight-path angle, and altitude for the different AOA curves. Obviously, the errors between the analytical solutions and numerical solutions are very small. In order to provide the quantitative description on the errors, the terminal states for analytical solutions and numerical solutions of α0=−10° are presented in Table 3. Table 3 shows that the relative errors are all less than 1%, which verifies the accuracy of analytical solutions. It is noted that the accuracy will be further improved as the total flight time decreases.Figure 4
Ascent state histories with various AOA curves.
(a)
AOA curve 1(b)
AOA curve 2(c)
AOA curve 3Table 3
Comparison of the analytical solutions and numerical solutions.
CaseMethodsVfm/sγfdeghfmComputing time (s)1Numerical solutions3225.576.438348217.100.0230Analytical solutions3227.146.449248279.288.03E−04Relative errors0.067%0.079%0.38%—2Numerical solutions3272.49-5.008934662.200.0267Analytical solutions3288.19-4.994834929.097.16E−04Relative errors0.48%-0.28%0.77%—3Numerical solutions3250.870.582238319.760.0234Analytical solutions3259.150.452838297.287.86E−04Relative errors0.35%0.67%-0.36%—Besides, the computing time of the analytical and numerical solutions obtained by fourth-order Runge-Kutta integral with the step length of 0.1 s is presented in Table3. The results indicate that the computing time of the analytical solutions is less than 1 ms, which is about 3% of that of the numerical solutions. Obviously, the computational load of the trajectory planning method based on the numerical solutions is heavy because thousands of numerical trajectory integrations in practice are needed. However, if the analytical solutions are applied to trajectory planning, only one calculation is needed so that the consumed time will be greatly reduced. Conclusively, the analytical trajectory prediction has high accuracy and efficiency.


## 6.2. Nominal Case with Different Terminal Constraints

In this subsection, optimality and applicability of the proposed guidance law are evaluated by several different terminal constraints. It is noted that the initial conditions are the same as the last subsection and the coefficientn is selected to be 0.5.


### 6.2.1. Optimality Verification

In this subsection, the optimality of the proposed guidance law is verified by a comparison with the results generated by the world famous nonlinear optimal control solver software, GPOPS-II [41, 42]. There are four cases of terminal constraint which are listed in Table 4 are considered. It should be noted that all the results in this subsection are generated by open-loop optimal guidance.Table 4
Terminal constraints for different cases.
ConstraintsCase 1Case 2Case 3Case 4γfcdeg0055hfckm40454540Figures5 and 6 depict the altitude vs. flight-path angle histories and AOA profiles obtained by the two methods, respectively. It is obvious that the all the curves are the same, which implies that the proposed method is capable of providing the optimal trajectory. The terminal errors of both methods for different cases are presented in Table 5. It can be seen from the table that although the terminal errors of the proposed method are slightly larger than that of GPOPS-II, it is still quite small, just less than 0.2% except case 4. In case 4, the terminal relative errors are a little large that about 3%. This is because the AOA is larger than the other three cases, which results in the relatively larger errors for the analytical integration. The errors absolutely met the requirement of the online guidance. Moreover, the accuracy will be improved as the time-to-go decreases.Figure 5
Altitude vs. flight-path angle histories.Figure 6
AOA profiles.Table 5
Terminal relative errors and computation times for different cases.
CaseMethodFlight-path angle errorAltitude errorComputing time (s)1GPOPS-II0.0025%-0.0053%0.171Proposed method0.063%-0.042%2.2E-032GPOPS-II-0.0036%0.0042%0.119Proposed method-0.069%0.11%2.5E-033GPOPS-II0.0044%-0.0032%0.280Proposed method-0.078%0.15%2.8E-034GPOPS-II0.0034%-0.0038%0.6772Proposed method2.42%-2.97%3.0E-03Table5 also lists the computing time of both methods. Though the accuracy of the proposed method is slightly lower than that of GPOPS-II, its computing time is less than 3 ms, which is much less than that of GPOPS-II.


### 6.2.2. Applicability of Different Cases

In this subsection, the applicability of the proposed guidance law is tested by different initial conditions and terminal constraints, which are listed in Table6. All the cases can be divided into two groups. The first one is cases 1-4, in which the initial flight-path angle and altitude are same and the terminal constraints are different. And the second one is cases 5-8, in which the initial flight-path angle and altitude are different but the terminal constraints are same. Noted that the initial velocity is the same and set to 875 m/s2 in all these cases. In these simulations, the real thrust of vehicle is regarded as the expression of (5), and the aerodynamic coefficients, Clα, Cd0, and Cdα, are approximated by the polynomials of Mach number as shown in (7). Additionally, International Standard Atmosphere (ISA) mode, in which the density and pressure are the piecewise functions of altitude, is adopted.Table 6
Initial conditions and terminal constraints for different cases.
Caseγ0degh0kmγfcdeghfckm1203204022032045320325404203254551537040625370407152704082527040The flight-path angle and altitude histories are depicted in Figure7. It is obvious that the proposed guidance is able to guide the vehicle to the desired terminal states with high accuracy. The terminal errors are presented in Table 7. The maximum error for all cases is 2.04E−04 deg for fight-path angle and 8.8906e−04 m for altitude. It is very small and meets the requirements on the subsequent successful handover. Note that the flight-path angle curves for cases 1 and 2 are monotonous, but the ones for cases 3 and 4 decrease at the beginning and then increase at the end. The reason is that the terminal altitude constrains for cases 3 and 4 are relatively low so that the flight-path angles need to decrease at the beginning so as to meet the altitude constraints. That is consistent with the requirement of the energy management. As for cases 5 and 6, the altitude increases at the beginning and then decreases at the end. This is because the initial altitudes for these two case are large; the vehicle needs to climb and then descend to meet the terminal altitude constraints.Figure 7
Flight-path angle and altitude histories for different cases.
(a)
Flight-path angle histories(b)
Altitude historiesTable 7
Terminal errors with close-loop guidance.
CaseFlight-path angle error (deg)Altitude error (m)CaseFlight-path angle error (deg)Altitude error (m)1−1.15E−041.71E−035−1.7930e−051.9401e−042−2.04E−04−4.94E−0464.2814e−058.8906e−043−8.29E−05−3.28E−0471.0245e−041.5770e−0441.04E−046.61E−0582.7622e−052.6255e−05The AOA profiles obtained by proposed guidance are presented in Figure8. Obviously, the terminal AOA for all cases converge to zero, which result from the weighted squared sum performance index whose denominator is flight time to go. Additionally, all AOA profiles vary smoothly and its maximal value happens at the beginning. Moreover, as seen from Figure 9, let us focus on the consuming time for the proposed method at each guidance period. It is apparent that most of them distribute from 1 ms to 3 ms, and its maximum is only 5 ms. Therefore, the proposed guidance law has high computational efficiency and guidance accuracy.Figure 8
AOA profiles for different cases.Figure 9
Consuming time at each guidance period.


## 6.2.1. Optimality Verification

In this subsection, the optimality of the proposed guidance law is verified by a comparison with the results generated by the world famous nonlinear optimal control solver software, GPOPS-II [41, 42]. There are four cases of terminal constraint which are listed in Table 4 are considered. It should be noted that all the results in this subsection are generated by open-loop optimal guidance.Table 4
Terminal constraints for different cases.
ConstraintsCase 1Case 2Case 3Case 4γfcdeg0055hfckm40454540Figures5 and 6 depict the altitude vs. flight-path angle histories and AOA profiles obtained by the two methods, respectively. It is obvious that the all the curves are the same, which implies that the proposed method is capable of providing the optimal trajectory. The terminal errors of both methods for different cases are presented in Table 5. It can be seen from the table that although the terminal errors of the proposed method are slightly larger than that of GPOPS-II, it is still quite small, just less than 0.2% except case 4. In case 4, the terminal relative errors are a little large that about 3%. This is because the AOA is larger than the other three cases, which results in the relatively larger errors for the analytical integration. The errors absolutely met the requirement of the online guidance. Moreover, the accuracy will be improved as the time-to-go decreases.Figure 5
Altitude vs. flight-path angle histories.Figure 6
AOA profiles.Table 5
Terminal relative errors and computation times for different cases.
CaseMethodFlight-path angle errorAltitude errorComputing time (s)1GPOPS-II0.0025%-0.0053%0.171Proposed method0.063%-0.042%2.2E-032GPOPS-II-0.0036%0.0042%0.119Proposed method-0.069%0.11%2.5E-033GPOPS-II0.0044%-0.0032%0.280Proposed method-0.078%0.15%2.8E-034GPOPS-II0.0034%-0.0038%0.6772Proposed method2.42%-2.97%3.0E-03Table5 also lists the computing time of both methods. Though the accuracy of the proposed method is slightly lower than that of GPOPS-II, its computing time is less than 3 ms, which is much less than that of GPOPS-II.


## 6.2.2. Applicability of Different Cases

In this subsection, the applicability of the proposed guidance law is tested by different initial conditions and terminal constraints, which are listed in Table6. All the cases can be divided into two groups. The first one is cases 1-4, in which the initial flight-path angle and altitude are same and the terminal constraints are different. And the second one is cases 5-8, in which the initial flight-path angle and altitude are different but the terminal constraints are same. Noted that the initial velocity is the same and set to 875 m/s2 in all these cases. In these simulations, the real thrust of vehicle is regarded as the expression of (5), and the aerodynamic coefficients, Clα, Cd0, and Cdα, are approximated by the polynomials of Mach number as shown in (7). Additionally, International Standard Atmosphere (ISA) mode, in which the density and pressure are the piecewise functions of altitude, is adopted.Table 6
Initial conditions and terminal constraints for different cases.
Caseγ0degh0kmγfcdeghfckm1203204022032045320325404203254551537040625370407152704082527040The flight-path angle and altitude histories are depicted in Figure7. It is obvious that the proposed guidance is able to guide the vehicle to the desired terminal states with high accuracy. The terminal errors are presented in Table 7. The maximum error for all cases is 2.04E−04 deg for fight-path angle and 8.8906e−04 m for altitude. It is very small and meets the requirements on the subsequent successful handover. Note that the flight-path angle curves for cases 1 and 2 are monotonous, but the ones for cases 3 and 4 decrease at the beginning and then increase at the end. The reason is that the terminal altitude constrains for cases 3 and 4 are relatively low so that the flight-path angles need to decrease at the beginning so as to meet the altitude constraints. That is consistent with the requirement of the energy management. As for cases 5 and 6, the altitude increases at the beginning and then decreases at the end. This is because the initial altitudes for these two case are large; the vehicle needs to climb and then descend to meet the terminal altitude constraints.Figure 7
Flight-path angle and altitude histories for different cases.
(a)
Flight-path angle histories(b)
Altitude historiesTable 7
Terminal errors with close-loop guidance.
CaseFlight-path angle error (deg)Altitude error (m)CaseFlight-path angle error (deg)Altitude error (m)1−1.15E−041.71E−035−1.7930e−051.9401e−042−2.04E−04−4.94E−0464.2814e−058.8906e−043−8.29E−05−3.28E−0471.0245e−041.5770e−0441.04E−046.61E−0582.7622e−052.6255e−05The AOA profiles obtained by proposed guidance are presented in Figure8. Obviously, the terminal AOA for all cases converge to zero, which result from the weighted squared sum performance index whose denominator is flight time to go. Additionally, all AOA profiles vary smoothly and its maximal value happens at the beginning. Moreover, as seen from Figure 9, let us focus on the consuming time for the proposed method at each guidance period. It is apparent that most of them distribute from 1 ms to 3 ms, and its maximum is only 5 ms. Therefore, the proposed guidance law has high computational efficiency and guidance accuracy.Figure 8
AOA profiles for different cases.Figure 9
Consuming time at each guidance period.


## 6.3. Monte Carlo Simulations

In this subsection, 1000 times Monte Carlo simulations are performed to further verify the robustness of the guidance law. Dispersions and uncertainties including aerodynamics error, atmospheric uncertainties, thrust uncertainties, and wind interference are fully considered. In these cases, the percentage deviations of the lift and drag coefficients typically vary with Mach number and AOA. And the percentage deviations of the thrust typically vary with altitude. For simplicity, the linear dispersion model is used in the simulations.(120)δP=δP0+kδPh91000−h91000,δCl=δCl0+kδClMaMa−78+kδClα18πα−π18,δCd=δCd0+kδCdMaMa−78+kδCdα18πα−π18.Dispersion parameters used in simulations are listed in Table8. VWind represents the wind speed, and δρ represents the percentage deviation of the atmospheric density. The first four cases listed in Table 6 are simulated in this subsection.Table 8
The dispersion parameters of the Monte Carlo simulations.
Variables3σDistributionVariables3σDistributionδP015%GaussiankδPh15%GaussianδCl015%GaussianδCd015%GaussiankδClMa30%GaussiankδCdMa30%GaussiankδClα45%GaussiankδCdα45%GaussianVwind300 m/sGaussianδρ30%GaussianThe dispersion of terminal states are presented in Figure10. The terminal altitude errors are mostly near 0.01 m, and the maximum is no more than 0.02 m. The terminal flight-path angle errors are all within 6×10-3 deg. The statistical distribution is further demonstrated in Table 9. The means and standard deviations of the terminal constraints are all within a small range, which reveals that the proposed phase guidance has great robustness and performance in guidance accuracy.Figure 10
Terminal deviations of the Monte Carlo simulations.
(a)
Case 1(b)
Case 2(c)
Case 3(d)
Case 4Table 9
The statistics on the terminal deviations of the Monte Carlo simulations.
CaseMeanδγfdegMeanδhfmStdδγfdegStdδhfm14.90E−041.28E−039.36E−058.86E−042−7.68E−04−2.47E−049.58E−053.86E−033−2.91E−051.02E−031.73E−044.62E−0342.61E−057.10E−041.62E−043.4E−03


## 7. Conclusions

In this paper, an endoatmospheric ascent optimal guidance law within the framework of predictor-corrector algorithm is proposed. This method uses a regular perturbation method and pseudospectral collocation scheme to successfully derive the precise analytical nonlinear prediction of the ascent trajectory with arbitrary AOA, which significantly improves the computational efficiency in comparison with the traditional method. Then, an iterative scheme is used to solve the nonlinear optimal control problem around the analytical trajectory so as to provide the optimal AOA profile with high accuracy. Several nominal and Monto Carlo simulations are carried out to evaluate the performance and robustness of the proposed guidance law. The results show that this method is applicable for various cases with large difference, the consuming time at each guidance period is less than 5 ms, and the maximum terminal error is less than1.0E−3 deg for flight-path angle and 1.0E−2 m for altitude even in large dispersions and uncertainties. The optimality of the proposed method is verified by comparing the AOA profiles with the optimal ones given by famous nonlinear optimal control software GPOPS-II. Conclusively, it has high guidance accuracy, strong robustness, and low computational burden, which has great potential to be applied as baseline algorithm for online guidance.



---

*Source: 5729335-2022-01-18.xml*