1. Introduction
With the increasing demand for spacecraft autonomy, not only the accuracy and optimality of the result trajectory but also the stability and efficiency of the solving algorithm are of crucial importance for real-time guidance. For the hypersonic entry problem, trajectory optimization is a popular method due to its extensive application prospect, yet it is still very challenging because of the highly nonlinear and complicated dynamics and constraints involved [
1,
2,
3,
4,
5].
For the trajectory optimization problem, the existing methods can be classified as indirect methods and direct methods [
6]. Based on Pontryagin’s maximum principle [
7], indirect methods derive the necessary conditions of optimality and solve a two-point boundary value problem (TPBVP) to obtain the result trajectory [
8,
9]. Due to the inherent nonlinearity in dynamics of the hypersonic entry problem, the solution of TPBVP presents significant challenges and falls short of real-time requirements [
10]. Despite theoretical guarantees on the optimality of the solution, practical considerations limit its feasibility. The direct methods convert the original problem into an approximated finite parameter optimization problem and use a nonlinear programming (NLP) algorithm to solve it [
11,
12,
13]. Although NLP-based methods have been successful in many applications, the time-consuming solution process and the lack of convergence guarantee are major challenges [
10]. Furthermore, an appropriate initial guess should be selected to achieve a high-quality solution, especially for the hypersonic entry problem.
In recent years, the application of convex optimization in trajectory optimization has shown great potential in many problem, including powered descent guidance [
14] and spacecraft rendezvous and proximity operations [
15]. If a problem can be formulated in a convex form, it can be solved in polynomial time with a strong convergence guarantee while obtaining the global optimal solution [
16,
17]. There are many state-of-the-art solvers based on the interior-point method (IPM) [
18], such as GUROBI [
19], MOSEK [
20], and ECOS [
21], to solve convex optimization problems efficiently. However, the majority of the trajectory optimization problems are non-convex. To deal with the nonlinearity in dynamics and constraints, the sequential convex programming (SCP) technique is developed to solve a sequence of convex subproblems to approximate the original problem [
22]. The sub-problems are formulated by linearization about a reference trajectory and subsequently solved via iterative refinement of the reference trajectory until convergence of the solution is achieved. In order to ensure the convex approximation is accurate, a trust region is imposed to make sure that the solution is not far from the reference trajectory. There are both hard and soft methods to address the trust region [
23,
24]. Recently, this SCP technique has been applied to hypersonic re-entry problems and has shown its effectiveness to obtain high-accuracy solutions. Liu formulated the entry problem as a second-order cone programming (SOCP) problem and applies the successive convexification method to solve it [
25]. Wang and Grant proposed an improved SCP algorithm and introduced a new control input for the entry problem [
26,
27]. Wang and Cui developed a rapid trajectory optimization algorithm with the pseudospectral method [
28]. However, multiple hundreds of temporal nodes are usually required to maintain the accuracy due to the nonlinearity of the hypersonic entry problem and the long flight duration. Achieving a balance between real-time performance and accuracy represents a significant challenge, as reducing the number of discrete points may result in a loss of precision.
To deal with this issue, Kamath and Açıkmeşe et al. [
29,
30,
31] propose an inverse-free precise discretization based on first-order hold (FOH) discretization. Considering the consistency of the original non-convex dynamics with the reference trajectory and addition of the
stitching condition, this discretization would guarantee high accuracy with few temporal nodes and has been effectively applied to various problems, including powered descent guidance [
29], multi-phase rocket landing [
30], and hypersonic entry guidance [
31]. In [
31], the amount of the temporal nodes is only 40 to achieve the commensurate accuracy in [
27], which necessitated an excess of 200 nodes. Nevertheless, for uniform discretization, a reduced number of temporal nodes would generate a sparse time grid, which may lead to constraint violation between the temporal nodes since the constraints are only imposed at discrete points in the SCP subproblems, resulting in a new issue. In [
30], an non-uniform discretization with additional time interval dilation variables is introduced in multi-phase rocket landing to adaptively decide the turning points of different phases. A similar idea is extended to the hypersonic entry problem and the penalized trust region (PTR) algorithm, a soft trust region method, is used to construct the SCP process [
31]. However, the resulting trajectory would still experience constraint violation. In addition, in our experiment for the hypersonic entry problem, both hard and soft trust region methods with additional time interval variables showed worse convergence compared with those with uniform discretization. One of the reasons is that the dynamics are time-sensitive, and the other is that an effective search direction should be given to achieve a purposeful distribution of discrete points.
In this paper, we propose a novel adaptive non-uniform discretization method to handle the above issues.
An inverse-free precise discretization is adopted to obtain high accuracy with few temporal nodes for real-time performance.
An adaptive non-uniform discretizaiton is proposed to construct the SCP subproblem with additional penalty terms. This would give the solver a search direction to cluster the temporal nodes more purposefully and the propagated trajectory would satisfy all the path constraints as a result, which is the main contribution of this paper.
The validity of proposed method is substantiated through a numerical experiment compared with other SCP methods.
This paper is organized as follows.
Section 2 presents the model of the hypersonic entry trajectory optimization problem including the dynamics and constraints.
Section 3 introduces the adaptive non-uniform discretization and constructs the SCP subproblem and iteration process. Simulation and results analysis are shown in
Section 4. The conclusions are summarized in
Section 5.
3. Improved SCP Method with Adaptive Non-Uniform Discretization
In this section, we introduce the improved SCP algorithm using adaptive non-uniform discretization for the hypersonic entry problem. In the interest of completeness, we provide a brief introduction to the non-uniform scheme and precise discretization technique, both of which, as in [
31], are utilized in the proposed method. Further details will be presented subsequently. In order to seek an appropriate search direction and achieve a purposeful distribution of temporal nodes, additional penalty terms with respect to the nonlinear term of path constraints and the distance term from the trajectory to the no-fly zone center are considered in the SCP subproblem construction.
3.1. Time Interval Dilation
To introduce the non-uniform discretization, we consider the original nonlinear dynamics in the sub-interval
,
where
, and define an affine map to normalize the orginal time interval (may not be equal) to a fixed interval,
:
where
is the length of the
th time interval and can be referred to as the time interval dilation [
30]. So far, the dynamics Equation (
8) can be rewritten with respect to the normalized time
:
By treating as additional decision variables, the solver is allowed to decide the adaptive time grids rather than a uniform temporal grid.
What is more, some exact constraints should be added in practical implementation to ensure physical meaning:
to ensure the time order
and adjacent temporal nodes are not far away, and
to ensure the fixed flight time.
3.2. Convexification and Discretization
A convex approximation of the dynamic (
10) can be obtained by the first-order Tylor expansion about a reference trajectory
. The approximate equation is a linear time-varing (LTV) system as follows:
where
are the Jacobians of the dynamics with state control and time dilation, respectively.
For discretization, a precise inverse-free discretization technique based on first-order hold (FOH) is adopted [
29,
30,
31]. In the FOH case, the control input signal is considered as a piecewise affine function; thus, the control variables are only defined at the discrete time nodes and the control signal in the sub-interval can be parameterized as follows:
where
and
as given in (
9). Thus, The LTV dynamics (
13) can be easily rewritten with respect to the deviations from the reference trajectory.
where
denotes the reference quantity,
denotes the deviations from the reference trajectory, i.e.,
and
, and the coefficient matrixes
are the same as (
14). It can be considered that the reference trajectory in the sub-interval
is in accordance with the original dynamics (
10), rather than the convex approximation (
13), like the typical FOH discretization in [
32].
According to the knowledge of the linear system [
33], the unique solution of (
16) for
and
is
where
where
is called the state transition matrix (STM) with the following properties:
,
, and
.
In order to eliminate the inversion operation to avoid numerical problems, the
in (
18) have the closer forms, as shown in Thereom 1, which is not proven in [
29,
30,
31].
Theorem 1. The coefficient matrixes (
18)
of the linear time-varying system (
17)
have the inverse-free form: : Proof. Choosing the
as an example, then taking the derivative and invoking the chain rule yields
The first and second equal signs come from the properties of STM, while the last one is a simplification of the original form of
from (
18). With
, the inverse-free form of
is obtained as shown in (
19).
and
can be obtained by the same process. □
For simplicity, we define
and
as 0 and 1, respectively, which denote that
t is in the sub-interval
. Then evaluating the LTV system (
13) at
, we obtain
Since the reference trajectory may not satisfy the original dynamics (
1) in the sub-interval
, Equation (
22) would give
trajectory segments, which makes a discontinuity occur in the temporal time nodes
. The
stitching condition is introduced to obtain a continuous trajectory over the time horizon, as shown in
Figure 1.
As a result, the discretized dynamics in terms of the deviations are as follows:
where
,
and
, which denotes the integration of the original dynamic in
from the reference trajectory.
Thus the discretized dynamics with respect to absolute variables are recovered from Equation (
23):
With the idea of consistency of the original dynamics of the reference trajectory in sub-intervals and the addition of the stitching condition, the accuracy of the solution would be improved even over a sparse time grid.
3.3. Additional Penalty Terms in SCP Subproblem
In our experiments for the non-uniform scheme, it is observed that regardless of a hard trust region method with additional fixed constraints or a soft trust region method by augmenting the objective function with penalty terms, the method may not converge or converge very slowly, while the resulting propagated trajectory may violate the path constraints between temporal nodes as well.
As mentioned above, one reason is that an effective search direction should be given. Inspired by the PTR algorithm [
24], we augment the objective function with additional penalty terms with respect to the nonlinear term of path constraints and the distance term of the no-fly zone, which would give a more purposeful direction to distribute the temporal nodes. As the path constraints exhibit high levels of nonlinearity, the logarithm transformation is used to mitigate this issue.
3.3.1. Log-Tranforms of Path Constraints
Consider the typical path constraints (
5) of the hypersonic entry problem: heat rate, dynamic pressure, and normal load. Since the Tylor expansion of the original path constraints (
5) would obtain complicated nonlinear terms, we take the logarithm transformation of both sides of (
5):
Due to the monotonicity of logarithmic transformations, Equation (
25) is equivalent to (
5), while the transformed constraints are linear to orbital radius
r and only nonlinear to velocity
V. Thus, the Tylor expansion of the transformed constraints has a simpler nonlinear term
than that of the original constraints.
Consider the second-order Tylor series expansion of
:
Then we replace
in (
25) with the primary term
as the linearization of the transformed constraints and augment the objective function with the nonlinear term
, since the quadratic term is convex. With the variable time-step scheme and additional nonlinear penalty term, this would give the optimizer a target or direction to incentivize the temporal nodes to cluster around the highly nonlinear region.
3.3.2. Distance Penalty Term of No-Fly Zone
The propagated trajectory between temporal nodes could be within the no-fly zone because of the sparse time interval, as shown in
Figure 2. In order to address the above issue, prior work is to set a dense time grid around the no-fly zone in advance. Instead of that, our method can allow discrete points to adaptively cluster around the no-fly zone during iteration, which would take full advantage of the non-uniform scheme.
The nonlinear no-fly zone constraints (
6) can be linearized with first-order Tylor expansion:
where
.
We penalize the distance from the trajectory to the center of the no-fly zone
in the objective function, which would give the solver a search direction with physical significance. With the linear no-fly zone constraints (
27), the discrete points would tend to cluster around the no-fly zone and disperse beyond it.
Note that the penalty term of the no-fly zone should be restricted to several temporal nodes to prevent all nodes from gathering around the no-fly zone.
3.4. Discrete Convex Subproblem and Iteration Algorithm
After the appeal discussion, we summarize the discrete convex sub-problem as shown in Problem 2. The control difference is regarded as the control rate constraint because of the FOH approximation of control. A hard trust region on time dilation is enforced in constraints, while the objective function is augmented with trust region terms of state and control, and two additional penalty terms: the nonlinear term of path constraints and distance penalty term of the no-fly zone. The trust term of time dilation can be omitted since the penalty terms of path constraints and the no-fly zone have shown a good iterative performance in our experiment. Thus, the SCP iteration process is given in Algorithm 1 as follows:
Algorithm 1: Improved SCP algorithm with adaptive non-uniform discretization |
|
Problem 2. Discrete convex subproblem in SCP iteration
where
and
can be obtained from (
25) and (
27), and
are the temporal nodes set to cluster around the no-fly zone. Note that the augmented objective with
and
is the general PTR algorithm for non-uniform scheme.
4. Numerical Results
In this section, the effectiveness of the proposed method is verified compared with a different SCP algorithm, as shown in
Table 2. Two cases are considered to focus on distinct instances of constraint violation. Case 1 does not include the no-fly zone constraint and focuses on potential violations of path constraints as the number of discrete points decreases. With an additional no-fly zone constraint, Case 2 focuses on the phenomenon that the propagated trajectory may pass through the no-fly zone.
The reference area and the mass of the vehicle are
m
2 and
kg. The aerodynamic coefficients depend on the attack angle
(in degrees), while the angle-of-attack profile depends on the vehicle’s velocity:
The remaining simulation parameters are shown in
Table 3, which is the same as [
27].
The subproblem is constructed in YALMIP [
34], a MATLAB modeling toolbox, and solved by ECOS [
21], an open-source convex optimization solver. All of the numerical simulations are running on a personal desktop with an Intel Core i9 3.1 GHz processor.
The initial reference is the trajectory obtained by integrating the original dynamics (
1) with the given initial control, as in [
27]. The convergence condition is selected as
and
, where
and
are the solution of the subproblem.
4.1. Iterative Performance
Comparisons of the state and control profiles for Case 1 and Case 2 are shown in
Figure 3,
Figure 4,
Figure 5,
Figure 6,
Figure 7 and
Figure 8. The account of temporal nodes for the proposed method are as follows. SCP1 and SCP2 are each 40, in which case the constraint violation is observed, and SCP3 is 300 in order to maintain the same accuracy as the above methods.
It can be seen that the solutions of the proposed method are similar to those from SCP2 and SCP3. Note that the solutions of SCP1 are quite different due to its poor convergence. As shown in
Figure 5 and
Figure 8, one shortcoming of the proposed method is that the control jitter in the segments where points cluster is obvious because of the dense time grid around those points.
The iterative performance of the proposed method and comparative methods is shown in
Figure 9. The iteration number of the proposed method is less than the comparative methods. Note that the comparative methods would require more iterations to meet the convergence condition when the objective function is near the optimal value, while the proposed method required fewer iterations, demonstrating its fast convergence. What is more, for the non-uniform scheme, SCP1 with light penalty weighting
in
reached the maximum number of iterations, while that with heavy weighting
would show negligible alterations for the distribution of temporal nodes; thus, the addition of time interval dilation is not necessary in this case.
For the non-uniform scheme, the discrete point distribution and the change of time dilation,
, with iterations are shown in
Figure 10. It can be observed that there is no obvious clustering rule for SCP1, while the position of discrete points always changes with the number of iterations. In contrast, the discrete points of the proposed method would cluster after a few iterations. In addition, the results of time interval change with iterations show that the proposed method has stable convergence performance, since the change between two adjacent iterations decreases progressively, which means the result trajectory would become increasingly similar.
In order to guarantee that the result trajectory of the SCP process is feasible to meet the original dynamics, the residual error between the optimized results and the trajectory obtained by integrating the original dynamics is measured. The residual error results for Case 1 and Case 2 are shown in
Table 4 and
Table 5, averaged over 50 simulation runs.
It can be observed that for SCP2 and SCP3, the precise discretization can guarantee commensurate accuracy with fewer temporal nodes, while more than 200 nodes are needed to achieve the same result in [
27], which demonstrates the effectiveness of the precise discretization. Note that, with the non-uniform scheme, SCP1 showed worse convergence performance and low accuracy as mentioned above, which means that the feasibility of the result trajectory is not guaranteed. In contrast, the proposed method overcomes the shortcomings, maintains the advantage of the precise discretization, and shows better convergence performance, as shown in
Table 4 and
Table 5.
4.2. Constraint Satisfaction Performance
As mentioned above, for the uniform precise discretization, the propagated trajectory may violate the path constraints between temporal nodes due to the sparse time grid. This phenomenon was observed to occur when the account of temporal nodes decreased to 40. Thus, we choose the case of 40 nodes for presentation. The path constraints of the proposed method for Case 1 are shown in
Figure 11 and contrasted with SCP1 and SCP2.
It can be observed that the heat load would touch the boundary during the initial state of flight. The constraint violations occur for the propagated trajectory of SCP1 and SPC2 due to the sparse time grid. In contrast, the propagated trajectory of the proposed method satisfies the path constraints over the time horizon. In addition, the discrete points would cluster around the peak value of heat load, as shown in
Figure 11.
Furthermore, the same phenomenon may occur for the no-fly zone constraints as well. Case 2 focuses on the no-fly zone constraint violation, and the numerical results are shown in
Figure 12. In this case, all the path constraints were satisfied for all methods in our experiments. It can be observed that the propagated trajectories of SCP1 and SCP2 both pass through the no-fly zone. Note that, for the non-uniform scheme, the discrete points of SCP1 do not cluster around the no-fly zone. The 3D trajectory of tbe proposed method is shown in
Figure 12b; the points in red are set to cluster around the no-fly zone. The propagated trajectory of the proposed method skimmed over the no-fly zone, while the time interval between adjacent red points is the minimum set in (
11), which indicated that the point distribution of the proposed method is better than SCP2.
5. Conclusions
This paper proposes an improved SCP algorithm for the hypersonic entry problem using a novel adaptive non-uniform discretization. The proposed method has advantages in performance of path constraint satisfaction and convergence. Firstly, the proposed method employs an inverse-free precise discretization to ensure high accuracy and real-time performance. Then, an adaptive non-uniform scheme is developed to distribute discrete points adaptively by adding additional penalty terms in the SCP subproblem, which would guarantee constraint satisfaction. Finally, numerical results show that the proposed method achieves a fast convergence while maintaining high accuracy with few temporal nodes. More importantly, the discrete points of the proposed method would cluster around the segment where the constraints may be violated, and the propagated trajectory satisfies all the path constraints over the time horizon even for a small number of discrete points.
Future work will focus on the following points: (1) Due to the non-uniform scheme, a similar idea can extend to the hypersonic entry problem with the waypoint constraint and other problems; (2) The simulation will be carried out on an embedded platform to verify the effectiveness of the proposed method for a limited-power environment; (3) High-performance solvers are considered to further improve the solving speed, such as the proportional integral projected gradient method [
35], a first-order method for the conic convex problem.