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 are the coefficient matric expressed as where is the -th row of the identity matrix with the dimension . Then, combining (79) with (81), we get
By using Chebyshev polynomials to approximate the integrands, the analytical solution to can be obtained as
Similar to (76), the first term of (84) can be rewritten as where
Similarly, substituting (69) and (76) into the second term of (84), it can be expressed as where
Then, combining (85) with (87), the analytical solution to is
Similar to (78) and (84), the analytical solution to can be approximated as
Substituting (89) into (90), the first term of (90) can be rewritten as where
Similar to (91), the second term of (90) can be rewritten as where
Substituting (69) and (72) into the third term of (90), it can be expressed as where
Therefore, combining (91), (93), and (95), the analytical solution to is expressed as where
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 where and are the partitioned diagonal matrices expressed as (100), respectively.
The coefficients are 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 where is the analytical solution to the trajectory with zero AOA, which is expressed as
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: where is the value of at the terminal time. Then, substituting (104) into (3), the terminal constraint functions can be converted to the functions about the AOA vector, :
Furthermore, using Chebyshev polynomial approximation, the performance index expressed in (4) can also be converted into where is a positive definite diagonal matrix whose elements along the diagonal are
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, . 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, . 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: where is normal AOA vector whose elements are the value of nominal AOA profile at interpolating points, i.e.,
is the deviation of the real AOA from the nominal AOA; is the partial derivative of the terminal constraint function with respect to the terminal state ; and is the partial derivative of the terminal state with respect to the vector of AOA at the interpolation points, , which can be expressed as where and . is the terminal value of trajectory with the nominal AOA profile:
As for the performance index (106), using the as the optimal variables, it can also be rewritten as
Considering the terminal constraints, the augmented performance index is where is the Lagrange multiplier vector. According to numerical optimization theory, the corresponding KKT conditions are
It is easy to find that the KKT conditions are a set of linear algebraic equations which can be expressed as where is the requested variables; The elements of and are defined as
Note that 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
Then, using the previous calculation process, we iteratively calculate until , 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 can be approximated by Chebyshev interpolating polynomial: where is the value of at the time , which can be obtained by
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 Figure 3 and summarized in detail as follows. (1)In the initialization, set initial parameters and select zero AOA as the initial nominal AOA (2)Calculate the analytical solutions to and the coefficients using the current states(3)Apply the method of linearization around the nominal AOA to formulate the matrices and , and calculate control deviation by solving linear algebraic equations(4)Generate the new nominal AOA, and repeat step 3 until is within the prescribed tolerance, then go to Step 5(5)Update the current guidance command by (118). Then, go to next step(6)Apply 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 step
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 to , , and . And three types of AOA curve listed in Table 2 are considered.
Firstly, the cases of three types of AOA curve with 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 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.
(a) AOA curve 1
(b) AOA curve 2
(c) AOA curve 3
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 Table 3. 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 coefficient 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.
Figures 5 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.
Table 5 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 Table 6. 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, , , and , 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.
The flight-path angle and altitude histories are depicted in Figure 7. 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 deg for fight-path angle and 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.
(a) Flight-path angle histories
(b) Altitude histories
The AOA profiles obtained by proposed guidance are presented in Figure 8. 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.
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.
Dispersion parameters used in simulations are listed in Table 8. 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.
The dispersion of terminal states are presented in Figure 10. 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 -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.
(a) Case 1
(b) Case 2
(c) Case 3
(d) Case 4
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 than deg for flight-path angle and 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.
Appendix
A. Appendix
As for the following nonlinear integral initial value problem: it is difficult to obtain the analytical solution directly. The regular perturbation technique [43, 44] is an effective method to derive its approximate analytic solution. Firstly, (A.1) needs be rewritten 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 the regular perturbation technique, can be expanded in terms of a parameter as follows: where the superscripts are used to denote the th-order solutions. Substituting (A.3) into (A.2) yields
By expanding (A.4) into the Taylor series expansion and treating both sides of the equation according to the equal coefficient of to the same power, the expanded-form dynamics can be obtained, where the zeroth-order and first-order dynamics are where is the partial derivative of to . It can be found from (A.5) that (1)the zero-order dynamics is only concerned with zero-order items. By selecting proper , the zero-order items can be solved analytically(2)The first coefficient and the second term on the right side of the first-order dynamics are functions of the zero-order items. Under the premise that the analytical solution to zero-order items have been obtained, the first-order dynamics is a linear time-varying system, which can be solve analytically using some numerical technologies
Therefore, the solution of state can be approximated as (A.6) with high accuracy.
It is noted that the accuracy of analytical solution depends on the difference between and . The more is larger than ; the higher the accuracy of the analytic solution will be, and vice versa. Therefore, the key of using regular perturbation method is how to select the appropriate functions, and , which can not only ensure the accuracy of solution but also facilitate analytically integrating.
B. Appendix
As for the integral term , let where is the value of at the -th interpolation point and is the interpolation basis function defined as where the interpolation points are the roots of Chebyshev polynomials.
is the dimensionless variable with , where and . It is obvious that is a -th order polynomial, which can be rewritten as where the coefficients are where and meet the following recurrence relation:
Therefore, the coefficient is
Then, the analytical solution to is derived as where is
Substituting (A.13) into (A.14), can be rewritten as where is the integral weight coefficient expressed as
Data Availability
The simulation data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported in part by the National Natural Science Foundation of China (NSFC) under Grant 62003019.