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.
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: Section 2 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 where , , , , and are the mass, velocity, flight-path angle, altitude, and AOA, respectively; is the distance from the center of Earth to the vehicle, where is the average radius of the Earth, and is the gravitational acceleration; is the thrust of the vehicle; and are the aerodynamic drag and lift expressed as where is the atmospheric density and is the reference area. The terms, and , 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 where and are the ascent phase terminal flight-path angle and altitude and and 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.
It is noted that the denominator, , is used to shape the AOA profile and makes the terminal AOA zero if is selected.
2.3. Vehicle Description
For rocket-powered launch vehicle, its thrust can be formulated as where is specific impulse and is the mass flow which can be treated as a constant for the solid launch vehicle; is the area of nozzle; and 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 where , , and are all related to Mach number, which can be fitted as the polynomial of Mach number, i.e.,
In the solving of the analytical solutions, these coefficients are regarded as the constants, , , and , 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.
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, , where is the sea-level atmospheric density and , is adopted for the solving of the analytical solutions.
In order to avoid confusion, the variables corresponding to zero AOA are expressed with subscript “.” For the case of zero AOA, it is obvious that and ; then, the dynamics is simplified to
For simplifying the derivation, a variable related to is introduced.
Replacing with as the state variable, the equation of motion presented in (8) can be rewritten as
The increment of the velocity, flight-path angle, and altitude with nonzero AOA relative to that with zero AOA are defined as , , and , 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.
Then, the dynamics equations of increments are established as
Because , the altitude increment on the value of 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 can be divided into the following two parts: where is the major axial force, which is the sum of thrust, average atmospheric drag, and average axial component of gravity, and is the perturbation axial force: where ,, , , and 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 can be written as follows: where is the constant coefficient, which can be obtained by the initial and terminal values of and :
In the boost phase, the flight-path angle varies slightly; therefore, the linear term is much larger than the rest terms. Therefore, according to the regular perturbation method presented in Appendix A, (10) can be modified as where is a small parameter for the subsequent regular perturbation order and is a constant, which is set to be equal to . According to (A.5), the zeroth-order dynamics is
It is noted that, in order to decouple the differential equations of and , the terms related to are expanded at initial point since the change of is much smaller than its initial value, . The first-order dynamics is where is the zeroth-order solution to flight-path angle, which can be derived from .
It should be noted that the terms related to in the differential equations of are ignored due to . Additionally, according to the boundary conditions, the initial values of zeroth-order and first-order states are
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 , , and can be derived by integrating these three differential equations in sequence.
Integrating the first equation of (18) from to , the analytical solution of at time can be easily obtained as where and . Obviously, the differential equations of other variables are all related to ; thus, the analytical solutions of these variables are the functions of . For the sake of simplifying the expression, a dimensionless variable, , is employed, which is regarded as an independent variable in the subsequent derivation. From (22), the expression of is where is initial value of . Substituting (22) into the first equation of (18), the differential equation of can be expressed as a simple function.
Dividing the second and third equations of (18) by (24) yields where and . Integrating (25) from to , can be obtained as where is
is the integral term expressed as
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 where is the value of at the -th interpolation point ; is the interpolation basis function; and is the dimensionless variable with , where and . Then, can be approximated as where 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.
In order to derive the analytical solution to , an integral related to is introduced.
Using the formula of partial integration, the analytical solution to is expressed as where is the integral shown below:
It is obvious that the analytical solution to is similar to that of . Substituting (25) into (34), can be expressed as
Then, integrating (26) from to , the analytical solution to is where
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, , as the dependent variable, the first-order dynamics can be formulated as where and are the dimensionless variables and ,, and are the constant coefficients.
The analytical solution to can be derived by integrating the first equation of (38). But the term, which is the integration of , cannot be solved analytically because is not the linear function of . Similar to (30), let where . As shown in Figure 2, the approximation polynomials are accurate enough in comparison with the original functions . Then, the analytical solution to is where according to (33), the expression of is
By integrating the second equation of (38), the analytical solution to can be obtained as where , , and are all the functions of , which are expressed as where the expression of , , and 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 where the coefficient is
Then, the analytical solutions to the two integration terms are where
By the formula of partial integration, the analytical solution to and can be derived as where the expressions of and are
Then, the analytical solution to is
As for the analytical solution to , it is the integral of from to . It is noted that all the variables in , such as , , , and , are expressed as the function of . Therefore, we can use Chebyshev interpolation polynomial to approximate the integrand. Let
Then, the solution to can be obtained as
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, , the differential equations of trajectory increments should be formulated as where , , and are all the functions of , which are expressed as
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 time and , which can be easily expressed as
In the differential equation of , it is easy to find that the term, , is much larger than the sum of other terms, , which is the dominant term. Similarly, the term, , is the dominant term of the differential equation of . As for the differential equation of , due to and , it is obvious that the term, , is the dominant term. Therefore, the dynamics can be modified as
According to the regular perturbation method presented in Appendix B, the zeroth-order and first-order dynamics are
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 , , and can be obtained as
Although the integral terms in the above formulas are quite complex, the integrands are related only to the independent variable, . Here, Chebyshev polynomials are also used to approximate the integrands. Let
Substituting (62) into (59), the analytical solution to the zeroth-order flight-path angle increment can be obtained as where is the vector consisting of the value of AOA at interpolating points, i.e.,
is the coefficient matrix expressed as where
is the diagonal matrix whose diagonal elements are the value of at interpolating points, i.e.,
Similarly, , , and 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 where is the vector consisting of the square of AOA at interpolating points, i.e.,
and are the coefficient matrices expressed as where is
Substituting (65) and (66) into (61), the analytical solution to the zeroth-order altitude increment can be obtained as where
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 is
Similar to (76), the first term of (78) can be rewritten as where
Substituting (69) and (76) into the second term of (78), it can be expressed as where and