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

Authors: Shiyao Li; Yushen Yan; Kun Zhang; Xinguo Li
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

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 300 ms within the MATLAB implementation. Consequently, the proposed method has the potential to generate optimal trajectories on-board in real time.


---

## Body

## 1. 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–3] and the direct method [4–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–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–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–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
(1)yt=gt+∑nk=1Φktpkt,where gt is a free function and Φkt are the switch functions composed of the support function skt 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Φkt=sktαk. The projection functionals pkt 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. Section2 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.


## 2. Theory of Functional Connections

In this section, we present an outline of the TFC mathematical framework and a method for solving second-order TPBVP with TFC.


### 2.1. TFC for TPBVP

In general, trajectory optimization problems are second-order TPBVP [28], which is expressed as
(2)Ft,yt,ẏt,y¨t=0subjectto:yt0=y0,ytf=yf,ẏt0=ẏ0,ẏtf=ẏf,where t0 and tf represent the initial time and the terminal time, respectively, and y0,yf,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
(3)yt=gt+Φ1tp1+Φ2tp2+Φ3tp3+Φ4tp4,where pkt is expressed as
(4)p1t,git=y0i−git0,p2t,git=yfi−gitf,p3t,git=ẏ0i−ġit0,p4t,git=ẏfi−ġitf.The algorithm to derive the term of the constrained expression is given as follows:(1)
Chooseskt, which are k linearly independent support functions(2)
Calculate switching functionsΦkt 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ΦktThe support functions are defined asskk=tk−1, according to [28]. The switching function is then obtained by solving equations. When the first switch function is activated, the equations are
(5)Φ1t0=1,Φ1tf=0,Φ1′t0=0,Φ1′tf=0.The equations of the first switch functions are combined into the matrix form as(6)1t0t02t031tftf2tf3012t03t02012tf3tf2α11α21α31α41=1000.Similarly, application of this method to the other three switch functions yields the following matrix:(7)1t0t02t031tftf2tf3012t03t02012tf3tf2α11α21α31α41α21α22α32α42α31α23α33α43α41α24α34α44=1000010000100001.The coefficientαk is obtained by matrix inversion, and the expression of the switch functions is obtained as
(8)Φ1t=1tf−t03−tf23t0−tf+6t0tft−3t0+tft2+2t3,Φ2t=1tf−t03−t02t0−3tf−6t0tft+3t0+tft2−2t3,Φ3t=1tf−t02−t0tf2+tf2t0+tft−t0+2tft2+t3,Φ4t=1tf−t02−t02tf+t0t0+2tft−2t0+tft2+t3.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 Ft,yt,ẏt,y¨t, the constrained TPBVPs are transformed into an unconstrained problem. According to (3), yt is replaced by gt; thus, the original vector differential equation Ft,yt,ẏt,y¨t is transformed into F^t,gt,ġt,g¨t, which is only a function of t, the free functions gt, and their derivatives
(9)F^t,gt,ġt,g¨t=0.As mentioned earlier,F^t,gt,ġt,g¨t is unconstrained because the boundary is represented by the switch function Φkt and the projection function pkt.After determining the switch function and the projection function, we next discuss the free functiongt.


### 2.2. 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(10)gt=ξTh,where ξ are m×1 unknown coefficients and h are m basis functions.Next, the problem domaint 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
(11)z=z0+zf−z0tf−t0t−t0↔t=t0+tf−t0zf−z0z−z0.By using (11), the derivatives of gt are computed as
(12)dkgdtk=ckξTdkhzdzk,wherek∈0,n,where c≔dz/dt=zf−z0/tf−t0.


### 2.3. Domain Discretization

For solving TPBVPs numerically, the domaint∈t0,tf 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
(13)zk=−coskπNfork=1,2,⋯,N.Thus, the new vector differential equationF^t,gt,ġt,g¨t becomes F~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
(14)Liξ=F~iz0,ξ⋮F~izd,ξ⋮F~izN,ξ.By settingLi=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)Jξ=∂F~iz0,ξ∂ξ1∂F~iz0,ξ∂ξ2⋯∂F~iz0,ξ∂ξj⋮⋮⋮∂F~izd,ξ∂ξ1∂F~izd,ξ∂ξ2⋯∂F~izd,ξ∂ξj⋮⋮⋮∂F~izN,ξ∂ξ1∂F~izN,ξ∂ξ2⋯∂F~izN,ξ∂ξj.The estimation is updated by(16)ξk+1=ξk−Δξ,where Δξ=JξTJξ−1JξTLξ.The iterative process stops when the convergence tolerance is met:(17)L2Lξ<ε,where ε is the stopping criterion that is defined by the user.Figure1 shows the outline of the TFC framework.Figure 1
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 gt 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.


## 2.1. TFC for TPBVP

In general, trajectory optimization problems are second-order TPBVP [28], which is expressed as
(2)Ft,yt,ẏt,y¨t=0subjectto:yt0=y0,ytf=yf,ẏt0=ẏ0,ẏtf=ẏf,where t0 and tf represent the initial time and the terminal time, respectively, and y0,yf,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
(3)yt=gt+Φ1tp1+Φ2tp2+Φ3tp3+Φ4tp4,where pkt is expressed as
(4)p1t,git=y0i−git0,p2t,git=yfi−gitf,p3t,git=ẏ0i−ġit0,p4t,git=ẏfi−ġitf.The algorithm to derive the term of the constrained expression is given as follows:(1)
Chooseskt, which are k linearly independent support functions(2)
Calculate switching functionsΦkt 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ΦktThe support functions are defined asskk=tk−1, according to [28]. The switching function is then obtained by solving equations. When the first switch function is activated, the equations are
(5)Φ1t0=1,Φ1tf=0,Φ1′t0=0,Φ1′tf=0.The equations of the first switch functions are combined into the matrix form as(6)1t0t02t031tftf2tf3012t03t02012tf3tf2α11α21α31α41=1000.Similarly, application of this method to the other three switch functions yields the following matrix:(7)1t0t02t031tftf2tf3012t03t02012tf3tf2α11α21α31α41α21α22α32α42α31α23α33α43α41α24α34α44=1000010000100001.The coefficientαk is obtained by matrix inversion, and the expression of the switch functions is obtained as
(8)Φ1t=1tf−t03−tf23t0−tf+6t0tft−3t0+tft2+2t3,Φ2t=1tf−t03−t02t0−3tf−6t0tft+3t0+tft2−2t3,Φ3t=1tf−t02−t0tf2+tf2t0+tft−t0+2tft2+t3,Φ4t=1tf−t02−t02tf+t0t0+2tft−2t0+tft2+t3.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 Ft,yt,ẏt,y¨t, the constrained TPBVPs are transformed into an unconstrained problem. According to (3), yt is replaced by gt; thus, the original vector differential equation Ft,yt,ẏt,y¨t is transformed into F^t,gt,ġt,g¨t, which is only a function of t, the free functions gt, and their derivatives
(9)F^t,gt,ġt,g¨t=0.As mentioned earlier,F^t,gt,ġt,g¨t is unconstrained because the boundary is represented by the switch function Φkt and the projection function pkt.After determining the switch function and the projection function, we next discuss the free functiongt.


## 2.2. 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(10)gt=ξTh,where ξ are m×1 unknown coefficients and h are m basis functions.Next, the problem domaint 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
(11)z=z0+zf−z0tf−t0t−t0↔t=t0+tf−t0zf−z0z−z0.By using (11), the derivatives of gt are computed as
(12)dkgdtk=ckξTdkhzdzk,wherek∈0,n,where c≔dz/dt=zf−z0/tf−t0.


## 2.3. Domain Discretization

For solving TPBVPs numerically, the domaint∈t0,tf 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
(13)zk=−coskπNfork=1,2,⋯,N.Thus, the new vector differential equationF^t,gt,ġt,g¨t becomes F~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
(14)Liξ=F~iz0,ξ⋮F~izd,ξ⋮F~izN,ξ.By settingLi=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)Jξ=∂F~iz0,ξ∂ξ1∂F~iz0,ξ∂ξ2⋯∂F~iz0,ξ∂ξj⋮⋮⋮∂F~izd,ξ∂ξ1∂F~izd,ξ∂ξ2⋯∂F~izd,ξ∂ξj⋮⋮⋮∂F~izN,ξ∂ξ1∂F~izN,ξ∂ξ2⋯∂F~izN,ξ∂ξj.The estimation is updated by(16)ξk+1=ξk−Δξ,where Δξ=JξTJξ−1JξTLξ.The iterative process stops when the convergence tolerance is met:(17)L2Lξ<ε,where ε is the stopping criterion that is defined by the user.Figure1 shows the outline of the TFC framework.Figure 1
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 gt 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.


## 3. 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.


### 3.1. 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:(18)ṙ=v,v̇=Tm+ag,ṁ=TIspg0,where r is the inertial position, which is normalized by the radius of the Earth R0=6378145m. v is the velocity, which is normalized by R0g0, in which g0=9.81m/s2 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=TIb, where Ib is the unit vector of the body axis satisfying
(19)Ib≡1.For most launch vehicles, the mass flow is uncontrollable; thus, the thrust magnitudeT=T is constant and uncontrollable during the same flight phase. The gravity acceleration ag=−r/r3, where r=r. The specific impulse of the engine is denoted by Isp. The differentiation of the equations in (18) is with respect to dimensionless time normalized by R0/g0.As the mass flow is constant, the optimal ascent problem is treated as a minimum-time problem:(20)minJ=tf,subject to(21)ṙ=v,v̇=Tm+ag,ṁ=TIspg0,r0=r0,v0=v0,m0=m0,rtf=rf,vtf=vf.


### 3.2. First-Order Necessary Conditions

On the basis of the Pontryagin Minimum Principle, the mass of the vehicle is treated as a prescribed function of time instead of a state variable. The Hamiltonian function is written as(22)H=1+λrTv+λvTTmIb−rr3,where Ib=−λv/λv is called the primer’s vector, according to Lawden’s theory [31]. Thus, (22) is rewritten as
(23)H=1+λrTv−Tmλv−λvTrr3.The first-order necessary conditions for optimality then give the differential equations of the costate variables:(24)λ̇r=−∂H∂r=λvTrr5,λ̇v=−∂H∂v=−λr.The transversality condition is expressed as(25)Htf=0.


## 3.1. 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:(18)ṙ=v,v̇=Tm+ag,ṁ=TIspg0,where r is the inertial position, which is normalized by the radius of the Earth R0=6378145m. v is the velocity, which is normalized by R0g0, in which g0=9.81m/s2 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=TIb, where Ib is the unit vector of the body axis satisfying
(19)Ib≡1.For most launch vehicles, the mass flow is uncontrollable; thus, the thrust magnitudeT=T is constant and uncontrollable during the same flight phase. The gravity acceleration ag=−r/r3, where r=r. The specific impulse of the engine is denoted by Isp. The differentiation of the equations in (18) is with respect to dimensionless time normalized by R0/g0.As the mass flow is constant, the optimal ascent problem is treated as a minimum-time problem:(20)minJ=tf,subject to(21)ṙ=v,v̇=Tm+ag,ṁ=TIspg0,r0=r0,v0=v0,m0=m0,rtf=rf,vtf=vf.


## 3.2. First-Order Necessary Conditions

On the basis of the Pontryagin Minimum Principle, the mass of the vehicle is treated as a prescribed function of time instead of a state variable. The Hamiltonian function is written as(22)H=1+λrTv+λvTTmIb−rr3,where Ib=−λv/λv is called the primer’s vector, according to Lawden’s theory [31]. Thus, (22) is rewritten as
(23)H=1+λrTv−Tmλv−λvTrr3.The first-order necessary conditions for optimality then give the differential equations of the costate variables:(24)λ̇r=−∂H∂r=λvTrr5,λ̇v=−∂H∂v=−λr.The transversality condition is expressed as(25)Htf=0.


## 4. 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:(26)ṙ=v,(27)v̇=−Tmλvλv−rr3,(28)ṁ=TIspg0,(29)λ̇r=−∂H∂r=λvTrr5,(30)λ̇v=−∂H∂v=−λr,(31)Htf=1+λrTtfvtf−Tmλvtf−λvTtfrtfr3tf=0,where (26) and (27) are subject to
(32)r0=r0,v0=v0,m0=m0,rtf=rf,vtf=vf.Additionally, there is a redundant equation and can be removed by the TFC constraints. The derivative ofrt is exactly the function vt because TFC constraints are analytical expressions; thus, (26) can be disregarded. The problem is now reduced. The new equations are expressed as
(33)Lat=a−λ̇r−λvTrr5,Lrt=λ̇r−λvTrr5,Lvt=λ̇v+λr,LH=1+λrTv−Tmλv−λvTrr3.To solve the above equations, the TFC constraints withrt,λrt,λvt need to be constructed by the TFC method. The unknown coefficients in the TFC constraint expressions are expressed as ξa,ξλr,ξλv, and vt,at,λ̇rt,λ̇vt can be obtained by taking the derivative of the TFC constraint expression rt,λrt,λvt, respectively.The initial and terminal constraints of the problem discussed in this paper are position and velocity constraints, respectively. The TFC constraint expressions ofrt,vt,at are written as
(34)rt=gt+Φ1tp1+Φ2tp2+Φ3tp3+Φ4tp4,vt=gt+Φ̇1tp1+Φ̇2tp2+Φ̇3tp3+Φ̇4tp4,at=gt+Φ¨1tp1+Φ¨2tp2+Φ¨3tp3+Φ¨4tp4,where the projection function is written as
(35)p1t=r0−gt0,p2t=rf−gtf,p3t=v0−ġt0,p4t=vf−ġtf.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 tf. Combined with the TFC method, the new unknown variable ξt is used to represent the parameter c, and the optimal time is obtained by solving ξt. In addition, to ensure the solved final time is positive, b is used instead of b, where b2=c. Then, rt is rewritten as
(36)rt=rz=Φ1zr0+Φ2zrf+Φ3zv0b2+Φ4zvfb2+hz−Φ1zh0−Φ2zhf−Φ3zh0′−Φ4zhf′Tξa,vt=b2vz=b2Φ1′zr0+Φ2′zrf+Φ3′zv0b2+Φ4′zvfb2+h′z−Φ1′zh0−Φ2′zhf−Φ3′zh0′−Φ4′zhf′Tξa,at=b4rz=b4Φ1″zr0+Φ2″zrf+Φ3″zv0b2+Φ4″zvfb2+h″z−Φ1″z0−Φ2″zhf−Φ3″zh0′−Φ4″zhf′Tξa,where h0=hz0andhf=hzf.The expression of the switch function is similar to (8):
(37)Φ1z=1zf−z03−zf23z0−zf+6z0zfz−3z0+zfz2+2z3,Φ2z=1zf−z03−z02z0−3zf−6z0zfz+3z0+zfz2−2z3,Φ3z=1zf−z02−z0zf2+zf2z0+zfz−z0+2zfz2+z3,Φ4z=1zf−z02−z02zf+z0z0+2zfz−2z0+zfz2+z3.Next, the TFC constraint expressions of the costates are constructed by the same steps as above:(38)λrt=λrz=hzTξr,λ̇rt=b2λr′z=b2h′zTξr,λvt=λvz=hzTξv,λ̇vt=b2λv′z=b2h′zTξv.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
(39)ξ=ξa1Tξa2Tξa3Tξλr1Tξλr2Tξλr3Tξλv1Tξλv2Tξλv3TξtT∈ℝ9m+1.The loss function can be expressed as(40)L=La1TLa2TLa3TLλr1TLλr2TLλr3TLλv1TLλv2TLλv3TLHT.To use the nonlinear least-square method, the partial derivative of the loss function needs to be calculated as(41)at=b4rt=b4Φ1′′zr0+Φ2′′zrf+Φ3′′zv0b2+Φ4′′zvfb2+h″z−Φ1′′zh0−Φ2′′zhf−Φ3′′zh0′−Φ4′′zhf′Tξa.The terms of (41) are defined by
(42)Ja,ξa=∂La1∂ξa1∂La1∂ξa2∂La1∂ξa3∂La2∂ξa1∂La2∂ξa2∂La2∂ξa3∂La3∂ξa1∂La3∂ξa2∂La3∂ξa3,Ja,ξv=∂La1∂ξv1∂La1∂ξv2∂La1∂ξv3∂La2∂ξv1∂La2∂ξv2∂La2∂ξv3∂La3∂ξv1∂La3∂ξv2∂La3∂ξv3,Ja,ξt=∂La1∂ξt∂La2∂ξt∂La3∂ξtT,Jλr,ξa=∂Lr1∂ξa1∂Lr1∂ξa2∂Lr1∂ξa3∂Lr2∂ξa1∂Lr2∂ξa2∂Lr2∂ξa3∂Lr3∂ξa1∂Lr3∂ξa2∂Lr3∂ξa3,Jλr,ξr=∂Lr1∂ξr1000∂Lr2∂ξr2000∂Lr3∂ξr3,Jλr,ξv=∂Lr1∂ξv1000∂Lr2∂ξv2000∂Lr3∂ξv3,Jλr,ξb=∂Lr1∂ξt∂Lr2∂ξt∂Lr3∂ξt,Jλv,ξr=∂Lv1∂ξr1000∂Lv2∂ξr2000∂Lv3∂ξr3,Jλv,ξv=∂Lv1∂ξv1000∂Lv2∂ξv2000∂Lv3∂ξv3,Jλv,ξt=∂Lv1∂ξt∂Lv2∂ξt∂Lv3∂ξtT,JH,ξa=∂LH∂ξa1∂LH∂ξa2∂LH∂ξa3,JH,ξr=∂LH∂ξr1∂LH∂ξr2∂LH∂ξr3,JH,ξv=∂LH∂ξv1∂LH∂ξv2∂LH∂ξv3,JH,ξt=∂LH∂ξt,where all partial derivatives are provided in the Appendix.


### 4.2. Initialization

When using the iterative least-square method to calculate, some parameters need to be estimated reasonably, so that the iterative process of the algorithm can go on without violating the basic mathematical principles. The simplest initialization is to set them equal to zero. This is equivalent to connecting the boundary value problem with the simplest interpolating expression. However, theλv is related to the thrust direction; initialization of λv=0 will cause issues in the TFC method because Ib=−λv/λv and ∥Ib∥≡1. Thus, the coefficient ξv is initialized using
(43)ξv0=v0v0,ξvf=−r0r0.In addition, the first guess ofξr,ξa is set equal to zero, and for setting the first guess of ξt, an estimate of the final time tf is needed; in this paper, the initial tf is set to be 300, and the initial ξt can be expressed as
(44)ξt=2tf−t0.These are uniformly discrete according to the number of polynomialsm.


## 4.1. TPBVP in TFC Framework

As mentioned in the previous section, to find the optimal state, the following nonlinear TPBVPs must be solved:(26)ṙ=v,(27)v̇=−Tmλvλv−rr3,(28)ṁ=TIspg0,(29)λ̇r=−∂H∂r=λvTrr5,(30)λ̇v=−∂H∂v=−λr,(31)Htf=1+λrTtfvtf−Tmλvtf−λvTtfrtfr3tf=0,where (26) and (27) are subject to
(32)r0=r0,v0=v0,m0=m0,rtf=rf,vtf=vf.Additionally, there is a redundant equation and can be removed by the TFC constraints. The derivative ofrt is exactly the function vt because TFC constraints are analytical expressions; thus, (26) can be disregarded. The problem is now reduced. The new equations are expressed as
(33)Lat=a−λ̇r−λvTrr5,Lrt=λ̇r−λvTrr5,Lvt=λ̇v+λr,LH=1+λrTv−Tmλv−λvTrr3.To solve the above equations, the TFC constraints withrt,λrt,λvt need to be constructed by the TFC method. The unknown coefficients in the TFC constraint expressions are expressed as ξa,ξλr,ξλv, and vt,at,λ̇rt,λ̇vt can be obtained by taking the derivative of the TFC constraint expression rt,λrt,λvt, respectively.The initial and terminal constraints of the problem discussed in this paper are position and velocity constraints, respectively. The TFC constraint expressions ofrt,vt,at are written as
(34)rt=gt+Φ1tp1+Φ2tp2+Φ3tp3+Φ4tp4,vt=gt+Φ̇1tp1+Φ̇2tp2+Φ̇3tp3+Φ̇4tp4,at=gt+Φ¨1tp1+Φ¨2tp2+Φ¨3tp3+Φ¨4tp4,where the projection function is written as
(35)p1t=r0−gt0,p2t=rf−gtf,p3t=v0−ġt0,p4t=vf−ġtf.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 tf. Combined with the TFC method, the new unknown variable ξt is used to represent the parameter c, and the optimal time is obtained by solving ξt. In addition, to ensure the solved final time is positive, b is used instead of b, where b2=c. Then, rt is rewritten as
(36)rt=rz=Φ1zr0+Φ2zrf+Φ3zv0b2+Φ4zvfb2+hz−Φ1zh0−Φ2zhf−Φ3zh0′−Φ4zhf′Tξa,vt=b2vz=b2Φ1′zr0+Φ2′zrf+Φ3′zv0b2+Φ4′zvfb2+h′z−Φ1′zh0−Φ2′zhf−Φ3′zh0′−Φ4′zhf′Tξa,at=b4rz=b4Φ1″zr0+Φ2″zrf+Φ3″zv0b2+Φ4″zvfb2+h″z−Φ1″z0−Φ2″zhf−Φ3″zh0′−Φ4″zhf′Tξa,where h0=hz0andhf=hzf.The expression of the switch function is similar to (8):
(37)Φ1z=1zf−z03−zf23z0−zf+6z0zfz−3z0+zfz2+2z3,Φ2z=1zf−z03−z02z0−3zf−6z0zfz+3z0+zfz2−2z3,Φ3z=1zf−z02−z0zf2+zf2z0+zfz−z0+2zfz2+z3,Φ4z=1zf−z02−z02zf+z0z0+2zfz−2z0+zfz2+z3.Next, the TFC constraint expressions of the costates are constructed by the same steps as above:(38)λrt=λrz=hzTξr,λ̇rt=b2λr′z=b2h′zTξr,λvt=λvz=hzTξv,λ̇vt=b2λv′z=b2h′zTξv.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
(39)ξ=ξa1Tξa2Tξa3Tξλr1Tξλr2Tξλr3Tξλv1Tξλv2Tξλv3TξtT∈ℝ9m+1.The loss function can be expressed as(40)L=La1TLa2TLa3TLλr1TLλr2TLλr3TLλv1TLλv2TLλv3TLHT.To use the nonlinear least-square method, the partial derivative of the loss function needs to be calculated as(41)at=b4rt=b4Φ1′′zr0+Φ2′′zrf+Φ3′′zv0b2+Φ4′′zvfb2+h″z−Φ1′′zh0−Φ2′′zhf−Φ3′′zh0′−Φ4′′zhf′Tξa.The terms of (41) are defined by
(42)Ja,ξa=∂La1∂ξa1∂La1∂ξa2∂La1∂ξa3∂La2∂ξa1∂La2∂ξa2∂La2∂ξa3∂La3∂ξa1∂La3∂ξa2∂La3∂ξa3,Ja,ξv=∂La1∂ξv1∂La1∂ξv2∂La1∂ξv3∂La2∂ξv1∂La2∂ξv2∂La2∂ξv3∂La3∂ξv1∂La3∂ξv2∂La3∂ξv3,Ja,ξt=∂La1∂ξt∂La2∂ξt∂La3∂ξtT,Jλr,ξa=∂Lr1∂ξa1∂Lr1∂ξa2∂Lr1∂ξa3∂Lr2∂ξa1∂Lr2∂ξa2∂Lr2∂ξa3∂Lr3∂ξa1∂Lr3∂ξa2∂Lr3∂ξa3,Jλr,ξr=∂Lr1∂ξr1000∂Lr2∂ξr2000∂Lr3∂ξr3,Jλr,ξv=∂Lr1∂ξv1000∂Lr2∂ξv2000∂Lr3∂ξv3,Jλr,ξb=∂Lr1∂ξt∂Lr2∂ξt∂Lr3∂ξt,Jλv,ξr=∂Lv1∂ξr1000∂Lv2∂ξr2000∂Lv3∂ξr3,Jλv,ξv=∂Lv1∂ξv1000∂Lv2∂ξv2000∂Lv3∂ξv3,Jλv,ξt=∂Lv1∂ξt∂Lv2∂ξt∂Lv3∂ξtT,JH,ξa=∂LH∂ξa1∂LH∂ξa2∂LH∂ξa3,JH,ξr=∂LH∂ξr1∂LH∂ξr2∂LH∂ξr3,JH,ξv=∂LH∂ξv1∂LH∂ξv2∂LH∂ξv3,JH,ξt=∂LH∂ξt,where all partial derivatives are provided in the Appendix.


## 4.2. Initialization

When using the iterative least-square method to calculate, some parameters need to be estimated reasonably, so that the iterative process of the algorithm can go on without violating the basic mathematical principles. The simplest initialization is to set them equal to zero. This is equivalent to connecting the boundary value problem with the simplest interpolating expression. However, theλv is related to the thrust direction; initialization of λv=0 will cause issues in the TFC method because Ib=−λv/λv and ∥Ib∥≡1. Thus, the coefficient ξv is initialized using
(43)ξv0=v0v0,ξvf=−r0r0.In addition, the first guess ofξr,ξa is set equal to zero, and for setting the first guess of ξt, an estimate of the final time tf is needed; in this paper, the initial tf is set to be 300, and the initial ξt can be expressed as
(44)ξt=2tf−t0.These are uniformly discrete according to the number of polynomialsm.


## 5. 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. Table1 lists the parameters of the launch vehicle and mission in the numerical simulations.Table 1
Parameters for launch vehicle.
ParameterValueLongitude of launch point (°)110.95Latitude of launching point (°)19.61Launch azimuth (°)90The specific impulse,Isp (m/s)3365Thrust,T (N)2843599.98Mass rate,ṁ (kg/s)845.052Table2 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.Table 2
Parameters of the boundary conditions.
ParameterValueInitial position,r0 (m)[371973.739, 6493779.849, -13899.978]Initial velocity,v0 (m/s)[3652.033, 556.843, -2.666]Initial mass,m0 (kg)350306Semimajor (m)6595487Eccentricity0.0053Inclination (°)20.009Ascending node (°)17.250Argument of perigee (°)32.390True anomaly (°)77.969Terminal position,rf (m)[1912866.558, 6304148.648, 2551.256]Terminal velocity,vf (m/s)[7457.930, -2220.619, 178.661]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 Figures2 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 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, Figure 7 shows the residual of acceleration, and it can be seen that the residual of TFC is about 10−14 or less for the whole solution domain.Figure 2
Vector of velocity solution by TFC, GPOPS-II, and single shooting.Figure 3
Vector of position solution by TFC, GPOPS-II, and single shooting.Figure 4
Vector of thrust solution by TFC.Figure 5
Height solution by TFC, GPOPS-II, and single shooting.Figure 6
Velocity solution by TFC, GPOPS-II, and single shooting.Figure 7
Residual by TFC.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 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. Table3 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.Table 3
Solutions of different discrete points and polynomial numbers.
Discrete pointsState and costate polynomial numbersFirst-order optimalityCost time (s)Iteration number10010, 40.4140.156110030, 40.001380.161410060, 42.09e-060.231510080, 40.005260.291410080, 102.13e-050.422910060, 6050.3—Max5010, 40.3490.16545030, 45.94e-050.23455050, 40.0005550.2526


## 6. Conclusion

In this study, we proposed a new approach to solve the fuel-optimal 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 algorithm can be used to obtain the solution of the problem(3)
The residual of the solution is about10−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



---

*Source: 2734230-2021-09-13.xml*