Fuel-Optimal Ascent Trajectory Problem for Launch Vehicle via Theory of Functional Connections

In this study, we develop a method based on the Theory of Functional Connections (TFC) to solve the fuel-optimal problem in the ascending phase of the launch vehicle. The problem is first transformed into a nonlinear two-point boundary value problem (TPBVP) using the indirect method. Then, using the function interpolation technique called the TFC, the problem’s constraints are analytically embedded into a functional, and the TPBVP is transformed into an unconstrained optimization problem that includes orthogonal polynomials with unknown coefficients. This process effectively reduces the search space of the solution because the original constrained problem transformed into an unconstrained problem, and thus, the unknown coefficients of the unconstrained expression can be solved using simple numerical methods. Finally, the proposed algorithm is validated by comparing to a general nonlinear optimal control software GPOPS-II and the traditional indirect numerical method. The results demonstrated that the proposed algorithm is robust to poor initial values, and solutions can be solved in less than 300ms within the MATLAB implementation. Consequently, the proposed method has the potential to generate optimal trajectories on-board in real time.


Introduction
With the recent development in space exploration, launch vehicles are very important as they are the only means for humans to explore space from the earth. In general, a launch vehicle mission has been planned over a long period, and the trajectory was designed in advance, and it cannot be updated during flight, which means it is not robust or flexible. The fast launch and trajectory reconstruction are the main research of the guidance system, and both need rapid trajectory planning technology. Rapid trajectory planning can shorten the launch mission cycle and quickly update the trajectory in case of thrust failure during the flight of the vehicle to ensure the success of the mission.
The primary aim of the trajectory planning algorithm is to solve the optimal control problem that is generally based on nonlinear dynamics, which achieves specific performance indicators under the constraints of state and control variables. The solution of such problems is mainly achieved using the indirect method [1][2][3] and the direct method [4][5][6]. The direct method transforms the optimal control problem of continuous space into a nonlinear programming problem and uses a numerical method to directly optimize the performance index [7][8][9]. Although the direct method, represented by a sequential quadratic programming algorithm with the pseudospectral discrete, has advanced a lot over a period of time, it still has considerable issues between the on-board application. The algorithm for the general nonlinear programming, for example, the famous sequential quadratic programming algorithm, had low algorithm efficiency and low sensitivity to the initial value and was unable to guarantee convergence in the past, but now, these issues have easily resolved by using a convex optimization method. Ralph Rockafellar, a renowned mathematician, pointed out that the key determinant of the performance of a numerical optimization algorithm is neither the linearity nor nonlinearity of the problem, but the convexity or nonconvexity of the problem [10]. In 2007, JPL proposed lossless convexity technology for the dynamic descent guidance of the Mars lander [11]. After that, a systematic summary of the research and development of lossless convexity technology is presented in [12]. The applicability of the lossless convexity method is extended to general linear systems with multiple state constraints using the concepts of control theory. Unfortunately, only a few nonconvex constraints can be used for lossless convexification, and there is no analytical convexification method for the system dynamic constraints in the trajectory optimization problem of aircraft. A method based on the Newton-Kantorovich/pseudospectral and sequence convexification is used for the ascending phase of the launch vehicle [13]. However, the sequential convex optimization algorithm is a convexification method based on linearization, which increases the dependence on the reference trajectory. This not only offers higher requirements for the reference trajectory but also annihilates the advantage of the convex optimization algorithm that does not rely on the good initial value. Nevertheless, considering the rapidity of the convex optimization algorithm in solving convex problems, in recent years, trajectory planning based on the convex optimization algorithm, such as planetary landing [14][15][16], rocket ascent guidance [17], and entry guidance [18], has been widely studied.
The indirect method solves the optimal control problem by using the classical variational method and the Pontryagin Minimum Principle, derives the first-order necessary conditions of the optimal control, and transforms the optimal control problem into a two-point boundary value problem (TPBVP) [19] that is comprised of initial conditions, Hamiltonian differential equations, optimal conditions, and terminal boundary conditions (including terminal transversal conditions and terminal constraints). However, since the convergence radius of the indirect method is small, and the convergence of the numerical iteration is extremely sensitive to the initial value estimation, which requires a higher accuracy of the initial value estimation, determination of the initial value is highly difficult. To overcome this problem, the deep learning algorithm is used to obtain a higher accuracy of initial value estimation and a better target shooting success rate [20][21][22]. In general, the higher sensitivity of TPBVP to the initial value makes the problem difficult to solve. Therefore, although the indirect method yields a more accurate solution, it is rarely used in practice.
Recently, a mathematical framework called the Theory of Functional Connections (TFC) has been proposed in [23] to derive the expressions with embedded constraints. The expressions, called constrained expressions, are composed of functionals and functions of functions. The constrained expression is written as where gðtÞ is a free function and Φ k ðtÞ are the switch functions composed of the support function s k ðtÞ with unknown coefficients α k . The support function is a set of linearly independent functions. If one of the switch functions is equal to 1, the constraint it is referencing is evaluated; otherwise, that is, if it is equal to 0, all other constraints are evaluated. The switch functions can be expressed asΦ k ðtÞ = s k ðtÞα k . The projection functionals p k ðtÞ are derived by constraint functions. The constrained optimization problem is transformed into an unconstrained one using the constrained expressions, which reduces the search space of the solution to the admissible solutions that satisfy all constraints. Finally, using common basis functions such as Chebyshev polynomials or Legendre polynomials to express the free function and then using the least-square method to find its unknown coefficient, the solution of the problem can be found. In [24,25], the TFC algorithm quickly solves the nonlinear differential equations and obtains high-precision solutions. In [26,27], it is applied to fixed-time asteroid landing and optimal energy landing.
The results show that the solution time is basically less than 100 ms, which proves that the algorithm has real-time application potential. This paper is organized as follows. Section 2 gives a brief description of the TFC mathematical framework to solving TPBVP. In Section 3, the fuel-optimal ascent trajectory problem is described in detail, and the necessary conditions are derived. In Section 4, the fuel-optimal problem in the ascending phase of the launch vehicle is formulated using the TFC framework. Finally, the result and discussion are provided in Section 5.

Theory of Functional Connections
In this section, we present an outline of the TFC mathematical framework and a method for solving secondorder TPBVP with TFC.
2.1. TFC for TPBVP. In general, trajectory optimization problems are second-order TPBVP [28], which is expressed as F t, y t ð Þ, _ y t ð Þ, € y t ð Þ ð Þ = 0 subject to : where t 0 and t f represent the initial time and the terminal time, respectively, and y 0 , y f , _ y 0 , _ y f are the initial and terminal constraints, respectively. As mentioned earlier, the constraints are expressed by (1), and then, (2) is simply rewritten as where p k ðtÞ is expressed as International Journal of Aerospace Engineering The algorithm to derive the term of the constrained expression is given as follows: (1) Choose s k ðtÞ, which are k linearly independent support functions (2) Calculate switching functions Φ k ðtÞ as a linear combination of the support functions with k unknown coefficients (3) Formulate a system of equations to solve for the unknown coefficients based on Φ k ðtÞ The support functions are defined as s k ðkÞ = t k−1 , according to [28]. The switching function is then obtained by solving equations. When the first switch function is activated, the equations are The equations of the first switch functions are combined into the matrix form as The boundary conditions of (2) are effectively embedded within the constrained expression by substituting (4) and (8) into (3). Then, by substituting (3) into the vector differential equation Fðt, yðtÞ, _ yðtÞ, € yðtÞÞ, the constrained TPBVPs are transformed into an unconstrained problem. According to (3), yðtÞ is replaced by gðtÞ; thus, the original vector differential equation Fðt, yðtÞ, _ yðtÞ, € yðtÞÞ is transformed intoFðt, gðtÞ, _ gðtÞ, € gðtÞÞ, which is only a function of t, the free functions gðtÞ, and their derivativeŝ As mentioned earlier,Fðt, gðtÞ, _ gðtÞ, € gðtÞÞ is unconstrained because the boundary is represented by the switch function Φ k ðtÞ and the projection function p k ðtÞ.
After determining the switch function and the projection function, we next discuss the free function gðtÞ.

Definition of the Free Function.
In selecting a free function, we are essentially looking for the best function approximator. A natural choice for the free function is a linear combination of basis functions, as this is capable of spanning the entire function space that the basis spans, as the number of basis functions approaches infinity. The free function is expressed as where ξ are m × 1 unknown coefficients and h are m basis functions.
Next, the problem domain t is mapped to the domain of the basis functions z, and Chebyshev and Legendre polynomials are commonly used, the domains of which are defined in ½−1, 1. To implement the basis functions, a map between t and z is defined as By using (11), the derivatives of gðtÞ are computed as

Domain Discretization.
For solving TPBVPs numerically, the domain t ∈ ½t 0 , t f must be discretized by N + 1 points. The common method is uniform distribution, but the advantage of Chebyshev-Gauss-Lobatto collocation points is that when the number of basis functions increases, the condition number should also increase slowly, which is useful for improving computational efficiency. The Chebyshev-Gauss-Lobatto collocation points are defined as 3 International Journal of Aerospace Engineering Thus, the new vector differential equationFðt, gðtÞ, _ g ðtÞ, € gðtÞÞ becomesFðz, ξÞ, where the unknown coefficient ξ is the variable that needs to be solved.Fðz, ξÞ is expressed in the form of loss functions at each discrete point By setting L i = 0, the unknown coefficient ξ is solved using optimization schemes such as iterative least-squares.
To solve the nonlinear least-square problem, we need the Jacobian matrix of the loss function, which is written as : ð15Þ The estimation is updated by where Δξ = ðJðξÞ T JðξÞÞ −1 JðξÞ T LðξÞ. The iterative process stops when the convergence tolerance is met: where ε is the stopping criterion that is defined by the user. Figure 1 shows the outline of the TFC framework.
Recently, the position of numerical calculation in the current guidance and control field is emphasized in [29], which, based on numerical calculation, is called as the Computational Guidance and Control (CG&C). The CG&C replaces offline planning and closed-loop guidance with on-board computing, which is more robust, more accurate, and more flexible and can adapt to more complex environments and missions, but offers high requirements for computational efficiency. As mentioned earlier, the algorithms used in the CG&C are basically divided into two: direct and indirect. In this study, we used the indirect method, because the optimal control problem is transformed into TPBVP, and then, the TFC method is used to transform and solve the problem.
We found that the TFC method is quite similar to the collocation method in which all three methods using orthogonal polynomials over the global domain, according to the different method. In the collocation method, the states and costates are expanded by using orthogonal polynomials, and the boundary conditions are considered as part of the optimization scheme [30]. It is similar that the pseudospectral method used orthogonal polynomials like Chebyshev or Legendre polynomials to approximate the states and costates, and the boundary conditions are also considered as part of the optimization scheme. The TFC method may use a similar operation, but the fundamental difference between the TFC method and the other two methods lies in handling the constraints of the problem. The TFC method uses orthogonal polynomials to expand the free function gðtÞ in a constrained expression and then expresses the problem constraints analytically step by step as mentioned above, which can reduce the search space of the solution and thus improve the computational efficiency. In fact, the advantages of the TFC method are presented in [27]; the results in [27] show that the TFC method is two orders of magnitude faster than the pseudospectral method in a fixed time optimal control problem and one order of magnitude faster than the pseudospectral method in a free time optimal control problem.

Fuel-Optimal Problem in the Ascending
Phase of the Launch Vehicle In this section, the problem is transformed into a TPBVP, and the first-order necessary conditions and transversal conditions of the problem are derived using the Pontryagin Minimum Principle.

Dynamical Model.
In this section, the last stage of the launching vehicle is studied. The dimensionless equations of motion of a three-dimensional (3-D) launch vehicle can be expressed in the Earth Center Inertial Coordinate System as follows: where r is the inertial position, which is normalized by the radius of the Earth R 0 = 6378145 m. v is the velocity, which is normalized by ffiffiffiffiffiffiffiffiffi ffi R 0 g 0 p , in which g 0 = 9:81 m/s 2 represents the gravitational acceleration magnitude on the surface of the Earth. The mass of the launch vehicle is denoted by m. The thrust is denoted by T = TI b , where I b is the unit vector of the body axis satisfying For most launch vehicles, the mass flow is uncontrollable; thus, the thrust magnitude T = kTk is constant and uncontrollable during the same flight phase. The gravity acceleration 4 International Journal of Aerospace Engineering a g = −r/r 3 , where r = krk. The specific impulse of the engine is denoted by I sp . The differentiation of the equations in (18) is with respect to dimensionless time normalized by ffiffiffiffiffiffiffiffiffiffiffi R 0 /g 0 p . As the mass flow is constant, the optimal ascent problem is treated as a minimum-time problem: subject to where I b = −λ v /λ v is called the primer's vector, according to Lawden's theory [31]. Thus, (22) is rewritten as The first-order necessary conditions for optimality then give the differential equations of the costate variables: The transversality condition is expressed as

Solution via TFC
4.1. TPBVP in TFC Framework. As mentioned in the previous section, to find the optimal state, the following nonlinear TPBVPs must be solved: Projection function p k (t) Support function s k (t) Switch function k (t) Free function

International Journal of Aerospace Engineering
where (26) and (27) are subject to Additionally, there is a redundant equation and can be removed by the TFC constraints. The derivative of rðtÞ is exactly the function vðtÞ because TFC constraints are analytical expressions; thus, (26) can be disregarded. The problem is now reduced. The new equations are expressed as To solve the above equations, the TFC constraints with rðtÞ, λ r ðtÞ, λ v ðtÞ need to be constructed by the TFC method. The unknown coefficients in the TFC constraint expressions are expressed as ξ a , ξ λ r , ξ λ v , and vðtÞ, aðtÞ, _ λ r ðtÞ, _ λ v ðtÞ can be obtained by taking the derivative of the TFC constraint expression rðtÞ, λ r ðtÞ, λ v ðtÞ, respectively.
The initial and terminal constraints of the problem discussed in this paper are position and velocity constraints, respectively. The TFC constraint expressions of rðtÞ, vðtÞ, a ðtÞ are written as where the projection function is written as Next, consider constructing a free function. According to (10)- (12), the time domain is mapped to the Chebyshev domain. However, it should be noted that in the time-free TPBVP, the parameter c is a function of t f . Combined with the TFC method, the new unknown variable ξ t is used to represent the parameter c, and the optimal time is obtained by solv-ing ξ t . In addition, to ensure the solved final time is positive, b is used instead of b, where b 2 = c. Then, rðtÞ is rewritten as where h 0 = hðz 0 Þ and h f = hðz f Þ.
The expression of the switch function is similar to (8): Next, the TFC constraint expressions of the costates are constructed by the same steps as above: Substitution of the above TFC constraint expression into the loss function yields the loss function with respect to unknown coefficient ξ. Then, the solution of the problem is obtained using the nonlinear least-square method. The unknown coefficient ξ is expressed as The loss function can be expressed as International Journal of Aerospace Engineering To use the nonlinear least-square method, the partial derivative of the loss function needs to be calculated as The terms of (41) are defined by , , These are uniformly discrete according to the number of polynomials m.

Simulations
In this section, we apply the proposed algorithm to the ascent problem of the launch vehicle to verify the feasibility of the algorithm, and the results are compared with those of the GPOPS-II and classical indirect method solutions of the other two algorithms. All numerical results are obtained on a desktop with Intel Xeon E3-1230 3.4 GHz. Table 1 lists the parameters of the launch vehicle and mission in the numerical simulations. Table 2 shows the initial and terminal parameters of the experiment, where the orbital elements corresponding to the terminal position and velocity are also given because the launch vehicle generally uses the orbital elements for the target.
To prove the validity and effectiveness of the algorithm proposed in this study, the results of the algorithm proposed in this paper are compared with those obtained by the traditional indirect method and GPOPS-II. The final time calculated by GPOPS-II is 300.97 s, that obtained by the single shooting method is 301.01 s, and the final time obtained by TFC is 301.25 s. The locations of the vehicle and the velocity vector are provided in Figures 2 and 3, respectively. Figure 4 shows the thrust vector of the launch vehicle. In the bottom part of Figure 4, the purple line representing the sum of squares of thrust directions is equal to 1, which also indicates the validity of the TFC solution, and the other three lines are the vector of the body axis. Considering that the pitch angle and yaw angle are generally used as the guidance command of the launch 7 International Journal of Aerospace Engineering vehicle, the results of pitch angle and yaw angle are also given here. Figures 5 and 6 show the results of the height and the velocity, respectively. It can be seen that the results solved by the three methods are basically the same, which again shows the validity of the TFC algorithm proposed in this study. To quantify the accuracy of the TFC method,   In this simulation, the cost time of the proposed method is 0.23 s and those of GPOPS-II and the single shooting method are 3.88 and 0.34 s, respectively. It is also known that  9 International Journal of Aerospace Engineering the MATLAB programming language is 10 times slower than C++; thus, the algorithm proposed in this study has online application potential. In terms of solution accuracy, the TFC method is not dominant among the three methods, but combined with the analysis of calculation efficiency, it shows that the proposed method is an effective method.
The above results show the comparison between the results of the proposed algorithm, GPOPS-II, and single shooting method, which verifies the validity of the proposed algorithm. Next, the effect of the number of discrete points and polynomials on the algorithm is studied. Table 3 shows that the excessive number of discrete points will not only reduce the calculation efficiency, but also reduce the accuracy of the solution. In addition, the selection of the number of polynomials is also analyzed in this paper. In [28], the number of state polynomials and costate polynomials is the same and not studied separately. In our simulation, the results show that the selection of the number of state polynomials and costate polynomials can be different, and a better result can be obtained. If the number of the state polynomials is selected too much, then not only the calculation performance will be degraded but also the accuracy of the solution will be greatly affected. It is seen from Table 3 that when the number of costate polynomials is as large as the number of state polynomials, the allowed iteration number is reached and the iteration progress will not converge; it means the solution of the problem cannot be solved. In addition, appropriately increasing the number of costate polynomials can increase the accuracy of the solution, but it will reduce the computational efficiency. Thus, the number of discrete points and polynomials should be select carefully.

Conclusion
In this study, we proposed a new approach to solve the fueloptimal problem in the ascending phase of the launching vehicle using TFC. The main conclusions can be summarized as follows: (1) The first-order necessary condition of the optimization problem is constructed by the indirect method; the problem's constraints are embedded in the expression by using the TFC method (2) The constrained optimization problem is transformed into an unconstrained optimization problem by using the TFC method, which reduces the search space of the solution, and a simple root-finding algo-rithm can be used to obtain the solution of the problem (3) The residual of the solution is about 10 −14 or less; for obtaining more accurate numerical solutions, the number of discrete points and polynomials should be selected carefully (4) The proposed algorithm has the potential for online application. The calculation time of the algorithm is within 300 ms with MATLAB programming