Abstract
This paper presents a trajectory optimization algorithm for super-synchronous-transfer-orbit (SSTO) large launch systems by convex optimization. The payload of SSTO launch systems is typically a geostationary equatorial orbit (GEO) satellite, and the time and position of orbital injection are constrained, which is quite different from the case of general satellites. In this paper, the optimal control problem of SSTO large launch systems is formulated considering the terminal constraints including orbital elements and the time-position equation. To improve the computational performance of the algorithm, the terminal orbital element constraints are expressed in the perifocal coordinate system with second-order equations. And then, several convexification techniques and their modified strategies are applied to transform the original trajectory optimization problem into a series of convex optimization problems, which can be solved iteratively with high accuracy and computational efficiency. Considering the time-position constraint of the payload, the flight time updater design method is proposed to correct the error of time during the flight, which lays solid foundation for the subsequent flight phase, guaranteeing that the GEO satellite settles into the required position. Finally, simulation results indicate the high efficiency and accuracy and strong robustness of the proposed algorithm in different special situations including engine failure and time delay. The algorithm proposed in this paper has great development potential and application prospect in onboard trajectory optimization of SSTO launch missions and similar situations.
1. Introduction
Geostationary equatorial orbit (GEO) is a special circular earth orbit, in which satellites remain relatively stationary with the earth, and the ground stations do not have to track the satellites [1]. The GEO attitude is 35786 kilometers, where the visible portion of the earth’s surface is very large. Because of these meaningful properties, the satellites for communication, global weather, radio, etc. are always placed in GEO [1]. However, in order to avoid interference with other satellites in GEO, and achieve predetermined functions, the satellites should stay above the predetermined point on the earth’s equator. In other words, the time and position of orbital injection for GEO satellites are strictly constrained, which is different from the case of general satellites [2, 3]. Because of the high attitude of the GEO, the scale of the launch system is always very large, and the special situations such as engine failure of large launch systems should be considered.
In general, there are three flight phases from the satellite launching to the satellite settling into the GEO. Firstly, the launch system put the satellite and upper stage into a low earth orbit (LEO), which is a circular orbit locating 200 km~400 km away from the surface of the earth. After that, the launch system will park in the LEO for a few minutes before turning on the engine of the upper stage and sending the satellite into a geostationary transfer orbit (GTO). GTO is a highly eccentric orbit, of which the perigee is near the parking LEO, and the apogee is near the GEO. Finally, the satellites transfer into the target GEO. In this paper, super-synchronous-transfer-orbit (SSTO) [4] is adopted to transfer the satellites from LEO to GEO instead of GTO. SSTO has a higher apogee and eccentricity, and the speed increment required by transferring from SSTO to GEO is greater than that required by transferring from GTO to GEO, which is beneficial for saving the propellant in the process of transferring orbit plane to GEO. The mission of launch systems for GEO satellites is to put the payload into the predetermined SSTO accurately, during which the guidance algorithm plays an important role.
As for the ascent phase of launch systems, the iterative guidance method has been widely and successfully used in the past few decades [5, 6]. However, for most rocket engines, the thrust magnitude is uncontrollable, and the time and position of orbital injection cannot be constrained based on the iterative guidance method. When the deviation of the flight time is large, the iterative guidance method cannot correct the error, which will cause more propellant consumption for orbit transfer, or even lead to the failure of settling into the predetermined point of GEO.
In recent years, with the development of computational technologies, the online trajectory optimization method has been developed rapidly, which provides new ideas for solving the optimal guidance and control problems in aerospace applications [7]. Various optimization theories and algorithms have been developed for online or onboard trajectory optimization of different vehicles [8]. There are three important performance indexes for online trajectory algorithms: accuracy, computational efficiency, and robustness [7]. For guidance algorithms, accuracy is certainly the most important index. Considering the limited onboard computing resources, the computation amount and memory occupation of the algorithm should be reduced. The algorithm should also be robust enough to handle the deviation or other situations when the real trajectory is quite different from the nominal one.
Among the trajectory optimization algorithms developed in recent years, including indirect methods and direct methods [9–14], convex optimization algorithms have great advantages in onboard aerospace applications because the convex optimization problem can be solved in polynomial time with no need for initial guesses supplied by the user [15]. However, most of the original trajectory optimization problems are infinite and nonconvex, which cannot be solved by convex optimization methods directly [16]. Thus, the convexification techniques for trajectory optimization problems are widely studied. In general, lossless convexification and successive convexification are two effective methods to convexify the nonconvex term of the trajectory optimization problem. Lossless convexification was proposed to solve general optimal control problems [17–19]. For aerospace applications, lossless relaxation and convexification are always applied to convexify the thrust magnitude constraints [20]. Lossless convexification has been successfully applied to solve the optimization problem of landing vehicles [21], launch vehicles [20], missiles [22], etc. Successive convexification has a broader range of applications. With successive linearization, all the nonlinear and nonconvex terms can be converted into linear ones based on a known solution, which is very simple and of practical significance [16]. On this basis, the original optimization problem can be transformed into a series of convex optimization subproblems, which can be solved iteratively until the solution converges. The convex optimization has been successfully and widely applied to launch vehicles [23], unmanned aerial vehicles [24], hypersonic glide vehicles [25], etc. [26–30].
To solve the trajectory optimization problem of SSTO launch systems accurately and rapidly, a convex-optimization-based algorithm is applied in this work. In general, the terminal constraints of launch missions are expressed as six orbital elements [1], including the semimajor axis, eccentricity, inclination, longitude of the ascending node, argument perigee, and true anomaly. However, the calculating formulas of these elements are very complex and strongly nonlinear, which negatively affects the computational efficiency and convergence of the algorithm during the iteration. To improve the computational performance, the trajectory optimization problem, the terminal constraints are given in the perifocal coordinate system. The perifocal coordinate system is defined in the target orbital plane based on the geometric feature of the elliptical orbit, and the terminal constraints of orbital elements can be formulated simply as linear or second-order equations. After that, several convexification techniques are applied to transform the original trajectory optimization problem into a series of second-order cone programming (SOCP) problems, which can be solved by the primal-dual interior-point method (IPM). IPM is a typical and widely used algorithm for SOCP. For any given initial guess and accuracy, a globally optimal solution can be found by IPM within the predetermined upper bound of iteration times on condition that the feasible solution exists [10]. To ensure the problem can be solved by IPM, the flip-Radau pseudospectral discretization method is adopted to convert the continuous and infinite problem into a finite one, and the thrust magnitude constraint is relaxed based on lossless convexification. Other nonconvexity of the problem is handled by successive convexification. The main procedure of successive convexification is described as follows: linearize the nonlinear part of the solution obtained by the initial guess or the previous iteration and then solve the linearized (convexified) problem iteratively. If the initial guess is not accurate, the solution cannot be found at the beginning of the iteration. To avoid this, the relaxation method is applied [31], but the computational efficiency declines seriously with the additional relaxation variables. In this paper, the successive convexification is modified to improve both the robustness and computational efficiency of the algorithm. Considering the accuracy of linearization and equality constraints, the parameters of optimization and relaxation variables are changed or removed adaptively with the convergence of the solution during the iteration. In this way, the trajectory optimization problem can be solved based on convex optimization with satisfactory computational performance.
Another research focus of this paper is the terminal time and position constraints. As for the GEO satellites, the position of the subastral point is stationary and strictly constrained. In general, the low-thrust stage is often applied to transfer the satellite from SSTO to GTO [32], and the ability of the stage to correct the deviation is weak. Therefore, the launch mission of SSTO must ensure that the satellite arrives at the point of orbit transfer at a predetermined time, which lays the foundation for orbit transfer and the subsequent on-orbit missions. In other words, the trajectory optimization algorithm should concentrate on the orbital elements and the time/position of orbital injection. In this paper, the terminal constraint of the orbital injection point is formulated as a function of time and eccentric anomaly equivalently. However, as the thrust magnitude of a rocket engine is uncontrollable, it is difficult to constrain the time and position of orbital injection during the flight from LEO to SSTO. To this end, the flight time of parking in the LEO updater design algorithm is studied in this paper. With the adjustment of the time parking in LEO, the time-position of orbital injection can meet the requirement. The algorithm proposed in this paper has strong robustness and can solve the online trajectory optimization problem even under partial engine failure or launch time delay.
This paper is organized as follows. In Section 2, the trajectory optimization problem of SSTO launch systems is formulated. In Section 3, the original trajectory optimization problem is transformed into SOCP subproblems by several modified convexification techniques. In Section 4, the flight time of parking in the LEO updater design algorithm is studied. In Section 5, simulation is carried out to compare the proposed optimization method with the traditional optimization method under different conditions. In Section 6, some conclusions are given.
2. Problem Formulation
Firstly, we formulate the dimensionless equations of motion of SSTO launch systems in the Earth Center Inertial Coordinate System [1] as follows: where and are the dimensionless inertial position and velocity vectors, respectively; is the mass of the system. is the thrust vector and also represents the attitude angle of the system. is the gravitational acceleration magnitude on the surface of the Earth. is the specific impulse of the engine. The distance is normalized by the radius of the Earth at the equator , the time by , and the velocity by [33].
When the payload is boosted into an LEO, the rocket engine turns off and the thrust magnitude . After a few minutes of unpowered flight, the rocket engine of the upper stage turns on, and the payload is boosted into an SSTO. For SSTO launch systems, the thrust magnitude of a rocket engine is uncontrollable, and the magnitude of the thrust vector is constrained:
Considering the strict constraints of the terminal position and velocity in the noninertial coordinate system for GEO satellites, the terminal constraints for SSTO launch systems are also important. Traditionally, the terminal constraints of launch systems are expressed as orbital elements: the semimajor axis, eccentricity, inclination, longitude of the ascending node, argument perigee, and true anomaly [1]. As for GEO satellite launch missions, the position of the subastral point is strictly constrained; so, the time and position of the SSTO injection also need to be considered. In other words, the payload should settle into the nominal SSTO, and the position as a function of time should be the same as the nominal one.
For convenience, the terminal constraints of SSTO launch systems can be expressed in the perifocal coordinate system [34]. is the center of the earth, the -axis points towards the perigee, and the -axis is parallel to the normal of the orbital plane (along the positive direction of the normal). The -axis completes the right-hand coordinate system. And then the terminal constraints of SSTO launch systems can be expressed in the perifocal coordinate system as follows:
Considering the accuracy of the semimajor axis, eccentricity, inclination, longitude of the ascending node, and argument perigee, the following five terminal equality constraints must be satisfied: where is the semiminor axis and ; is the distance from the center of earth to the center of ellipse orbit and . is the required magnitude of angular momentum. The subscript “” represents the final value of the parameters.
In addition, the system must settle into SSTO at a certain time and position. As an equivalent transformation, the time and position of injection should satisfy the following condition:
As shown in Figure 1, and are the expected time and position of injection, and and are the actual time and position of injection. is the unpowered flight time from to . is the eccentric anomaly. If satisfies, then The payload can be considered to settle into SSTO at the expected time and position equivalently. According to the Kepler’s equation, can be calculated by where is the flight time starting at the apogee, which can be calculated as a function of eccentric anomaly [1]:
And eccentric anomaly can be calculated simply in the perifocal coordinate system as follows:
Considering Eqs. (4) and (7), the error of time can be defined as
In conclusion, the terminal constraints of SSTO launch systems are Eqs. (3)–(8). As for most launch systems, because the thrust magnitude is uncontrollable, only five terminal constraints can be guaranteed; so, the time and position of the injection are uncontrollable. In Section 4, the strategy to decrease the error of time is proposed in detail.
Considering the following mission of orbit transfer, the trajectory optimization problem of SSTO launch systems is defined as an optimal control problem to achieve the minimum fuel consumption: subject to where the variable , , is the dry mass of the launch system. Equation Eq. (11) is inertial constraints. It is obvious that the original trajectory optimization problem (Eqs. (9) and (11)) is nonconvex. In the next section, the trajectory optimization problem is transformed into a series of SOCP subproblems, which can be solved directly by the primal-dual interior point method rapidly with good accuracy.
3. Convexification
In this section, the original trajectory optimization in Section 2 is transformed into a series of discrete SOCP subproblems by pseudospectral discretization, lossless convexification, and successive convexification. For better computational efficiency, a second-order correction algorithm and an improved relaxation method for successive convexification are proposed.
3.1. Pseudospectral Discretization and Lossless Convexification
To meet the requirement of the convex optimization method, the continuous infinite dynamical constraints are always converted into a finite set of equality constraints by discretization. Considering the accuracy of the terminal constraints in different phases, the flip-Radau pseudospectral discretization method, by which the collocation points discreted within the domain , is adopted in this paper: where is the flip-Radau pseudospectral differentiation matrix [35], is the right side of the differential dynamic equations Eq. (1), and is the collocation points within the domain . is the number of collocation points. is state variables including , , and .
On the other hand, the constraint of the thrust magnitude in Eq. (2) is a nonconvex equality constraint. For launch systems, the most commonly used method to handle the nonconvexity of the thrust magnitude constraint is feasible domain relaxation. The feasible domain can be relaxed from a spherical-shell region to a solid sphere, and the relaxed constraint of the thrust magnitude is
The relaxed optimization problem has the same optimal solution as the original one, which means the transformation is equivalent, and the convexification technique is called lossless convexification [17–19]. The equivalence of the transformation can be proved based on optimal control theory. The detailed proof can be found in [20].
After discretization and lossless convexification, the original trajectory optimization problem formulated in Section 2 can be expressed as subject to
The equality constraints of terminal conditions and discretized dynamical equations are still nonlinear (nonaffine), which is not suitable for the SCOP-based method. This problem can be handled by successive convexification in the next subsection.
3.2. Improved Successive Convexification
By linearization and the successive solution procedure, successive convexification has been successfully applied to converting the nonlinear equality constraints into affine ones. It is a popular and simple technique to handle the residual nonconvexity of the optimization problem in Section 4.1. In this paper, the nonlinear equality constraints Eq. (12) and the terminal constraints Eq. (3) are linearized repeatedly at the previous iteration with a known solution. For convenience, all nonlinear equality constraints are expressed as , and they can be linearized by first-order Taylor series expansion [16]: where .
Considering the accuracy of linearization, the update variable in Eq. (21) is added as a penalty term in the performance index function. And the convex optimization problem is formulated as COP1: subject to where is the penalty coefficient of .
In this way, the trajectory optimization problem of the SSTO launch system is converted into a series of SOCP subproblems, which can be solved by IPM iteratively. COP1 has optimal variables, including at every discrete point. It should be noted that the error of time is ignored in this section because it is hard to control the position and time of the injection when the thrust magnitude is uncontrollable. The terminal constraint is considered in Section 4.
However, because of the error of the initial guess, the feasible solution considering the strict constraint Eq. (17) may not exist at the beginning of the iteration, even if the solution to the original problem can be found. This phenomenon is called “artificial infeasibility,” and the detailed analysis can be found in [31]. This problem is solved by relaxation and the penalty strategy. Firstly, the constraint Eq. (17) is replaced by where is the relaxation variable for equality constraints. The dimension of is the same as the number of equality constraints.
And the relaxation variable is added as a penalty term in the performance index function. The relaxed convex optimization problem is called COP2: subject to where is the penalty coefficient of .
In the performance index function of COP2, there are two intercoupling penalty coefficients and . The accuracy of the equality constraints increases as increases and decreases, and at the same time, the accuracy of linearization decreases. However, if the error of linearization is too large, the high accuracy of the equality constraints becomes insignificance. The error of linearization can be defined as where is the number of equality constraints, . If is larger than , the accuracy of linearization is unacceptable, the penalty coefficient needs to be increased ,and the solution of this iteration should be given up. If is smaller than , can be increased to improve the accuracy of equality constraints and the rate of convergence.
By comparing COP2 with COP1, it can be found that considering the relaxation variable , the number of optimal variables increases to . It is well known that as the number of optimal variables grows, the computational time of the optimization problem grows exponentially. That is, the computational efficiency of COP2 is much worse than that of COP1. Actually, when the accuracy of equality constraints reaches a certain degree , artificial infeasibility will not happen, and the relaxation variable and penalty terms in the performance index can be ignored.
In conclusion, considering both accuracy and computational efficiency of the algorithm, the trajectory optimization problem can be solved iteratively by the following steps (as shown in Figure 2): (1)Initialization: input the initial variables including ,,, , , and , the permitted range of penalty coefficients and , and the accuracy requirement (2)Set (3)Solve the problem COP2 by IPM, obtain the optimal solution , and calculate by Eq. (25)(4)If , set , and give up the solution . Otherwise, update the optimal variable and set (5)If , set and . If , ignore the relaxation variable and the corresponding penalty term(6)If , solve the problem COP1 iteratively. Otherwise, return to step 3(7)Check the convergence condition:
If Eq. (26) is satisfied, the trajectory optimization problem is solved, and the optimal solution is .
In this way, the trajectory optimization problem is solved without considering the error of time in Eq. (8). In the next section, to ensure the injection accuracy of GEO satellites, the flight time in LEO for launch systems is calculated to meet the accuracy requirement of .
4. Flight Time Updater Design
As for the launch mission of GEO satellites, the systems always settle into a circular LEO, and then, SSTO is used to transfer satellites from LEO to GEO. Between these two flight phases, the launch systems have a few minutes of unpowered flight, and then the rocket engine turns on to settle the payload into SSTO. In Section 3, the trajectory optimization algorithm is proposed, but the error of time is ignored. In this section, the strategy of LEO flight time updater design is proposed, and the terminal constraint in Eq. (8) is considered.
By the analysis of the trajectory optimization problem of SSTO launch systems in Section 3, when an initial condition is given, the error of time can be calculated by Eq. (8) based on the optimal solution. In other words, the error of time can be treated as a function of the initial condition . As for the unpowered flight phase in LEO, the initial condition of the powered phase from LEO to SSTO can be easily calculated when the orbital elements of LEO are given. When the LEO is circular and the radius of the orbit is , the initial true anomaly is ,the flight time in LEO is , and the true anomaly at is where is the gravitational parameter of the Earth.
And then the initial condition can be easily calculated in the perifocal coordinate system as follows:
It should be noted that all variables in Eqs. (27) and (28) are dimensional to facilitate understanding, which is different from that in Section 3. When they are given as the initial condition of the trajectory optimization problem, all variables must be nondimensional.
Based on Eqs. (27) and (28), and the trajectory optimization algorithm in Section 3, the error of time can be calculated when the flight time in LEO is given. And can be calculated iteratively by the Newton method to solve the equation:
However, is calculated based on a complex optimization procedure, and the partial derivative is hard to calculate accurately by the numerical method when is too small. So, a modified Newton method is given, and the detailed calculation procedure is as follows (as shown in Figure 3): (1)According to the nominal trajectory, initialize the flight time in LEO . Set (2)Calculate the initial condition by Eqs.(27) and (28) and solve the trajectory optimization problem by the algorithm proposed in Section 3, and then the error of time can be calculated by Eq.(8) based on the optimal solution(3)Set , (4)Calculate in the same way as step 2(5)If ( is the accuracy requirement of ), the flight time in LEO is ; else, go to step 6(6)Update by(7)Set , and go to step 4
5. Simulation and Analysis
In this section, simulation experiments are carried out by taking the whole flight from LEO to SSTO. To verify the robustness of the algorithm proposed in this paper, two kinds of special conditions are considered: (1)One of the four engines breaks down and cannot work, which means the thrust magnitude and the rate of mass flow both decrease by 25%(2)Flight time delay caused by a fault in the previous flight phase, such as the launch delay, system fault, and flight deviation
In this paper, the convex optimization problems are solved by the MOSEK software [36], and the simulation results are compared with the solutions obtained by a traditional optimization method, which employs the hp-adaptive Radau pseudospectral method and the general NLP methods [31]. All numerical simulations in this paper are performed on a laptop with Intel Core i7 CPU 2.80GHz.
Parameters of the SSTO launch system are listed in Table 1. Parameters of the initial condition and target orbit (SSTO) are listed in Table 2. The variables contained in this section’s figures and tables are with respect to the Earth Centered Inertial Coordinate System.
5.1. Engine Failure
In this subsection, we assume that the thrust magnitude and the rate of mass flow both decrease by 25% at the same time. The number of the collocation points is , the penalty parameter is , and the permitted range of penalty coefficients is . , , , and the accuracy requirements .
As shown in Figures 4–9 and Table 3, the optimal solutions to the trajectory optimization problem obtained by both methods are very similar, and they are quite different from the nominal trajectory because of the failure of the engine. In other words, the algorithm proposed in this paper can solve the trajectory optimization problem without good initial guesses in case of engine failures. Moreover, the accuracy and optimality of the optimal solution are proved by comparison with the traditional optimization method. In Figure 9, the constraint of the thrust magnitude is active during the whole flight, which demonstrates the validity of the lossless convexification in Section 3.1. As shown in Table 3, the average CPU time is 2.23 s for the algorithm presented in Section 3.2, which is only 5.8% of the average CPU time for the traditional method. It takes 7 iterations and 15.61 s’ CPU time to calculate the flight time in LEO by the algorithm in Section 4. Due to the long time parking in the LEO, there is enough time to solve the trajectory optimization problem.
5.2. Time Delay
In this subsection, we assume that the flight time delays by 50 s because of the fault in the previous flight phase. All the parameters of the optimization algorithm are the same as those in Section 5.1. The simulation results are as follows:
Compared with the work in Section 5.1, similar simulation results (shown in Figures 10–15, Table 4) and research conclusions can be obtained. Another significant research result is that with the same optimization parameters, the trajectory optimization problem can be solved in different situations, which further proves the robustness of the algorithm. Similar experiment results can also be obtained based on different launch systems and missions. As for the time delay, the average CPU time is 2.48 s for the algorithm presented in Section 3.2, which is 5.7% of the average CPU time for the traditional method. It takes 8 iterations and 19.84 s’ CPU time to calculate the flight time in LEO by the algorithm in Section 4. The deviations of terminal orbital elements shown in Tables 5 and 6 verify the accuracy of the equation of terminal constraints Eq. (3) and the flight time updater design algorithm in Section 4.
The above simulation results demonstrate that the algorithm proposed in this paper has good robustness, accuracy, and computational efficiency, which indicates great development potential and application prospect in onboard trajectory optimization.
6. Conclusion
This paper presents a convex optimization algorithm for SSTO launch systems considering the orbital elements and time-position constraints. For convenience and better computational performance, the optimal control problem including the terminal constraints is given in the perifocal coordinate system. And then the flip-Radau pseudospectral method is adopted to convert the trajectory optimization problem into a finite problem. Lossless convexification is also utilized to convexify the constraint of the thrust magnitude. To improve the robustness and computational efficiency of the algorithm, successive convexification and its modified method are proposed. In this way, the trajectory optimization problem is transformed into a convex one, which can be solved by IPM accurately and rapidly. To correct the flight time deviation of SSTO launch systems, the LEO flight time updater design algorithm is proposed. Finally, the algorithm proposed in this paper is tested under two special conditions: engine failure and time delay. Compared with the traditional optimization method, the proposed algorithm demonstrates stronger robustness, higher accuracy, and higher computational efficiency. The improved convex approach and flight time design method proposed in this paper have great application potential in onboard trajectory optimization of SSTO launch systems and other similar systems, especially under nonnominal conditions. In our follow-up work, we will modify the proposed algorithm for the flight from SSTO to GEO and concentrate on the improvement of injection accuracy and robustness of the algorithm.
Data Availability
The numerical simulation data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grant 51809138 and China Postdoctoral Science Foundation under Grant 2019M651837.