2.0 OPTIMAL ENDO-ATMOSPHERIC ASCENT MODEL OF LAUNCH VEHICLES

2.1 Endo-atmospheric motion mathematical model

To increase stability and accuracy of the numerical calculation, a non-dimensional motion mathematical model was adopted. Non-dimensional motion equations of a launch vehicle in the Launch Inertial Coordinate System are as follows^{(Reference Huang, Wei, Gu and Cui13)}:

\begin{matrix}{\dot{r} = V} \cr {\dot{V} = - r/{r^3} + (T - A)1}_{b} + N{1_n}\end{matrix} (1)

where r and $$V \in {R^3}$ are the geocentric distance vector and velocity vector, respectively; T represents the non-dimensional engine thrust acceleration; A and N represent the non-dimensional aerodynamic acceleration in axial and normal directions, respectively; 1 _b is the unit vector in the axial direction of the launch vehicle; 1 _n is the unit vector perpendicular to 1 _b located within the longitudinal symmetry plane of the launch vehicle.

In order to ensure structure safety, the endo-atmospheric ascent flight of the launch vehicle needs to satisfy some path constraints. The aerodynamic bending moment constraint (product of attack angle and dynamic pressure) was taken into consideration in this study

\begin{equation}S = |q\alpha | - {Q_\alpha } \le 0\end{equation} (2)

where $$q = \rho V_r^2/2$ ; Q_α is the allowable maximum aerodynamic bending moment.

Normally terminal constraints are considered as standard burnout conditions, including geocentric distance $$ r_f ^*$ , velocity $$ V_f ^*$ , orbit inclination i* and flight path angle $$\gamma _f^*$ , which are equivalent to the given semi-major axis, eccentricity, orbit inclination and true anomaly at burnout time. Terminal constraints are represented as:

$$\matrix{{{1 \over 2}r_f^T{r_f} - {1 \over 2}r_f^{*2} = 0} \hfill\cr{} \hfill\cr{1_N^T({r_f} \times {V_f}) - \left\| {{r_f} \times {V_f}} \right\|\cos{i^*} = 0} \hfill\cr{} \hfill\cr{r_f^T{V_f} - {r_f}{V_f}\sin \gamma _f^* = 0} \hfill\cr{} \hfill\cr{{1 \over 2}V_f^T{V_f} - {1 \over 2}V_f^{*2} = 0} \hfill\cr} $$ (3)

Then, the optimal endo-atmospheric ascent problem can be equivalent to finding the optimal body-axis direction 1 _b and engine burnout time t_f according to the current initial conditions, so as to make the trajectory of the launch vehicle reach the given final conditions at burnout time with the path constraints being satisfied.

2.2 Optimal control model

Generally, a double-layer iterative method is adopted to solve the optimal endo-atmospheric ascent problem of a launch vehicle^{(Reference Huang, Wei, Gu and Cui13)}. For the internal layer, an optimal control algorithm is adopted to solve the ascent trajectory problem of acquired maximum terminal energy at given flight time; whereas for the external layer, the secant method is employed to iteratively solve the flight time to meet the terminal velocity constraint. A rapid, accurate and reliable solution of optimal ascent trajectory at the internal layer is the key to realising rapid optimal trajectory generation, as well as the research focus of this paper.

Performance index is selected as:

$$J = \phi ({r_f},{V_f},{t_f}) = 1/{r_f} - V_f^2$$ (4)

In the application of the optimal control theory, the Hamiltonian function considering the path constraints is represented as:

$$H = p_r^{\rm{T}}V + p_V^{\rm{T}}\left[ { - r/{r^3} + (T - A{{)1}_b} + N{1_n}} \right] + {\lambda _{q\alpha }}S$$ (5)

where p _r and $${ p_V } \in {R^3}$ are co-state variables; λ_qα is the scalar multiplier; when $$S \le 0,{\lambda _{q\alpha }} = 0$ and when $$S = 0,{\lambda _{q\alpha }} = (\partial H/\partial \alpha )/q$ .

According to the standard first order necessary conditions for the optimal solution, differential equations of co-state variables involve the calculation of partial derivatives of vectors, and the forms are relatively complex. According to control equations, the optimal body-axis expression can be derived as^{(Reference Huang, Wei, Gu and Cui13)}:

$$1_b^* = \left( {{{\sin \alpha } \over {\sin \Phi }}} \right){1_{{p_V}}} + \left[ {{{\cos\alpha - \cos\Phi \cos(\Phi - \alpha )} \over {\mathop{\sin }\nolimits^2 \Phi }}} \right]{1_{{V_r}}}$$ (6)

where Φ is the angle between p _V and $${ V_r }$ ; $${1_{{ p_V }}}$ is the unit vector of p_V ; α is obtained by solving $$\tan (\Phi - \alpha )(T - A + {N_\alpha }) - ({A_\alpha } + N) = 0$ .

According to the transversal conditions, three terminal constraint conditions are obtained:

$$\matrix{{(V_f^{\rm{T}}{p_{{r_f}}})r_f^2 - (r_f^{\rm{T}}{p_{{V_f}}})V_f^2 + (r_f^{\rm{T}}{V_f})(V_f^2 - r_f^{\rm{T}}{p_{{r_f}}}) = 0} \hfill\cr{} \hfill\cr{V_f^{\rm{T}}{p_{{V_f}}} - V_f^2 = 0} \hfill\cr{} \hfill\cr{(H_f^{\rm{T}}{p_{{r_f}}})\left[ {H_f^{\rm{T}}({r_f} \times {1_N}) + (H_f^{\rm{T}}{p_{{V_f}}})(H_f^{\rm{T}}({V_f} \times {1_N}))} \right] = 0} \hfill\cr} $$ (7)

where, $${ H_f } = { r_f } \times { V_f }$ .

Equation (7) and the first three constraints in Equation (3) constitute six terminal boundary conditions.

Based on the above derivation, the optimal ascent problem of a launch vehicle can be converted into a TPBVP consisting of differential equations of state variables and co-state variables, initial state conditions and terminal boundary conditions. By giving definitions of $$x = ({r^{\rm {T}}},{V^{\rm {T}}}{)^{\rm {T}}},p = ( p_r ^{\rm {T}}, p_V ^{\rm {T}}{)^{\rm {T}}}$ , and $$y = ({x^{\rm {T}}},{p^{\rm {T}}}{)^{\rm {T}}} \in {R^{12}}$ , the model can be expressed as

$$\matrix{{\dot y = f(t,y)} \hfill\cr{{B_0}({y_0}) = 0} \hfill\cr{{B_f}({y_f}) = 0} \hfill\cr} $$ (8)

where B ₀ is the given six initial state conditions x (t ₀); B _f is the six terminal boundary conditions. The TPBVP can be depicted as: to find initial values p(t₀) of co-state variables that make the system $$\dot y = f( t,y )$ meet boundary conditions Bf = 0 at the terminal.

3.0 RAPID MULTI-LAYER SOLVING METHOD

The model of the TPBVP provided in Equation (8) is very complex. The finite difference method is one of the classical solving methods currently used^{(Reference Lu, Sun and Tsai9)}, which discretizes nonlinear differential equations into linear algebraic equations through central finite difference. The principle of the method is as follows:

The time interval t_f - t ₀} is divided into M subintervals with an identical length of $$h = ({ t_f }{\kern 1pt} - {\kern 1pt} {t_0})/M$ . Then differential equations in Equation (8) can be obtained by central finite difference approximation:

$${1 \over h}({y_k} - {y_{k - 1}}) = f\left( {{t_{k - 1/2}},{{{y_k} + {y_{k - 1}}} \over 2}} \right)$$ (9)

The above equations can be rewritten as the finite difference equations (FDE) is Equation (10)

$${E_k}({y_k},{y_{k - 1}}) = {y_k} - {y_{k - 1}} - hf\left( {{t_{k - 1/2}},{{{y_k} + {y_{k - 1}}} \over 2}} \right) = 0$$ (10)

where $${E_0} = [B_0^{ {T}}, B_f ^{ {T}}{]^{ {T}}}$ , and the TPBVP is then converted into a problem of solving $$Y = (y_0^{ {T}},y_1^{ {T}}... y_M ^{ {T}}{)^{ {T}}} \in {R^{12(M + 1)}}$ , the solution of 12(M + 1) nonlinear algebraic equations $$E = (E_0^{{T}},E_1^{{T}}...E_M^{ {T}}{)^{ {T}}}$ .

When the finite difference method is used to solve optimal ascent trajectory of a launch vehicle, attention should be paid to the following points:

1. 1) A homotopy technique is needed (taking the atmospheric density or reference area as homotopy parameters) to ensure reliable convergence of the algorithm. The required homotopy time increases with the number of discrete intervals increasing;

2. 2) More discrete intervals, unknown node variables and density homotopy times are needed to obtain higher solving precision, which greatly reduces the solving efficiency of the algorithm;

3. 3) Fewer discrete intervals and unknown node variables are needed to improve solving efficiency, so as to ensure real-time performance of trajectory online generation and closed-loop guidance. However, this will cause the loss of solving precision, making it impossible to realise optimal closed-loop guidance in a real sense.

In order to remedy the weakness, this paper puts forward a rapid multi-layer solving method. First, the problem was solved with fewer initial discrete intervals; then, the interpolation values, which were obtained by interpolating convergent solutions at sparse nodes into doubly discrete nodes of higher layers, were taken as initial values to carry out higher-layer iterations. The process was repeatedly performed until the solving precision met requirements. The implementation steps are as follows:

1. (1) The initial value of the higher-layer iteration is obtained by solving the TPBVP based on fewer initial discrete intervals:

1. i. The initial values of engine burnout time and layer number K are set to t _f0 and 0, respectively; the number of discrete intervals is set as M = N ₀, wherein N ₀ depends on the engine burnout time and is set to be smaller than 10.

2. ii. The analytic vacuum solution of the non-constraint optimal ascent problem is taken as the iterative initial value Y ₀;

3. iii. By the finite difference method, the optimal ascent trajectory problem of acquired maximum terminal energy with fixed final time is transformed into FDE, which is iteratively solved by a numerical algorithm, such as the traditional Newton method with a relaxation factor.

For the lower-layer iteration, the homotopy technology is necessary to obtain a full atmospheric solution, whose continuation form is only the reference area;

$$\hat{S} = \varepsilon S,0 \le \varepsilon\le 1$$ (11)

The homotopic parameter ε is initially set to 0 for the vacuum solution and gradually increases to unity for the full atmospheric solution, that is, ε = 1. The convergent solution for the current value of ε serves as the starting point of the solution for the next value of ε till ε = 1.

1. iv. Engine burnout time t_f is adjusted by the secant method to meet the terminal velocity constraint. State variables, co-state variables and t_f at fewer discrete interval nodes are obtained.

1. (2) The TPBVP is solved with the multi-layer iteration method:

1. i. End if the desirable solving precision is obtained, otherwise double the time grid density, K = K + 1;

2. ii. State variables and co-state variables at layer K obtained in step 1) are allocated to newly increased nodes at layer K + 1 through Lagrange interpolation, as shown in Fig. 1;

3. iii. The state variables, co-state variables and t_f at nodes at layer K + 1 are obtained by using the method of iii and iv in step 1), wherein the homotopy technology is not applied;

4. iv. Return to step i.

Figure 1. Diagram of multi-layer recursion.

3.1 Improved numerical solution for finite difference equation

By introducing a multi-layer solving strategy, the solving efficiency and adaptability have been improved. Moreover, two improved numerical algorithms without and with fewer Jacobian calculations are proposed for solving the finite difference equation, in order to further improve the solving efficiency from the aspect of the numerical algorithm.

The Newton iteration method with a relaxation factor is usually used to solve the FDE through combined applications of the analytical vacuum initial guess technology and homotopy technology^{(Reference Huang, Wei, Gu and Cui13)}. The update of Y is given by:

$${{Y_j} = {Y_{j - 1}} + {\sigma _j}{d_j}} 0 \lt {\sigma _j} \le 1 {{d_j} = - {{\left[ {{{\partial E({Y_{j - 1}})} \over {\partial Y}}} \right]}^{ - 1}}E({Y_{j - 1}})} $$ (12)

where, the relaxation factor σ_j is determined by the following monotone criterion:

$${\sigma _j} = \mathop{\max }\limits_{0 \le i} \{ {{\beta ^i}|{E^{\rm{T}}}[ {{Y_{j - 1}} + {\beta ^i}{d_j}]E[ {{Y_{j - 1}} + {\beta ^i}{d_j}} ] \lt {E^{\rm{T}}}({Y_{j - 1}})E({Y_{j - 1}})} 0 \lt \beta \lt 1$$ (13)

A significant disadvantage of the Newton iterative method is that the Jacobian must be calculated in each iteration, which accounts for a large proportion of numerical computation load. Therefore, improved solving methods without or with fewer Jacobian calculations should be studied.

Traditional solving algorithms require gradient calculations, such as conjugate gradient methods, spectral gradient methods, and their deformations. Therefore, the key to the derivative-free algorithm is to approximate the gradient with a gradient-free search direction. Motivated by the idea, for solving the nonlinear equations F (x) = 0, La Cruz and Raydan^{(Reference Cruz, Martinez and Raydan15)} introduced a derivative-free algorithm named the Derivative-free Spectral Algorithm for Nonlinear Equations (DF-SANE). This algorithm extends the spectral gradient approach and uses the residual $${ d_k } = \pm F\left( {{ x_j }} \right)$ as search directions. Moreover, a merit function $$f\left( x \right) = {1 \over 2}{\left[ {F\left( x \right)} \right]^2}$ with the property that $$f\left( x \right) = 0 \ Leftrightarrow F\left( x \right) = 0$ is introduced. The first trial point at each iteration is $${ Y_j } = {Y_{j - 1}} + {\sigma _j}E({Y_{j - 1}})$ where σ_j is a spectral coefficient.

$$\matrix{{{\sigma _j} = {{s{'_{j - 1}}{s_{j - 1}}} \over {s{'_{j - 1}}{y_{j - 1}}}}} \cr{{s_{j - 1}} = {Y_j} - {Y_{j - 1}},{y_{j - 1}} = E({Y_j}) - E({Y_{j - 1}})} \cr} $$ (14)

The spectral coefficient σ_j is an appropriate Rayleigh quotient with respect to a secant approximation of the Jacobian.

$${{\left\langle{\left[ {\int_0^1 J\left( {{x_j} + t{s_j}} \right){\rm{d}}t} \right]{s_j},{s_j}} \right\rangle} \over {\left\langle{{s_j},{s_j}} \right\rangle}}$$ (15)

which makes it attractive for the numerical solution of (10).

Since $${ d_k } = \pm F\left( {{ x_k }} \right)$ is not necessarily a descent direction, the global convergence is guaranteed by an effective line search strategy. The majority of line search methods require a decrease in the objective function f _x at each iteration, which means that the corresponding sequence of function values monotonically decreases. This may result in a slow convergence. On the other hand, the nonmonotone line search strategies could converge faster and avoid the local convergence problem.

The well-known nonmonotone line-search technique was introduced by Grippo, Lampariello and Lucidi(GLL)^{(Reference Grippo, Lampariello and Lucidi16)}. Li and Fukushima^{(Reference Li and Fukushima17)} presented the LF method that avoided the necessity of descent directions to guarantee that each iteration was well defined.

The GLL condition can be written as follows:

$$f\left( {{x_k} + {\alpha _k}{d_k}} \right) \le {\bar f_k} + \gamma {\alpha _k}\nablaf{\left( {{x_k}} \right)^T}{d_k},{\bar f_k} = \mathop{\max }\limits_{0 \le j \lt M - 1} f\left( {{x_{k - j}}} \right)$$ (16)

The LF condition can be written as follows:

$$\left\| {F\left( {{x_k} + {\alpha _k}{d_k}} \right)} \right\| \le \left( {1 + {\eta _k}} \right)\left\| {F\left( {{x_k}} \right)} \right\| - \gamma \alpha _k^2\left\| {{d_k}} \right\|_2^2$$ (17)

The GLL strategy tolerates an increase in the objective function (nonmonotonic behavior), but requires the gradient and descent directions to guarantee global convergence. The LF strategy accepts nondescent directions but the nonmonotone behavior is insufficient.

By taking into account the advantages of both schemes, the new nonmonotone line search condition for DF-SANE can be written as:

$$f\left( {{x_{k + 1}}} \right){\kern 1pt} \le {\kern 1pt} {\kern 1pt} {\bar f_k} + {\eta _k} - \gamma \alpha _k^2{\kern 1pt} f\left( {{x_k}} \right),\;{\bar f_k} = \mathop{\max }\limits_{0 \le j \lt M - 1} f\left( {{x_{k - j}}} \right)$$ (18)

The GLL term $$\ mathop {\max }\limits_{0 \le j \lt M - 1} f\left( {{x_{k - j}}} \right)$ guarantees the sufficient nonmonotonicity of $$f\left( {{ x_k }} \right)$ ; η_k is responsible for the result that all the iterations are well defined; and the term $$ - \gamma \alpha _k^2f\left( {{ x_k }} \right)$ is responsible for improving the global convergence property.

Step 0. Set k = 0.

Step 1. Choose such that σ_k $$\left| {{\sigma _k}} \right| \in \left[ {{\sigma _{\min }},{\sigma _{\max }}} \right]$ (the spectral coefficient)

$${\bar f_k} = \max \left\{ {f\left( {{x_k}} \right), \cdots,{\kern 1pt} f\left( {{x_{\max \left\{ {0,k - M + 1} \right\}}}} \right)} \right\}$$

Set $$d = - {\sigma _k}F\left( {{ x_k }} \right), {\alpha _ + } = 1,{\alpha _ - } = - 1$

Step 2. If $$f\left( {{ x_k } + {\alpha _ + }d} \right) \le {\bar f_k } + {\eta _k} - \gamma \alpha _ + ^2f\left( {{ x_k }} \right)$ then define $${ d_k } = d,{\alpha _k} = {\alpha _ + },{x_{k + 1}} = { x_k } + {\alpha _k}{ d_k }$

Else if $$f\left( {{ x_k } - {\alpha _ - }d} \right) \le {\bar f_k } + {\eta _k} - \gamma \alpha _ - ^2f\left( {{ x_k }} \right)$ then define $${ d_k } = - d,{\alpha _k} = {\alpha _ - },{x_{k + 1}} = { x_k } + {\alpha _k}{ d_k }$

Else Choose $${\alpha _{ + new}} \in \left[ {{\tau _{\min }}{\alpha _ + },{\tau _{\max }}{\alpha _ + }} \right],{\alpha _{ - new}} \in \left[ {{\tau _{\min }}{\alpha _ - },{\tau _{\max }}{\alpha _ - }} \right]$ .

$${\alpha _{ + new\,}} = \left\{ {\matrix{ {{\tau _{\min }}{\alpha _{ + \,\,}}\,if\,{\alpha _t}\, \lt \,\,{\tau _{\min }}{\alpha _ + }} \hfill \cr {{\tau _{\max }}{\alpha _{ + \,}}\,if\,{\alpha _t}\,\, \gt \,\,{\tau _{\max }}{\alpha _ + }} \hfill \cr {{\alpha _{t\,\,\,\,\,\,\,\,\,}}\,\,\,\,otherwise} \hfill \cr } } \right.$$ (19)

where

$${\alpha _t} = {{\alpha _ + ^2f\left( {{x_k}} \right)} \over {f\left( {{x_k} + {\alpha _ + }d} \right) + \left( {2{\alpha _ + } - 1} \right)f\left( {{x_k}} \right)}}$$

Replace $${\alpha _{ + new}} = {\alpha _ + },{\alpha _{ - new}} = {\alpha _ - }$ and go to Step 2.

Step 3. If $$\left\| {F\left( {{x_{k + 1}}} \right)} \right\| \lt \ varepsilon $ , terminate the execution of the algorithm. Else, k = k + 1 and go to Step 1.

However, since the search direction is not the descent direction, the DF-SANE method may be trapped in local convergence problems. Therefore, an improved solving algorithm with fewer Jacobian calculations is also proposed in this paper. In order to reduce the calculation load of solving the FDE in Equation (10), a modified Newton method with a relaxation factor in combination with the inverse Broyden quasi-Newton method is applied, denoted as “MN-IBQ”. Through combining the advantages of the Newton method and the simple Newton method, the modified Newton method is capable of performing the Jacobian calculation once every two steps, which not only preserves the fast convergence property of the Newton method, but also reduces the number of Jacobian calculations. The formula for the modified Newton method with a relaxation factor is as follows:

The search direction d _j in the j-th iteration is determined by solving the linear algebraic equations

$$\matrix{{\matrix{{{Y_{j - 1,0}} = {Y_{j - 1}}} \hfill\cr{{Y_{j - 1,1}} = {Y_{j - 1}} - {\sigma _{j,1}}{{\left[ {{{\partial E({Y_{j - 1}})} \over {\partial Y}}} \right]}^{ - 1}}E({Y_{j - 1}})} \hfill\cr{{Y_j} = {Y_{j - 1,1}} - {\sigma _{j,2}}{{\left[ {{{\partial E({Y_{j - 1}})} \over {\partial Y}}} \right]}^{ - 1}}E({Y_{j - 1,1}})} \hfill\cr} } & {j = 1,2, \cdots} & {} \cr} $$ (20)

Such search step size ensures that the sequence $$\left\{ {\left\| {E({Y_{j - 1}})} \right\|} \right\}$ is monotonically decreasing. Convergence is achieved when $$\left\| {E({Y_{j - 1}})} \right\|$ is not greater than a preselected tolerance.

Next, the step size parameters $${\sigma _{j,1}},{\sigma _{j,2}}$ are determined by the following criterion:

$$\begin{array}{l} {\begin{array}{cc}{\sigma _{j,1} =\mathop{\max }\limits_{0\le i} \left\{\beta ^{i} |\textit{\textbf{E}}^{{\rm T}} \left[\textit{\textbf{Y}}_{j-1} +\beta ^{i} \textit{\textbf{d}}_{j,1} \right]\textit{\textbf{E}}\left[\textit{\textbf{Y}}_{j-1} +\beta ^{i} \textit{\textbf{d}}_{j,1} \right]<\textit{\textbf{E}}^{{\rm T}} (\textit{\textbf{Y}}_{j-1} )\textit{\textbf{E}}^{T} (\textit{\textbf{Y}}_{j-1} )\right\}\!,} & {0<\beta <1} \end{array}} \\[6pt] {\begin{array}{cc} {\sigma _{j,2} =\mathop{\max }\limits_{0\le i} \left\{\beta ^{i} |\textit{\textbf{E}}^{{\rm T}} \left[\textit{\textbf{Y}}_{j-1,1} +\beta ^{i} \textit{\textbf{d}}_{j,2} \right]\textit{\textbf{E}}\left[\textit{\textbf{Y}}_{j-1,1} +\beta ^{i} \textit{\textbf{d}}_{j,2} \right]<\textit{\textbf{E}}^{{\rm T}} (\textit{\textbf{Y}}_{j-1,1} )\textit{\textbf{E}}^{T} (\textit{\textbf{Y}}_{j-1,1} )\right\},} & {0<\beta <1} \end{array}} \\[6pt] {\textit{\textbf{d}}_{j,1} =\left[\frac{\partial \textit{\textbf{E}}(\textit{\textbf{Y}}_{j-1} )}{\partial \textit{\textbf{Y}}} \right]^{-1} \textit{\textbf{E}}(\textit{\textbf{Y}}_{j-1} ),\textit{\textbf{d}}_{j,2} =\left[\frac{\partial \textit{\textbf{E}}(\textit{\textbf{Y}}_{j-1} )}{\partial \textit{\textbf{Y}}} \right]^{-1} \textit{\textbf{E}}(\textit{\textbf{Y}}_{j-1,1} )}\hspace*{115pt}.\kern1.5pt.\kern1.5pt. \\[-15pt]\end{array}\end{align}$$ (21)

To further improve the solving efficiency, the inverse Broyden quasi-Newton method^{(Reference Kvaalen18)} that requires less computation load is introduced.

The inverse Broyden quasi-Newton method can provide a superlinear convergence rate. Meanwhile, it does not require the Jacobi matrix recalculation or matrix inversion operation, but only needs to determine the H-correction matrix. The iterative formula of the inverse Broyden quasi-Newton method is as follows:

$$\matrix{{{Y_j} = {Y_{j - 1}} - {H_{j - 1}}E({Y_{j - 1}})} \hfill\cr{{H_j} = {H_{j - 1}} + \Delta {H_{j - 1}}} \hfill\cr{\Delta {H_{j - 1}} = ({r_j} - {H_{j - 1}}{y_j}){{({r_j}){H_{j - 1}}} \over {({r_j}){H_{j - 1}}{y_j}}}} \hfill\cr} $$ (22)

where, $${ r_j } = { Y_j } - {Y_{j - 1}},{ y_j } = E({ Y_j }) - E({Y_{j - 1}})$ .

When applying Equation (22), we need good initial values and H matrix approximations in order to obtain better solutions quickly. Therefore, a hybrid solution can be used, in which the modified Newton iteration method with a relaxation factor is applied firstly until the iteration bias converges to a certain accuracy, and then the inverse Broyden quasi-Newton method is used until the final convergence accuracy is satisfied.

4.0 SIMULATION RESULTS AND DISCUSSION

Simulation verifications were carried out based on a suborbital reusable launch vehicle. The launch vehicle vertically rose for 5.0s after being launched from the ground, and then turned with exponential attack angle for 25.0s. An optimal regeneration of the ascent trajectory after 30.0s was needed to correct the bias caused by the wind disturbance, thrust deviation and aerodynamic parameter deviation. The aerodynamic bending moment was set to be ≤ 0.7 kPa•rad, and the required terminal constraints included a height of 60.0km, a velocity of 2700.0m/s, a flight path angle of 25.3° and an orbit inclination of 15.0°. All the algorithms involved were coded in C and executed in an industrial PC with a 1.0 GHz processor.

Firstly, in order to demonstrate the effectiveness of the multi-layer solving strategy, the multi-layer method was compared with the traditional single-layer method by using the traditional numerical solving algorithm. Secondly, on the basis of the multi-layer solving strategy, to validate the effectiveness of the improved numerical solving algorithm, two improved numerical algorithms were further compared with the traditional one. Thirdly, in order to verify the adaptability of the proposed multi-layer algorithm to thrust loss, the maximum fault degrees corresponding to different fault times were given. Finally, for a typical maximum thrust loss condition, a comparison simulation between the single-layer algorithm and the multi-layer method with improved numerical algorithms was conducted.

To investigate the advantages of the multi-layer solving method proposed over the traditional single-layer direct solving method, a comparison study was carried out for the nominal ascent guidance with different interval numbers N = 10, 20, 40. The initial number of discrete intervals N₀ was set to 5. For the two solving methods, the FDE was both solved by the traditional Newton method with a relaxation factor.

According to the optimal guidance command (burnout time and attitude angle) acquired online, the integral state variables of each solving method can be firstly obtained by integrating the dynamic equations in Equation (1), then the solving precision can be checked in terms of deviation between the integral state variables and optimised state variables. Table 1 shows the statistics of the homotopy times, positional accuracy, velocity accuracy and consuming time of CPU of each method.From Table 1, the following conclusions can be drawn.

1. (1) When the number of discrete intervals is the same, the multi-layer solution and traditional direction solution have consistent positional accuracy and velocity accuracy, while the solving efficiency of the multi-layer method is more than 3 times as high as that of the traditional direct solution.

2. (2) With the increase of the number of discrete intervals, the positional accuracy and velocity accuracy of optimal ascent trajectory gradually increase, and the homotopy time of the direct solution increases gradually, while that of the multi-layer solution remains unchanged.

3. (3) With the increase of the number of discrete intervals, the consuming time of CPU for the multi-layer solution increases gradually, while the consuming time of CPU for the direct solution increases rapidly.

Table 1 Comparison of simulation results between multi-layer solving method and direct solving method under normal conditions

The multi-layer solving method proposed in this paper has substantial improvements compared with the traditional one. The bigger the number of discrete intervals is, the higher the improved efficiency is. Figures 2–3 respectively show the height and velocity deviation curves when N is set to different values. Both height and velocity deviations reach the maximum at the terminal time. Figures 4–5 respectively show the angle of attack and pitch angle curves calculated by the multi-layer solving method.

Figure 2. Time-dependent variation of height deviation for normal condition.

Figure 3. Time-dependent variation of velocity deviation for normal condition.

Figure 4. Time-dependent variation of angle-of-attack for normal condition.

Figure 5. Time-dependent variation of pitch angle for normal condition.

On the basis of the multi-layer solving strategy, to validate the effectiveness of the improved numerical algorithm for solving FDE, a comparative analysis of two numerical solutions is conducted, namely the traditional Newton method with a relaxation factor, and the DF-SANE, which are denoted as “-Traditional” and “-DF-SANE”, respectively. DF-SANE was implemented with the following parameters: $${\sigma _{\min }} = 1e - 10$ , σ_max = 1e10, $${\tau _{\min }} = 0.05$ , $${\tau _{\max }} = 0.5$ , $${\eta _k} = \left\| {F\left( {{x_0}} \right)} \right\|/{\left( {1 + k} \right)^2}$ , γ = 1e3, M = 20, ε = 1e − 6.

The simulation results show that DF-SANE converged to the local minimum and failed to converge to the given ε. The residual curve and search step curve are as follows.

It can be seen from Fig. 6 that the residuals of the equations show obvious nonmonotonicity. The residuals gradually converge from 0.0014 on the overall trend. However, after converging to the local minimum 2.4e-4, the residuals remain stable instead of decreasing, thus failing to converge to the given ε. As can be seen from Fig. 7, after the number of iterations is greater than 160, the search step size is already less than 0.005. It would seem that the failure of DF-SANE is essentially due to the fact that the search direction in DF-SANE is not the descent direction and the non-monotone linear search method has insufficient global convergence.

Figure 6. Iteration number-dependent variation of residual for normal condition.

Figure 7. Iteration number-dependent variation of step size for normal condition.

Since the DF-SANE method without Jacobian calculations is trapped in the local convergence problem, the improved solving algorithm with fewer Jacobian calculations is proposed in this paper. It is the combination of the modified Newton method with a relaxation factor and the inverse Broyden quasi-Newton method, and denoted as “- MN-IBQ”. Meanwhile, the traditional Newton method with a relaxation factor is denoted as “-Traditional”. A comparison between MN-IBQ and Traditional is conducted. Table 2 shows the statistics of the homotopy time, positional accuracy, velocity accuracy and consuming time of CPU of each method.

Table 2 Comparison of simulation results between traditional and MN-IBQ numerical algorithms for multi-layer solving method under normal conditions

According to Table 2, the following conclusions can be drawn. Through the combination of the modified Newton method with a relaxation factor and the inverse Broyden quasi-Newton method, the solving rate can be increased by approximately 15 percent compared with the traditional Newton method with a relaxation factor.

To verify the adaptability of the proposed multi-layer solving method to thrust loss, a series of case studies with different engine failure times and degrees were carried out. Assuming that the engine fails at different times, the guidance system will perform an optimal trajectory replanning immediately after the failure is detected. The number of the discrete intervals N₀ was set to 40. The accuracy requirements included the position deviation $$\Delta r \le 20m$ and the velocity deviation $$\Delta v \le 0.5m/s$ . The results of the maximum thrust loss that the multi-layer solving method could tolerate were given, where the convergent solution would not be obtained in the homotopy process if the thrust loss was a little bigger. Corresponding deviations are shown in the following table.

From Table 3, the following conclusions can be drawn.

1. 1) The earlier the engine failure occurs, the smaller the maximum allowable thrust loss is;

2. 2) The allowable maximum thrust loss increases approximately linearly with the engine failure time.

Table 3 Results of multi-layer solving method under maximum allowable thrust deviations for different failure time

This experiment analyzed a typical worse situation in which a serious engine failure happened 70 seconds after the vehicle was launched and led to a 25% reduction in thrust. Table 4 shows the statistics of the homotopy time, positional accuracy, velocity accuracy and consuming time of CPU of each method for the case of 25% thrust reduction.

Table 4 Comparison between multi-layer solving method and direct solving method under engine failure conditions

Figures 8–9 respectively show height and velocity deviation curves when N is set to different values. Both height and velocity deviations reach the maximum at the terminal time. Figures 10–11 respectively show angle of attack and pitch angle curves calculated by the proposed multi-layer solving method.

Figure 8. Time-dependent variation of height deviation for engine failure.

Figure 9. Time-dependent variation of velocity deviation for engine failure.

Figure 10. Time-dependent variation of angle of attack for engine failure.

Figure 11. Time-dependent variation of pitch angle for engine failure.

From Table 4, conclusions can be obtained as follows:

1. (1) In engine failure situations, the traditional direct solving method can only converge when the node number N is less than 20, and its positional and velocity accuracies cannot meet the requirements;

2. (2) The proposed multi-layer solving method has strong robustness against the engine failure situation where the convergent solution cannot be obtained by the traditional direct solving method. That is, for N = 20 and N = 40, the convergent solution cannot be obtained by the direct solving method, but can be rapidly solved by the proposed multi-layer solving method;

3. (3) For the multi-layer solving method, the positional accuracy and velocity accuracy of optimal ascent trajectory gradually increase with the increase of the number of discrete intervals. N = 40 is enough to obtain the required accuracy $$\Delta r \le 20m$ ;

4. (4) Through the combination of the modified Newton method with a relaxation factor and the inverse Broyden quasi-Newton method, the solving rate can be further increased by approximately 15% compared with the traditional Newton method with a relaxation factor. The consuming time of CPU is 2.34 seconds, which can meet the requirement of the online computation speed.