$$\mathbf{M} = \arccos \frac{h\_{\mathbf{z}}}{h}, \begin{bmatrix} h\_{\mathbf{x}} \\ h\_{\mathbf{y}} \\ h\_{\mathbf{z}} \end{bmatrix} = \begin{bmatrix} \mathbf{y}V\_{\mathbf{z}} - zV\_{\mathbf{y}} \\ zV\_{\mathbf{x}} - xV\_{\mathbf{z}} \\ xV\_{\mathbf{y}} - yV\_{\mathbf{x}} \end{bmatrix}, \ h = \sqrt{h\_{\mathbf{x}}^2 + h\_{\mathbf{y}}^2 + h\_{\mathbf{z}}^2}, \tag{2.20}$$

where [*hx* , *hy* , *hz*] *T* is the vector product of velocity and position.

$$\sin\Omega = \frac{h\_x}{\sqrt{h\_x^2 + h\_y^2}}, \cos\Omega = -\frac{h\_y}{\sqrt{h\_x^2 + h\_y^2}},\tag{2.21}$$

where ∈ [0, 2π], and the quadrant is determined according to the symbols of sin and cos . Furthermore,

$$w = u - f,\tag{2.22}$$

where *f* is the true anomaly, *u* is the angular argument to the ascending node and calculated by the following equation:

$$\mu = \begin{cases} \arccos\frac{\frac{\mathbf{z}\_r \cos \Omega + \mathbf{y}\_r \sin \Omega}{r}}{r}, & (\mathbf{z}\_r \ge 0) \\\\ 2\pi - \arccos \frac{\mathbf{z}\_r \cos \Omega + \mathbf{y}\_r \sin \Omega}{r}, & (\mathbf{z}\_r < 0) \end{cases}, \tag{2.23}$$

*f* is used to characterize the position of the payload in orbit and is expressed as follows:

$$f = \begin{cases} \arccos\frac{a\left(1-e^2\right)-r}{er}, & (\nu \ge 0) \\ 2\pi - \arccos\frac{a\left(1-e^2\right)-r}{er}, & (\nu < 0) \end{cases} \text{.} \tag{2.24}$$

The terminal velocity and position are transformed into constraints of the OEs, as well as a terminal mass constraint, shown as follows:

$$\begin{cases} \left\| \begin{bmatrix} a\_s, \ e\_s, \ i\_s, \ \Omega\_s, \ w\_s, \ f\_s \end{bmatrix}^T - \begin{bmatrix} a, \ e, \ i, \ \Omega, \ w, \ f \end{bmatrix}^T \begin{pmatrix} t\_s \end{pmatrix} \right\| \le \Delta \operatorname{Orbit},\\ m \begin{pmatrix} t\_s \end{pmatrix} \le m\_{allow}, \end{cases} \tag{2.25}$$

where *ts* is the terminal time and the subscript *s* represents the nominal state at the departure time. *Orbit* is the maximum allowance of the six orbital elements, *mallo*w is the minimum allowable mass at the end of the ascent phase. (4) Objectives

The objective of the ascent guidance method of a launcher can be expressed as the weighted sum of maximizing the residual mass at the payload departure time while minimizing the terminal state deviations:

$$\min J = -m \left( t\_{\rm s} \right) + \lambda\_{\rm orbit} \left\| \left[ a\_{\rm s}, e\_{\rm s}, i\_{\rm s}, \mathfrak{Q}\_{\rm s}, w\_{\rm s}, f\_{\rm s} \right]^T - \left[ a, e, i, \mathfrak{Q}, w, f \right]^T \left( t\_{\rm s} \right) \right\|, \tag{2.26}$$

where λ*or bit* is the weight of the terminal state deviations.

In addition to ensuring a stable flight, the peak value of the hinge torque and the peak power consumption of the servomechanisms need to be minimized:

$$\min J = \lambda\_{MH} \, |MH| + \lambda\_{sv} \int\_{0}^{t\_f} |MH \times \omega\_{sv}| dt. \tag{2.27}$$

where λ*M H* , λ*s*v are the weight coefficients, ω*s*v is the angular velocity of servomechanisms.

This section constructs a complete motion model of a launch vehicle. It should be pointed out that the tracking guidance is still used when flying in the atmosphere, so the variables related to aerodynamics are interpolated according to the parameters of the nominal trajectory. The exo-atmospheric guidance methods are the focuses of the following discussion, and the rotational equations are not included in the AGMs due to the assumption that the attitude control can track the guidance commands well.

### **2.3 Exo-Atmospheric Analytical Guidance Methods**

### *2.3.1 Basic Closed-Loop Guidance Method for Long March Launch Vehicles (LMLVs)*

• The process of CLG

A closed-loop guidance method for LMLVs in a vacuum environment is summarized as follows.

Step 1: Release the expectation of the fixed-point injection, and take five OEs directly as terminal constraints.

Step 2: Find the most matching entry point according to the current state of the launcher.

The state includes the velocity, position, mass, specific impulse, and mass flow rate. The time-to-go and entry point are solved iteratively based on the above parameters, and the entry point is updated in each guidance cycle.

Step 3: Construct an optimization problem of the current guidance cycle in the orbital coordinate system (OCS).

The main features of this optimization problem are as follows:

#### 2 Autonomous Guidance Control for Ascent Flight 47

4. Make the following conversion:

$$\dot{W}\_{x1}(t) = \frac{P\_{x1}}{m(t)} = \frac{I\_{sp} \cdot \dot{m}}{m\_0 - \dot{m}t} = \frac{I\_{sp}}{\tau - t},\tag{2.28}$$

where *Px*1 is the axial thrust, *m*0 is the initial mass, and *m*˙ is the mass-flow rate. In this way, *W*˙ *x*1(*t*) is related to *Isp* and the mass flow ratio τ = *m*0/*m*˙ , rather than the parameters that are difficult to measure in real time, such as *Px*1, *m*0, and *m*˙ . It should be noted that τ can be determined by *W*˙ *x*1(0), then the analytical expression of *W*˙ *x*1(*t*) can be obtained.

5. Design an analytic expression to represent the thrust direction, i.e., the pitch and yaw commands.

Many simulations have shown that the optimal thrust direction in a vacuum environment can be approximated as a linear function of time, as follows:

$$\begin{cases} \varphi\_{cx}(t) = \tilde{\varphi} + (-k\_1 + k\_2 \cdot t) \\ \psi\_{cx}(t) = \tilde{\psi} + (-k\_3 + k\_4 \cdot t) \end{cases} . \tag{2.29}$$

6. Solve for the unknown variables in Eq. (2.29) based on the terminal velocity and position constraints, then obtain the guidance command of the current cycle for the real-time control.

Step 4: Repeat Steps 2 and 3 during each guidance cycle until the cutoff equations are met to shut down the engines.

The CLG is also known as the iterative guidance mode (IGM). The above processing assumes that the PTO lies within the performance scope of the launcher, so we can always find a matching entry point on the PTO, and the optimization problem is converted to the planning of fixed-point terminal constraints in each guidance cycle. In the following cycle, the entry point shall be updated again.

The CLG has the following advantages over the PGMs:

Seven parameters in the guidance law of Eq. (2.29), i.e., φ˜, ψ˜ , *k*1, ∼, *k*4, and *t*, need to be solved. Note that *t* represents the time-to-go, also denoted as *Tk* . The CLG

is adopted when the rocket flies out of the atmosphere, so only the thrust and gravity are considered as the external forces, and the aerodynamic drag is omitted, which makes the following analytical solution possible.

The OCS is labeled as *OE* − ξηζ , where *OE* η points from the center of the earth to the injection point, ξ *OE* η denotes the orbital plan, and the three coordinate axes follow the right-hand rule, as shown in Fig. 2.2.

Consider that the OEs are set as the terminal constraints and the flight has been out of the atmosphere, the planning problem described in Sect. 2.2 is revised in the OCS as follows:

$$\text{Objective}: J = \int\_0^{T\_k} d\tilde{t} = T\_k,\tag{2.30}$$

$$\begin{aligned} \textit{Dynamic}: \vec{X} &= \boldsymbol{f}(\boldsymbol{X}, \boldsymbol{u}, \boldsymbol{\tilde{t}}), \\ \begin{cases} \dot{V}\_{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) &= \dot{W}\_{\boldsymbol{x}1}(\boldsymbol{\tilde{t}}) \cdot \cos \boldsymbol{\varphi}^{\*}(\boldsymbol{\tilde{t}}) \cdot \cos \boldsymbol{\psi}^{\*}(\boldsymbol{\tilde{t}}) + \boldsymbol{g}\_{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) \\ \dot{V}\_{\boldsymbol{\eta}}(\boldsymbol{\tilde{t}}) &= \dot{W}\_{\boldsymbol{x}1}(\boldsymbol{\tilde{t}}) \cdot \sin \boldsymbol{\varphi}^{\*}(\boldsymbol{\tilde{t}}) \cdot \cos \boldsymbol{\psi}^{\*}(\boldsymbol{\tilde{t}}) + \boldsymbol{g}\_{\boldsymbol{\eta}}(\boldsymbol{\tilde{t}}) \\ \dot{V}\_{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) &= -\dot{W}\_{\boldsymbol{x}1}(\boldsymbol{\tilde{t}}) \cdot \sin \boldsymbol{\psi}^{\*}(\boldsymbol{\tilde{t}}) + \boldsymbol{g}\_{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) \\ \dot{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) &= V\_{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) \\ \dot{\boldsymbol{\eta}}(\boldsymbol{\tilde{t}}) &= V\_{\boldsymbol{\eta}}(\boldsymbol{\tilde{t}}) \\ \dot{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) &= V\_{\boldsymbol{\xi}}(\boldsymbol{\tilde{t}}) \end{aligned} \tag{2.31}$$

Constraints:

$$\mathbf{X} = \begin{bmatrix} V\_{\xi 0} \ V\_{\eta 0} \ V\_{\zeta 0} \ \xi\_0 \ \eta\_0 \ \xi\_0 \end{bmatrix}^T,\tag{2.32}$$

$$N\_1\left(\mathcal{X}\left(T\_k\right), T\_k\right) = 0.\tag{2.33}$$

#### 2 Autonomous Guidance Control for Ascent Flight 49

The performance index, Eq. (2.30), is fuel-efficient and can also be expressed as the shortest time for the liquid launcher considering a basically constant thrust and mass flow rate. Equation (2.32) is the initial condition, and Eq. (2.33) represents the terminal constraints, i.e., the target orbital elements.

The control variables are the directions of the thrust vectors, which can be expressed by the Euler angles, φ*cx* and ψ*cx* , as follows:

$$\mathfrak{u} = [\phi\_{c\mathbf{x}}(t), \psi\_{c\mathbf{x}}(t)]. \tag{2.34}$$

**Hypothesis 1:** A uniform gravitational field is introduced to simplify the state equations, i.e., the gravity is expressed by the average of the gravity of the current and entry point (the selection of the entry point will be introduced later):

$$\begin{cases} g\_{\xi}(\tilde{t}) = \bar{g}\_{\xi} \\ g\_{\eta}(\tilde{t}) = \bar{g}\_{\eta} \\ g\_{\zeta}(\tilde{t}) = \bar{g}\_{\zeta} \end{cases} . \tag{2.35}$$

The following Hamiltonian function is established:

$$\begin{split} H &= 1 + \lambda^{\prime} f = 1 + \lambda\_{V\xi} (\dot{W}\_{x1} \cos \phi^{\*} \cos \psi^{\*} + \bar{g}\_{\xi}) \\ &+ \lambda\_{V\eta} (\dot{W}\_{x1} \sin \phi^{\*} \cos \psi^{\*} + \bar{g}\_{\eta}) + \lambda\_{V\xi} (-\dot{W}\_{x1} \sin \psi^{\*} \\ &+ \bar{g}\_{\xi}) + \lambda\_{\xi} V\_{\xi} + \lambda\_{\eta} V\_{\eta} + \lambda\_{\xi} V\_{\xi}. \end{split} \tag{2.36}$$

To maximize the Hamiltonian function Eq. (2.36), the following conditions should be met:

$$\frac{\partial H}{\partial \phi^\*} = \dot{W}\_{x1} \cos \psi^\* \left( -\lambda\_{V\xi} \sin \phi^\* + \lambda\_{V\eta} \cos \phi^\* \right) = 0,\tag{2.37}$$

$$\frac{\partial H}{\partial \psi^\*} = \dot{W}\_{x1} \left( -\lambda\_{V\xi} \cos \phi^\* \sin \psi^\* - \lambda\_{V\eta} \sin \phi^\* \sin \psi^\* - \lambda\_{V\xi} \cos \psi^\* \right) = 0. \tag{2.38}$$

By solving the above equations, we obtain the following equations:

$$\phi^\* = \arctan \frac{\lambda\_{V\eta}}{\lambda\_{V\xi}},\tag{2.39}$$

$$
\psi^\* = -\arctan \frac{\lambda\_{V\xi}}{\lambda\_{V\xi}} \cos \phi^\*.\tag{2.40}
$$

The adjoint equations are as follows:

$$\begin{split} \dot{\lambda}\_{V\xi} &= -\frac{\partial H}{\partial V\_{\xi}} = -\lambda\_{\xi}, \dot{\lambda}\_{V\eta} = -\frac{\partial H}{\partial V\_{\mu}} = -\lambda\_{\eta}, \dot{\lambda}\_{V\xi} = -\frac{\partial H}{\partial V\_{\zeta}} = -\lambda\_{\zeta}, \\ \dot{\lambda}\_{\xi} &= -\frac{\partial H}{\partial \xi} = 0, \dot{\lambda}\_{\eta} = -\frac{\partial H}{\partial \eta} = 0, \dot{\lambda}\_{\zeta} = -\frac{\partial H}{\partial \zeta} = 0. \end{split} \tag{2.41}$$

The following solution is derived from Eq. (2.41):

$$\begin{aligned} \lambda\_{V\xi} &= \lambda\_{V\xi 0} - \lambda\_{\xi}\tilde{t}, \lambda\_{V\eta} = \lambda\_{V\eta 0} - \lambda\_{\eta}\tilde{t}, \lambda\_{V\xi} = \lambda\_{V\xi 0} - \lambda\_{\xi}\tilde{t}, \\\lambda\_{\xi} &= \lambda\_{\xi 0}, \lambda\_{\eta} = \lambda\_{\eta 0}, \lambda\_{\xi} = \lambda\_{\xi 0}. \end{aligned} \tag{2.42}$$

(1) First, only the velocity constraints are considered.

If only the terminal velocity constraints are considered and the terminal position constraints are relaxed, then

$$
\lambda\_{\xi} = \lambda\_{\eta} = \lambda\_{\zeta} = 0. \tag{2.43}
$$

Substituting Eqs. (2.43) and (2.42) into Eqs. (2.39) and (2.40), we obtain the optimal solutions of the control variables:

$$\phi^\* = \arctan \frac{\lambda\_{V\eta 0}}{\lambda\_{V\xi 0}} = \tilde{\phi},\tag{2.44}$$

$$
\psi^\* = -\arctan \frac{\lambda\_{V\xi0}}{\lambda\_{V\xi0}} \cos \tilde{\phi} = \tilde{\psi}.\tag{2.45}
$$

Thus, an important conclusion is drawn: the optimal control variables are constant if only the velocity constraints are taken into account. To determine this constant, we substitute Eqs. (2.44) and (2.45) into the first three terms of state equations given by Eq. (2.31):

$$\begin{cases} \dot{V}\_{\xi}(\tilde{t}) = \dot{W}\_{x1}(\tilde{t}) \cdot \cos \tilde{\varphi} \cdot \cos \psi + \bar{\mathfrak{g}}\_{\xi} \\ \dot{V}\_{\eta}(\tilde{t}) = \dot{W}\_{x1}(\tilde{t}) \cdot \sin \tilde{\varphi} \cdot \cos \tilde{\psi} + \bar{\mathfrak{g}}\_{\eta} \\ \dot{V}\_{\xi}(\tilde{t}) = -\dot{W}\_{x1}(\tilde{t}) \cdot \sin \tilde{\psi} + \bar{\mathfrak{g}}\_{\xi} \end{cases} \tag{2.46}$$

Assuming that the time-to-go, *Tk* , and the entry point are known, so the terminal velocity and position constraints are determined, then φ˜ and ψ˜ can be obtained by integration:

$$\begin{cases} V\_{\xi k} - V\_{\xi 0} = L \cdot \cos \tilde{\phi} \cdot \cos \tilde{\psi} + \bar{g}\_{\xi} \cdot T\_k \\\ V\_{\eta k} - V\_{\eta 0} = L \cdot \sin \tilde{\phi} \cdot \cos \tilde{\psi} + \bar{g}\_{\eta} \cdot T\_k \\\ V\_{\xi k} - V\_{\xi 0} = -L \cdot \sin \tilde{\psi} + \bar{g}\_{\xi} \cdot T\_k \end{cases} \tag{2.47}$$

where *L* = *Tk* 0 *W*˙ *x*1*dt* ˜. Then,

$$\tilde{\phi} = \arctan \frac{V\_{\eta k} - V\_{\eta 0} - \bar{\mathbf{g}}\_{\eta} \cdot T\_k}{V\_{\xi k} - V\_{\xi 0} - \bar{\mathbf{g}}\_{\xi} \cdot T\_k},\tag{2.48}$$

$$\tilde{\psi} = \arcsin \frac{-V\_{\xi k} + V\_{\xi 0} + \bar{\mathbf{g}}\_{\xi} \cdot T\_k}{L}. \tag{2.49}$$

2 Autonomous Guidance Control for Ascent Flight 51

$$L = \int\_0^{T\_k} \dot{W}\_{x1} d\tilde{t} = \int\_0^{T\_k} \frac{I\_{sp}}{\pi - \tilde{t}} d\tilde{t} = I\_{sp} \cdot \ln \frac{\pi}{\pi - T\_k}.\tag{2.50}$$

**Hypothesis 2:** For the guidance law of Eq. (2.29), φ˜ and ψ˜ are used to ensure the terminal velocity conditions, and they are determined by Eqs. (2.48) and (2.49), respectively. The expressions(−*k*1 + *k*2 · *t* ˜) and (−*k*3 + *k*4 · *t* ˜) could be used to meet the terminal position conditions.

(2) Second, the position constraints are also included.

Thus the expressions of Eqs. (2.44), (2.45) are updated as follows:

$$
\phi^\* = \tilde{\phi} - k\_1 + k\_2 t,\ \psi^\* = \tilde{\psi} - k\_3 + k\_4 t \tag{2.51}
$$

Based on the known entry point, the terminal conditions can be transformed into the following:

$$\mathbf{X}\ (T\_k) = \begin{bmatrix} V\_{\xi k} \\ V\_{\eta k} \\ V\_{\zeta k} \\ \xi\_k \\ \eta\_k \\ \xi\_k \end{bmatrix} = \begin{bmatrix} V\_k \cdot \cos \theta\_k \\ V\_k \cdot \sin \theta\_k \\ 0 \\ 0 \\ R\_k \\ 0 \end{bmatrix},\tag{2.52}$$

where *Vk* is the required injection velocity, θ*k* is the terminal velocity inclination, and *Rk* is the terminal geocentric vector diameter.

The state equation is converted to the following form:

$$\begin{cases} \ddot{\xi}(\tilde{t}) = \dot{W}\_{x1}(\tilde{t}) \cdot \cos \varphi^\*(\tilde{t}) \cdot \cos \psi^\*(\tilde{t}) + \bar{g}\_{\xi} \\ \ddot{\eta}(\tilde{t}) = \dot{W}\_{x1}(\tilde{t}) \cdot \sin \varphi^\*(\tilde{t}) \cdot \cos \psi^\*(\tilde{t}) + \bar{g}\_{\eta} \\ \ddot{\xi}(\tilde{t}) = -\dot{W}\_{x1}(\tilde{t}) \cdot \sin \psi^\*(\tilde{t}) + \bar{g}\_{\xi} \end{cases} \tag{2.53}$$

If (−*k*1 + *k*2 · *t* ˜) and (−*k*3 + *k*4 · *t* ˜) are small quantities, then

$$\begin{array}{l} \cos k\_i \approx 1, \\ \sin k\_i \approx k\_i, \\ \cos \phi^\* = \cos \tilde{\phi} + k\_1 \sin \tilde{\phi} - k\_2 \tilde{\iota} \sin \tilde{\phi}, \\ \sin \phi^\* = \sin \tilde{\phi} - k\_1 \cos \tilde{\phi} + k\_2 \tilde{\iota} \cos \tilde{\phi}, \\ \cos \psi^\* = \cos \tilde{\psi} + k\_3 \sin \tilde{\psi} - k\_4 \tilde{\iota} \sin \tilde{\psi}, \\ \sin \psi^\* = \sin \tilde{\psi} - k\_3 \cos \tilde{\psi} + k\_4 \tilde{\iota} \cos \tilde{\psi}. \end{array} \tag{2.54}$$

The expressions (−*k*1 + *k*2 · *t* ˜) and (−*k*3 + *k*4 · *t* ˜) are required not to have apparent effects on the terminal velocity, so

$$\int\_{0}^{t\_k} (-k\_1 + k\_2 \cdot t)dt = 0,\tag{2.55}$$

$$\int\_{0}^{t\_k} (-k\_3 + k\_4 \cdot t)dt = 0.\tag{2.56}$$

**Fig. 2.3** Solving of optimal entry point

The longitudinal position is mainly determined by *Tk* . Thus the transverse and normal position constraints are:

$$
\eta\_k = \eta + \upsilon\_\eta \cdot t\_k + \int\_0^{t\_k} \int\_0^t \ddot{\eta}(\tilde{t}) dt dt,\tag{2.57}
$$

$$
\zeta\_k = \zeta + v\_\zeta \cdot t\_k + \int\_0^{t\_k} \int\_0^t \ddot{\zeta}(\tilde{t}) dt dt. \tag{2.58}
$$

The four equations, i.e., Eqs. (2.55)–(2.58), are used to solve *k*1 ∼ *k*4. The integral processes are not complicated that the derivation is omitted, and interested readers can refer to [37].

(3) Finally, the time-to-go and the entry point are solved based on the geocentric angle.

*Tk* , the latest time-to-go, would be solved with the optimal entry point concurrently. This process is illustrated in Fig. 2.3.

*P*0 is the current position of the launcher. The terminal position *Pf* can be predicted according to the CLG command planned in the current cycle. According to the geocentric angle between *Pf* and the ascending node 0 + , the position *O f* on the target orbit can be determined with the same geocentric angle. If the velocities of *Pf* and *O f* are the same, *O f* is then regarded as the latest entry point, and *Tk* needs no compensation; otherwise, a correction time *t* should be found to ensure the new predicted terminal position after *Tk* + *t* is located in the PTO, as shown as *O*∗ *f* in Fig. 2.3. *O*∗ *f* is also considered to be a new entry point.

#### 2 Autonomous Guidance Control for Ascent Flight 53

According to the above analysis, *Tk* is updated by the following process:

$$\varepsilon v\_0 + L(T\_k) + \tilde{\mathcal{g}}(T\_k + \Delta t) + \frac{\partial L}{\partial t} \cdot \Delta t = f\_v(\mathbf{S}) + \frac{\partial f\_v(\mathbf{S})}{\partial t} \cdot \Delta t,\tag{2.59}$$

$$T\_k \leftarrow T\_k + \Delta t,\tag{2.60}$$

where v0 is the current velocity of the rocket, *L*(*Tk* ) denotes the apparent velocity increment, *g*˜(*Tk* + *t*)represents the gravitational effects on the velocity, *f*v denotes the velocity of position *S* on the PTO, and *S* has the same geocentric angle as the predicted terminal position of the launcher.

After *Tk* converges to a stable value during the iteration process, the orbit entry point, *O*∗ *f* , is also determined. Thus the terminal velocity and position are known, which are used by Eqs. (2.35), (2.48), (2.49), (2.57), (2.58).

At this point, all the variables in Eq. (2.29) have been solved, and the guidance law is then updated and applied for current control. The above solving process is carried out iteratively in each guidance cycle.

In the above treatment, some approximates are made, which would produce deviations. However, as the rocket approaches the entry point, the accuracy of the above processing is also continuously improved. If the flight arc is long, the gravitational effect can be processed in segments, or a high-order approximation can be substituted.

The CLG is very sensitive to thrust variations, including thrust drops. Thus, it has a certain fault-tolerance ability by adjusting the flight path in time. When performing trajectory planning at fault time *td* and taking *td* as the start time of the following flight, we obtain

$$\dot{W}\_{x1}(t\_d) = \frac{I\_{sp}}{\tau(t\_d)}, \text{ i.e., } \tau(t\_d) = \frac{I\_{sp}}{\dot{W}\_{x1}(t\_d)}. \tag{2.61}$$

That is, τ (*td* ) is updated during each planning cycle according to the apparent acceleration *W*˙ *x*1(*td* ) measured by the IMU. Then, the dropped thrust is reflected in the apparent acceleration, which causes τ (*td* ) to increase. Under the assumption that *Isp* and *Tk* remain unchanged, *L* decreases according to Eq. (2.50), so *t* increases therefrom according to Eq. (2.59), and then *Tk* increases. This means that the flight time is extended and the orbit entry point is re-determined.

An example of LM-7 is given to illustrate IGM's adaptability to thrust variations. The engines of its second stage start up with a thrust of 150 kN, and then are tuned to 180 kN within 7 s. Four failure modes are considered, and the guidance commands (the pitch Euler angles) are compared with those under nominal conditions, as shown in Fig. 2.4, where the guidance commands are all re-planned after the failure occurs. However, the premise of the fault tolerance is that there is sufficient remaining performance to reach the PTO, and the measures if the performance is greatly degraded are discussed in Sect. 2.4.

Nowadays, many enhanced algorithms have evolved from this basic method. The improvements are mainly concentrated in Step 3 of Sect. 2.3.1, i.e., revising the form

**Fig. 2.4** Adaptability of CLG to thrust reduction

of the thrust directions or correcting the terminal velocity and position constraints. These upgraded versions of CLGs are discussed in the following sections.

### *2.3.2 Evolutions of the Closed-Loop Guidance Methods*

The evolution methods provide acceptable suboptimal solutions under more complicated scenarios or constraints, which are proved to be feasible in real flights.

#### **<IP_ADDRESS> IGM Across Different Flight Phases**

The above discussion only considers single powered flight phase. However, the IGM is not always applied in the last stage. The earlier the IGM is introduced, the more robust it is to faults.

During planning, the accelerations of the next stages need to be integrated to obtain the terminal states, so the algorithm's complexity is closely related to the number of stages or segments, and each additional segment will require additional calculations and control branches in the software. However, the number of segments is also closely related to the change of the thrust-to-weight ratio, which is determined by the trajectory characteristics. Therefore, the moment when the CLG is introduced for the real-time control should be thoroughly studied.

Segmentation for continuous powered phases is applied in the LM-2F/Y8 mission. The CLG was introduced after the fairing was jettisoned. The second stage of the LM-2F was operating at that moment, which was equipped with five engines, i.e.,

one high-thrust main engine and four low-thrust swing engines. Thus, the flying was divided into two segments: all five engines were working, or only four swing engines were working after the main engine was shut down. During the working phase of the five engines, the equivalent impulse for the second flight segment is set based on the theoretical values of four swing engines; when the main engine shut down, it was updated based on the apparent acceleration measured by the IMU. Thus, the state equations were no longer continuous when facing different powered phases. An example including two powered phases (two burns) and a coasting phase in between is shown in Fig. 2.5.

The time-to-go has three components:

$$t\_k = t\_{k1} + t\_{k2} + t\_{k3}.\tag{2.62}$$

The corresponding apparent acceleration is as follows:

$$\dot{W}\_{x1}(t) = \begin{cases} \frac{I\_{sp1}}{\frac{I\_{p1}}{t\_1 - t}} \ t\_a \le t < t\_b\\ 0 \quad t\_b \le t < t\_c \quad . \end{cases} \tag{2.63}$$

$$\begin{cases} \frac{I\_{sp3}}{\frac{I\_{sp3}}{t\_3 - t}} \ t\_c \le t \le t\_d \end{cases} \tag{2.63}$$

As discussed above, τ˜3 is set as a theoretical value during the iterative computing of the first burn and then updated by the real flight data when the second burn initiates. Thus, the closed-loop guidance across different flight phases are implemented by replacing Eq. (2.28) with Eq. (2.63).

#### **<IP_ADDRESS> IGMs with Terminal Attitude Constraints**

The guidance command is realized by adjusting the attitude of the launcher. To meet the terminal velocity and position constraints simultaneously at the orbit entry point, the guidance law requires the attitude to be tuned to a certain state to satisfy the thrust vector requirements. Furthermore, the actual flight path would deviate from the planned trajectory because of the existence of interferences and model uncertainties. These lead to large dispersions between the real and nominal attitude when the payload is released, and the degree of dispersions is mainly related to the magnitude of the interferences and the flight accelerations.

The missions with terminal attitude constraints are not uncommon, and usually a reaction control system (RCS) is configured to regulate the attitude after the orbit injection. However, if the guidance method could satisfy the constraints, the RCS could be omitted to simplify the launcher, improve the reliability, and reduce the cost. This demands that the guidance method satisfies OEs and attitude constraints simultaneously only through the thrust vector control of the pitch, yaw, and cutoff sequences.

An upgraded quadratic time-to-go function expressing the thrust direction is proposed to meet the terminal attitude constraints, i.e., Eq. (2.29) is modified as follows:

$$\begin{cases} \varphi\_{cx}(t) = \tilde{\varphi} + \left( -k\_1 + k\_2 \cdot t + k\_5 \cdot t^2 \right) \\ \psi\_{cx}(t) = \tilde{\psi} + \left( -k\_3 + k\_4 \cdot t + k\_6 \cdot t^2 \right) \end{cases} \tag{2.64}$$

The two new variables, *k*5 and *k*6, can be obtained based on the terminal pitch and yaw attitude constraints. The solving of these new variables can be found in Ref. [15], and the method is applied in the LM-2F/ T3 mission.

However, the terminal attitude cannot be set arbitrarily in this way. If the angle between the thrust and the terminal velocity direction is too large, the assumption that *k*1-*k*6 are small does not hold. This is the premise for deriving an analytical guidance law by simplifying the trigonometric functions. In the next section, another prediction and correction algorithm is proposed and compared to tackle the same problem.

### *2.3.3 Prediction-Correction Iterative Guidance Method*

To avoid the singularity in the solving of the IGM when approaching the cutoff moment, the iterative calculation is terminated in advance before the engine shuts down. The variables of the guidance law then remain unchanged for the follow-up control. The errors arising therefrom should be compensated for, which is the initial purpose of the prediction-correction IGM. Its process is described as follows [38]: (1) Based on Eq. (2.29), calculate the pitch and yaw commands when the IGM is terminated, φ(*t f* 0) and ψ(*t f* 0), respectively, where *t f* 0 = *t f* − *t* , *t f* represents terminal time, and *t f* 0 is the moment when the IGM is terminated. (2) Construct the apparent acceleration model *W*˙ *x*1(*t*):

$$
\dot{W}\_{x1}(t) = \frac{I\_{sp}}{\tau(t) - t},\tag{2.65}
$$

where τ (*t*) = *m*(*t f* 0)/*m*˙(*t*), and *t* takes *t f* 0 as the starting point:

#### 2 Autonomous Guidance Control for Ascent Flight 57

$$
\dot{W}\_{x1}(t\_{f0}) = \frac{I\_{sp}}{\tau(t) - 0}.\tag{2.66}
$$

Then,

$$
\dot{W}\_{x1}(t) = \frac{I\_{sp}}{\frac{I\_{sp}}{\dot{W}\_{x1}(t\_{f0})} - t}. \tag{2.67}
$$

(3) Calculate the increases in the apparent velocity (or velocity) and position from *t f* 0 to *t f* , and then obtain the predicted terminal state **X**∗ *k* (*t f* ) under this condition, which can be expressed as a function of the following variables:

$$\mathbf{X}\_k^\*(t\_f) = f\_{\text{pts}}(\mathbf{X}\_k(t\_{f0}), \boldsymbol{\phi}(t\_{f0}), \boldsymbol{\psi}(t\_{f0}), \dot{W}\_{\text{x1}}(t)). \tag{2.68}$$

where *f*pts denotes the function to calculate the terminal state.

The thrust vector remains the same after *t f* 0.

(4) Calculate the compensation to the terminal constraints **X***k* (*t f* ) as follows:

$$
\Delta \mathbf{X}\_k(t\_f) = \mathbf{X}\_k(t\_f) - \mathbf{X}\_k^\*(t\_f). \tag{2.69}
$$

(5) Update the terminal constraints **X***k* (*t f* ):

$$\mathbf{X}\_k(t\_f) \leftarrow \mathbf{X}\_k(t\_f) + \Delta \mathbf{X}\_k(t\_f). \tag{2.70}$$

The terminal constraints, which include the velocity and position constraints, were formerly determined by the entry point, however Eq. (2.70) renews them based on error predictions, while the time-to-go remains the same. Then, the variables of guidance law is re-calculated based on Eq. (2.70) in the current guidance cycle. Although these updated terminal constraints are not strictly optimal, the simulations show that the resulting errors are acceptable.

For different application scenarios, Eq. (2.69) has various updates. Note that **X***k* (*t f* ) represents the terminal constraint compensation induced by various factors during the period from *t f* 0 (or other moments we are interested in) to *t f* . These compensations are due to the systematic errors caused by conditions such as the in-advance termination of the IGM, the deviations caused by the tracking control or the cutoff thrust, and the other processing such as the attitude regulation.

#### **<IP_ADDRESS> Direct Injection Under High Thrust**

Under high-thrust conditions before a payload is released to an orbit, the same attitude tracking error would result in large lateral or normal velocity deviations, and the disturbances and the uncertainties of the cutoff thrust also increase. All the above effects are adverse to the entry accuracy, and this is what LM-7 faces when launching cargo spacecrafts.

In general, a terminal velocity correction system (TVCS) could be installed to reduce the velocity errors after the main engines were shut down, but this auxiliary system would increase the complexity and cost of the launcher, reduce the carrying capacity and reliability. The final solution of LM-7 is to shut down two fixed engines in the last stage ahead of the cutoff of two other swing engines to reduce the overload before entering into an orbit. The time interval between these two cutoffs should not be too long to affect the launcher's performance, but a short interval would lead to a rapid time-varying thrust due to the coupling of the cutoff thrusts of four engines, and prominent variations in the guidance commands because the position constraints are very sensitive to the thrust. This leads to the increases in the attitude tracking errors and the injection deviations. To tackle this dilemma, the position constraints of the CLG are relaxed just before the cutoff of two fixed engines, letting the guidance commands quickly enter into a stable state. Although this strategy would produces systematic position errors, but a relatively accurate prediction of the terminal position is realized due to the stable guidance commands and high tracking accuracies, thus a prediction-correction scheme could be adopted before the shutdown of the two fixed engines.

The state vector at time *t f* , **X**(2) *k* (*t f* ), which indicates the scenario that only two swing engines operate while the other two fixed engines shut down, could be predicted through the motion equation of the last stage in the vacuum regime. It should be pointed out that after two fixed engines shut down, *t f* would not be updated again.

Similarly, **X***k* (*t f* ) represents the state vector predicted by the CLG if the four engines shut down simultaneously in the end. The deviation is solved as follows:

$$
\Delta \mathbf{X}\_k(t\_f) = \mathbf{X}\_k(t\_f) - \mathbf{X}\_k^{(2)}(t\_f). \tag{2.71}
$$

The above deviation is introduced into the IGM terminal constraint:

$$\mathbf{X}\_k(t\_f) \leftarrow \mathbf{X}\_k(t\_f) + \Delta \mathbf{X}\_k(t\_f). \tag{2.72}$$

The control variables of the IGM would not update when the two fixed engines shut down and remain unchanged until the end of the flight. More detailed discussion can refer to Ref. [13].

#### **<IP_ADDRESS> Error Correction of Terminal Velocity**

Even if the TVCS is configured, the process of velocity corrections is general openloop. The entry accuracy mostly depends on whether the state vectors at the cutoff time are consistent with the theoretical conditions. Owing to the influences of various disturbances and deviations (such as thrust deviations), it is most probable that the terminal state differs from the expected value. Thus, the open-loop velocity correction based on the prescribed command sequences will lead to non-negligible velocity errors, which is disadvantageous to the orbit entry accuracy. This problem can also be solved by the prediction-correction strategy, but at this time, only the terminal velocity would be corrected.

After the engines shut down, the velocity increments caused by the cutoff thrust can be predicted as follows in the target OCS:

$$\begin{cases} \Delta W\_{xcf} = \int\_{t\_f}^{t\_{cw}} \dot{W}\_{cf}(t) \cos \phi\_f^\* \cos \psi\_f^\* dt\\ \Delta W\_{ycf} = \int\_{t\_f}^{t\_{cw}} \dot{W}\_{cf}(t) \sin \phi\_f^\* \cos \psi\_f^\* dt\\ \Delta W\_{zcf} = \int\_{t\_f}^{t\_{cw}} \dot{W}\_{cf}(t) \sin \psi\_f^\* dt \end{cases} \tag{2.73}$$

where ϕ∗ *f* , ψ∗ *f* are the real pitch and yaw control commands at the cutoff moment, *W*˙ *cf* is the apparent acceleration of the cutoff thrust, *tcut* is the moment when the cutoff thrust ends, *Wxcf* , *Wycf* , *Wzcf* are the apparent velocity increments induced by the cutoff thrust under ϕ∗ *f* , ψ∗ *f* commands.

Similarly, the nominal apparent velocity increments based on the prescribed pitch and yaw angles can also be obtained as [*Wx f* , *Wy f* , *Wz f* ] *T* , and then,

$$
\Delta \mathbf{X}\_k(t\_f) = \left[\Delta W\_{xf} - \Delta W\_{xcf}, \,\Delta W\_{yf} - \Delta W\_{ycf}, \,\Delta W\_{zf} - \Delta W\_{zcf}, \,0, \,0, \,0\right]^T. \tag{2.74}
$$

This scheme plays a major role in the LM-8/Y1 mission. The deviation between the terminal pitch command of the IGM and the nominal condition was 10.5◦, and if no measures were adopted, the velocity errors caused by cutoff thrust would exceed the TVCS's correction ability, resulting in a timed shutdown of the TVCS (its working time was scheduled as 40 s). Thus, the velocity errors would not be fully compensated for, thereby affecting the injection accuracy. Benefiting from the above algorithm, the timed shutdown of the TVCS was avoided, and the accuracy of the semi-major axis was ensured and greatly improved.

#### **<IP_ADDRESS> Handling of Terminal Attitude Constraints**

Based on the prediction-correction strategy, a new approach different from that in Sect. <IP_ADDRESS> is discussed to handle terminal attitude constraints. IGM is terminated at time *t f* 0, which is close to the terminal time *t f* , and then an attitude adjustment phase is introduced to regulate the thrust vector from the current values φ(*t f* 0), ψ(*t f* 0) to the expected terminal states φ *f* , ψ *f* . Taking the pitch channel as an example, the expression of the attitude adjustment is as follows:

$$\phi(t) = \begin{cases} \frac{2(\phi\_f - \phi(t\_{f0}))}{\left(t\_{f1} - t\_{f0}\right)^2} (t - t\_{f0})^2 + \phi\_{f0}, t\_{f0} \le t \le \frac{t\_{f1} + t\_{f0}}{2} \\\\ -\frac{2(\phi\_f - \phi(t\_{f0}))}{\left(t\_{f1} - t\_{f0}\right)^2} (t - t\_{f1})^2 + \phi\_f, \frac{t\_{f1} + t\_{f0}}{2} < t \le t\_{f1} \\\\ \phi\_f, t\_{f1} < t \le t\_f \end{cases} \tag{2.75}$$

**Fig. 2.6** Comparisons of Euler angles

where *t f* 1 is the moment when the attitude adjustment ends, φ (*t*) represents the thrust vector at time *t*, and the expression ψ (*t*) in the yaw channel is similar and not repeated.

The apparent acceleration during the attitude adjustment can be expressed as

$$
\dot{W}\_{x1}(t) = \frac{I\_{sp}}{\pi - t}. \tag{2.76}
$$

Based on Eqs. (2.75) and (2.76), we can determine the terminal state as follows:

$$\mathbf{X}\_k^\*(t\_f) = f\_{\text{pts}}(\mathbf{X}\_k(t\_0), \phi(t), \psi(t), W\_{\text{x1}}(t)). \tag{2.77}$$

The terminal compensation is

$$
\Delta \mathbf{X}\_k(t\_f) = \mathbf{X}\_k(t\_f) - \mathbf{X}\_k^\*(t\_f). \tag{2.78}
$$

Three cases under different methods, i.e., the fundamental CLG (labeled as Case1\_0), the method introduced in this section (labeled as Case1\_1), and the method introduced in Sect. <IP_ADDRESS> (labeled as Case1\_2), are compared in the Fig. 2.6, where *FC X* and *PC X* represent the pitch and yaw Euler angles, respectively.

Compared with the method introduced in Sect. <IP_ADDRESS>, we can see an obvious attitude regulation process before entering into an orbit, and the guidance commands before the adjustment are more consistent with that of the fundamental method.

### **2.4 Joint Optimization of Target Orbit and Flight Path**

For most launch failures caused by a thrust drop, the engines can continue to operate without an explosion. If the engines deteriorate to a very risky level indicated by the sensed data, active shutdown is preferable to ensure flight safety. No matter what measures are taken, the carrying capacity will degrade. If still using the prescribed guidance law under this situation, whether an orbit is reachable depends on the remaining carrying capacity. Under severe failures, the propellant will be exhausted inflight, the terminal velocity and position will not be able meet the requirements of circling the Earth, and LV/payloads will crash to the ground. Thus, onboard decision making is required in the above situations to save missions.

The assumption for analytical CLGs is that we can always find an injection point on the PTO that matches the current flight states, and how to find the target orbit is not within the scope of guidance methods. However, this assumption does not hold if the PTO is beyond the performance capabilities of the rocket. A possible rescue orbit in which the rocket makes use of the remaining fuel should be found first, then the flight path should be planned or solved concurrently. A rescue orbit refers to a new target that is different from the PTO, where the payloads can enter as the starting point for the follow-up orbital transfer to avoid crashing. It has the same meaning as a parking orbit in most contexts, where the satellite can circle the Earth for many rounds. To consume as little fuel as possible by the payloads during the orbit transfer, an optimal rescue orbit becomes attractive. Under special conditions, a rescue orbit can also be a sub-orbit with a negative perigee height. The payload could not circle the Earth under this condition and should initiate the orbit transfer as soon as possible when released by the rocket.

It is difficult to find an analytical optimal rescue orbit, so a numerical method is usually adopted. This problem was first discussed in Refs. [40, 41]. In Ref. [40], the errors of different OEs were regulated through weights of element deviations in the objective. In Ref. [41], sequential optimization was conducted based on statetriggered-indices (STI), so as to gradually approach the optimal rescue solution. A convex optimization (COP) sub-problem was constructed, and its solution was taken as the initial value of the rescue planning problem. During the COP process, the geocentric angle of the injection point after the failure is estimated referring to the IGM process, then the COP sub-problem is transformed to the OCS to simplify the terminal constraints. This treatment greatly improves the calculation efficiency of the COP. Reference [39] proposed solutions for the maximum orbital radius optimization if the rescue orbit were confined to a circular orbit. However, none of the above methods can adapt to the flight scenario where a coasting phase is inserted. In Ref. [43], an autonomous mission reconfiguration algorithm considering the coasting phase was discussed to handle the typical failure modes that occur in real launchers, but the coasting orbit and the command sequences during coasting still refer to the prescribed planning results. This treatment can obtain a feasible solution when failures occur, but it may not make full use of the remaining performance. Thus the study of the multiple graded optimization (MGO) continues, while solving the MGO online is still very challenging.

**Fig. 2.7** Framework of STI method

### *2.4.1 State-Triggered-Indices (STI) Based Method for Continuous Powered Phases*

The process of the STI-based optimization is explained in Fig. 2.7 [41]. The impacts of the deviations of the OEs on the fuel consumed to correct these errors are closely related to the orbital characteristics and the launcher's current state, which exhibit strong nonlinear features. Thus, the objective in Fig. 2.7 cannot be solved directly because of the concerns on the convergence or local optimal solutions, and it is transferred to three sub-problems.

(1) For the sake of safety, the orbital height should be ensured first after a failure occurs, so the first reaction is to find a maximum height circular orbit (MCO).

If the height is less than a safety value, it means the rocket can hardly stay in any orbit, the rescue is then abandoned. If the height meets the safe threshold but is less than the perigee height of the PTO, the circular orbit is then taken as the rescue orbit (optimal circular orbit, OCO). However, if the height is much higher, it indicates that there is a certain performance margin used to adjust other OE errors. Then, the next planning is triggered.

(2) The orbital inclination and LAN are regulated while ensuring the height of the perigee to obtain the optimal eclipse orbit (OEO).

The deviations of the inclination and LAN are eliminated as much as possible while keeping the perigee height of the rescue orbit around the required value. If the rescue orbit can be coplanar with the PTO, the following planning will be triggered again.

(3) The argument of the perigee, semi-major axis, and eccentricity are regulated while maintaining the perigee height and orbital planar elements to obtain the optimal rescue orbit (ORO).

#### 2 Autonomous Guidance Control for Ascent Flight 63

During the optimization, the solution of the current sub-problem is taken as the initial guess of the next sub-problem. The initial guess can meet the equality constraints of all the motion equations, improving the convergence and efficiency of the numerical computation. However, to obtain a reasonable initial guess value of the first sub-problem, the nonlinear terminal constraints of the OEs are transformed to the OCS. This idea is inherited from the IGM, and the transformation matrix is as follows:

$$G\_O = \begin{bmatrix} -\sin\Omega\_0 \cos i\_0 \cos\Omega\_0 \cos i\_0 & \sin i\_0 \\ \cos\Omega\_0 & \sin\Omega\_0 & 0 \\ -\sin\Omega\_0 \sin i\_0 & \cos\Omega \sin i\_0 & -\cos i\_0 \end{bmatrix} \begin{bmatrix} \cos\Phi\_k - \sin\Phi\_k \ 0 \\ \sin\Phi\_k & \cos\Phi\_k \ 0 \\ 0 & 0 \end{bmatrix},\tag{2.79}$$

where *k* is the geocentric angle between the orbit entry point and the ascending node (see Fig. 2.3).

In the OCS labeled as *O* − ξηζ , the position components along the *O*ξ and *O*ζ axes are 0, and the velocity components along the *O*η and *O*ζ axes are 0. The terminal constraints are summarized as follows:

$$
\xi\_f = \xi\_f = 0, \ V\_\eta = V\_\zeta = 0, \ \mu = \eta\_f V\_{\xi f}^2. \tag{2.80}
$$

Compared with the constraints in the LICS, i.e., Eqs. (2.16)–(2.23), (2.80) is greatly simplified, where *k* in Fig. 2.3 can be predicted as follows [40]:

$$
\Phi\_k = \Phi\_0 + d\Phi,\tag{2.81}
$$

where 0 is the geocentric angle between the position of the launch vehicle at the current time and the ascending node, and *d* is the geocentric angle in the current orbital plane during the remaining flight range.

$$d\Phi \approx \frac{d\Phi\_{ref}}{\kappa},\tag{2.82}$$

where κ is the percentage of nominal thrust after a failure occurs.

$$d\Phi\_{ref} = \arccos\left(\frac{\mathbf{P}\_{ref} \cdot \mathbf{P}\_0}{\|\mathbf{P}\_{ref}\| \cdot \|\mathbf{P}\_0\|}\right). \tag{2.83}$$

An example is given below. For a PTO with *hp* = 200 km and *ha* = 300 km, the corresponding OEs are shown in Table 2.1. It is assumed that the thrust reduction is caused by the decrease in the mass flow rate at 118.2 s, and the remaining thrust is 77.94%. The results in Fig. 2.8 were obtained by the STI-based processing.

Trajref is the nominal flight path. At the fault moment, *k* is estimated as 171.05◦, and *dest* = 15.6475◦. The OCS can be established when *i*0 = 41.27◦ and 0 = 315.51◦ at the fault time. The result of the COP is shown by the blue line labeled TrajCVX, which can be taken as the initial guess for the OCO by the adaptive collo-

**Table 2.1** Orbital elements of prescribed target orbit

**Fig. 2.8** Trajectories of rescue results

**Table 2.2** Orbital elements of OCO, OEO and ORO

cation method. The OCO result is represented by TrajOCO, and *dact* = 15.6496◦, showing that the deviation from the estimated value is 0.0021◦.

Since the height of OCO is 208.8 km, which is greater than the *h p* of the PTO, then the following planning is triggered, where [λ*hp*, λ*a*, λ*i*, λ]=[10−3, 10−3, 1, 1]. The result is shown as TrajOEO with *hp* = 200.1 km and *ha* = 215.6 km.

By defining ε*i* = ε = 0.05◦, *i* and are both less than 0.05◦, triggering the next planning. The result of the ORO is represented as TrajORO, and the OEs are summarized in Table 2.2.

Because the first planning of the OCO needs to solve the COP problem to obtain the initial guess, the calculation time accounts for more than 60% of the total planning time. By taking the OCO as the initial value, the solutions of the OEO and ORO problems can converge quickly.

More detailed discussion can be found in Ref. [41].

### *2.4.2 Segmented Rescue Optimization Crossing Coasting Phase*

If a coasting phase is inserted into the optimization of the rescue orbit, the complexity of on-line planning will be further increased. Thus, a segmented planning strategy (SPS) is studied first, and the continuous solution is discussed in the next section. The SPS is similar to the IGM across the coasting phase, i.e., the coasting obit is taken as the target orbit of the first burn. Under nominal conditions, the IGM across the coasting arc does not lose optimality, because the terminal constraints of each flight phase are reachable. However, these constraints might not best match the remaining performance and guidance command sequences when a failure occurs, leading to the handover conditions between phases being unreachable.

However, the SPS relaxes the computational burden of online planning and demonstrates its effectiveness under the typical failure modes [42]. Under the background of launching satellites to the GTO using a two-stage rocket, its solutions are briefly explained as follows:

(1) Identify the fault mode first.

Three failure modes are considered. If the engine is going to explode, shut it down immediately, and let the subsequent stages make up for the performance loss of the premature cutoff. If an engine fails to start or shuts down by accident, restart it again if it has multiple burns. The restart may succeed or fail; even if it succeeds, it will make the engine unable to operate the following scheduled startup owing to the restriction on re-ignition times. However, the restarting scheme has the effect of postponing the fault moment and reducing the impact on the performance degradation. If only the thrust drops and there is no emergent risk, let the engine continue working. In the discussion in this section, it is assumed that there is no leakage, and all the remaining propellants can be utilized.

(2) Judge the flying regimes. If flying in the atmosphere and considering the landing area of the rocket debris, call the PGM for the tracing control until the fuel in the boosters is exhausted, then turn off the engines.

(3) If flying out of the atmosphere, evaluate the remaining performance by the ES-IGM algorithm [42].

If flying before coasting, first evaluate whether the prescribed coasting orbit is reachable; if flying after coasting, evaluate whether the PTO is reachable. The ES-IGM algorithm is based on the numerical integration and summarized in Ref. [42].

(4) Call the STI to optimize the new orbit and flight path if the prescribed orbit is beyond the remaining performance capabilities; otherwise, call the IGM for the guidance control.

The above algorithm is an approximate processing because the coasting phase is not optimized according to the fault state. If the coasting orbit is still reachable, the guidance command consequence is inherited during the coasting; if not, a new transfer orbit is planned by the STI method, and the triggering of the second burn is scheduled nearby the apogee of the new coasting orbit.

The discussion in Ref. [42] indicated that, if the PGM is adopted from the current point to the end or if the IGM is called only during the last burn, the LV/satellite may fall out of space under fault conditions with a high probability. In contrast, if the IGM is called as early as possible, the payload could be deployed into an orbit. This echoes the previous conclusions, the earlier the IGM is adopted inflight, the stronger the fault adaptability becomes. However, the IGM cannot guarantee a safe parking orbit, so the evaluation of the remaining performance is very important to support onboard decision making.

### *2.4.3 Multiple Graded Optimization*

The STI method specifies the minimum orbital height as a safety constraint, for example, not less than 150 km. Thus, the satellite could circle the Earth and then carry out an orbital transfer at an appropriate point. If taking the payload as the final stage of a launcher, the flight process of the LV/payload can be jointly optimized, which is the meaning of the E2E optimization. At this time, we can relax the safety restrictions on the orbital height, even plan a sub-orbit (the perigee height is negative) to increase the apogee height, and make the orbital transfer responsively when flying to the apogee. E2E optimization can reduce the propellant consumed during the orbital transfer.

With the increase in commercial launches and constellations, multiple-satellite ridesharing launches are becoming more and more common. The purpose of the MGO is to separate some payloads in advance during the coasting phase while sending the remaining payloads to the PTO if the performance of the launcher is greatly reduced.

To clearly explain the MGO problem, the trajectory planning problems of the powered-coasting-powered profiles are summarized in Table 2.3. Offline numerical optimization is applied to analyze and compare the features of IGM, autonomous coasting reconstruction (ACRC), and MGO under thrust drop failures.

In Table 2.3, *FT* is the nominal thrust, *F*1 *T* and *F*2 *T* are the nominal thrust magnitude of the 1st and the 2nd powered phase, respectively, and κ is the percentage of the remaining thrust to its nominal value; *t*0 and *t f* are defined as the fault time and the terminal time of the second stage, *t*1 and *t*2 are the engine cutoff time of the 1st powered phase and the start time of the 2nd powered phase, respectively; *tc*0, *tcf* are the initial and terminal times of the coasting phase; *t*1 max and *tcool* are the

#### 2 Autonomous Guidance Control for Ascent Flight 67

**Table 2.3** Description of typical optimization problems

maximum first burn time and the engine cooling time; *tc* max, *t*2 min are the maximum coast phase time and the minimum second burn time; *tcoast* is the standard coasting time, *m*min and *msep* are the minimum mass of the rocket and the separation mass off the payloads.

For the IGM, the coasting OEs [*ac*, *ec*,*ic*, *c*, w*c*] = *Fun* (*r* (*t*1), *V* (*t*1)) are introduced as the terminal constraints in the first powered flight phase. After entering the coasting phase, a timed schedule is applied as the startup condition of the second powered flight phase. Then, the IGM is called again to fly to the PTO.

For the ACRC method, the planning of the powered-coasting-powered profiles is optimized simultaneously while taking all payloads as a whole, so the coasting orbit will be re-planned, and there are no fixed OEs as the constraints of the first burn. The coasting time is planned onboard only considering the cryogenic propellant management and the precooling time required to restart engines. The terminal mass constraint is the same as that of the IGM.

Compared with the ACRC, the MGO method considers the solution of departing parts of the payloads during the coasting. Thus, *m* (*tc*0) = *m* (*t*1) − *msep*, and the terminal mass constraint of the second burn is reduced accordingly.

• Fault adaptability analysis

The following analysis is based on a two-stage launcher, and the launch site and PTO parameters are shown in Tables 2.4 and 2.5.

**Table 2.4** Main parameters

**Table 2.5** Parameters of PTO and launch site

We set *tcool* as 60 s, *tc* max as 850 s, and *t* 2 min as 50 s. The adaptive collocation method is used to plan the nominal trajectory of the launcher off-line, as shown in Fig. 2.9. The superscript '1st' represents the flight state during the flight phase when the side and core boosters are working, and '2nd', '3rd', and '4th' represent the first burn, coasting, and second burn of the second stage, respectively. According to the optimization results, the performance of the rocket is 5840 kg without considering the orbital height constraints of the coasting orbit.

According to the nominal trajectory, *tcoast* in the IGM is defined as 528.5 s. It is assumed that the launcher carries 10 identical satellites, each weighing 584 kg. The failure time is introduced in the time interval of 200–350 s of the first burn, and the thrust after failure occurs is represented by a factor κ. The simulation results are shown in Fig. 2.10, where S3 represents the fault adaptation range of the IGM, S2 represents the range of the ACRC, which is more than that of the IGM, and S1 represents the range of the MGO, which is greater than that of the ACRC. Method 1 and Method 2 represent the lower limits of fault adaptation ranges of the IGM and ACRC, respectively, and Method 3-1 and Method 3-2 represent the lower limits of the MGO corresponding to departing 5 or 9 satellites during coasting, respectively.

For Method 1, the IGM can only endure the thrust dropping by 10% if the failure occurs at 200 s. With the delay of the failure, the dropping tolerance increases exponentially, and 33% of the total thrust can still send the payloads to the PTO if the failure occurs at 350 s. For Method 2, if the coasting could be re-planned online, the allowable dropping thrust could be extended to 64% at 200 s and 6% at 350 s. If the failure state of the thrust were deteriorated beyond the lower limit, for example, thrust dropping to 45% at 200 s, all the payloads could not enter into the PTO. However, for Method 3-1 under this condition, the MGO could send half of the payloads to the PTO by releasing the other half during coasting, avoiding the complete loss of the mission. The more payloads released during the coasting, the more severe thrust drop failure could be endured, but fewer payloads would be sent to the PTO.

**Fig. 2.9** Parameters of the nominal flight trajectory

**Fig. 2.10** Fault adaptation ranges

Compared with the IGM, the ACRC could adaptively adjust the coasting orbital inclination and LAN by extending the flight time of the first burn, regulate the shape of the coasting orbit, and elevate its perigee height to reduce the fuel consumption of the second burn. For the MGO method, the main mission of the second burn was to elevate the apogee height by boosting the speed, so the acceleration of the second burn could be improved by releasing parts of the payloads in advance.

### • Case analysis

A test case is provided shown in Table 2.6.

If the minimum departure mass of 2625.5 kg could be determined onboard, five payloads should be separated in advance. With *msep* as 5 × 584 kg in the MGO, the optimization results are shown in Fig. 2.11.

The flight time of each phase by the MGO is shown in Table 2.7. The coasting OEs are shown in Table 2.8.

In the above analysis, the coasting orbit is optimized as a sub-orbit, and the satellites released during the coasting will inevitably crash to the ground. Another solution is to constrain a minimum safe perigee height of the coasting orbit, so the departed satellites could still circle the Earth and wait for rescue, but the number of satellites that could be put into the PTO would greatly decrease. No matter which solution was adopted, the MGO could avoid the complete loss of payloads for rideshare launches.

**Table 2.6** Failure state of test case

**Fig. 2.11** MGO planning results

**Table 2.8** Coasting orbital elements

However, the discussion in this section is based on off-line plannings, and onboard ACRC or MGO plannings are still challenging. A special issue is analyzed in Ref. [43], where the payload could still enter into the PTO by adjusting the coasting orbit.

### **2.5 Conclusions**

Ascent guidance methods are a basic, fully studied, and seemly mature technology. The off-line planning and on-line tracking strategy were widely applied in the early stage and have achieved good results; they are even still used currently. Considering the wind load relief and the restrictions of the debris landing area, the OLG or tracking guidance is still playing a major role in endo-atmospheric flight.

CLGs perform well for exo-atmospheric guidance. No structural load constraints and atmospheric disturbances are considered. The optimization problem is then simplified to obtain analytical solutions based on optimal control theories. Compared with the tracking guidance, the CLGs are more adaptive to model uncertainties and interferences, and they are capable of satisfying multiple terminal constraints such as six OEs to obtain higher injection accuracies. If mild thrust drop failures occur, they can also be taken as disturbances handled by the CLGs.

It's assumed that model uncertainties and disturbances are bounded. If the faults' effect is far beyond the limits, the CLGs no longer work. Thus, the AGMs are attracting more interest. AGMs cannot ensure an entry into the PTO, because they cannot violate physical laws under severe failures, but they may reconstruct the mission to avoid the complete loss of the payloads.

In conclusion, the AGMs need to solve the following sub-problems: (1) onboard model identification or reconstruction, which mainly occurs in the case of abnormal conditions, such as loss of thrust; (2) evaluation of the remaining performance, which is to simplify the decision making: whether to use the CLGs to the PTO or to reconstruct the mission; and (3) determination of the optimization objectives, i.e., keeping the payload in an orbit, the end-to-end planning, or the graded optimization for rideshare launches.

Although there is no perfect or groundbreaking method to solve all the above problems, the collocation method with smart initial guesses provides a strategy for complex onboard planning. The convergence is not guaranteed, but it is better than doing nothing to let the LV/payload assembly fall from space. Any solution, although not optimal, is acceptable under these failure scenarios. The study of the AGMs to reconstruct the mission is just beginning, while the study of the analytical guidance is still important because it is often the first choice of the initial guess.

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 3 Ascent Predictive Guidance for Thrust Drop Fault of Launch Vehicles Using Improved GS-MPSP**

**Xiaodong Yan and Cong Zhou**

### **3.1 Introduction**

Increasing complex space missions require launch vehicles to be with greater loadcarrying capacity, better orbit injection accuracy and higher reliability. Such demands also cause the increased complexity of the vehicle, leading to a higher probability of fault, especially for the propulsion system. To remedy this issue, an advanced and robust ascent guidance capable of fault-tolerant is critical for the success of mission. Iterative guidance method [1] (IGM) and powered explicit guidance [2] (PEG) are two commonly used methods for the ascent phase of launch vehicles. These two guidance methods work well in the nominal condition and can adapt to many offnominal conditions [3]. However, they lack of strong adaptive capacity, which cannot guarantee the reliability when the dynamic model or parameters change significantly. Alternatively, numerical approaches based on the optimal control theory may be the better choice. The existing algorithms can be divided into direct methods and indirect methods. Using the indirect methods, the guidance problem is transformed into Hamilton two-point boundary value problems [4] (TPBVP), but the solving process of this Hamilton two-point boundary value problem is complicated and highly sensitive to the initial guess. Using the direct method, the guidance problem is transformed into a nonlinear programming problem [5] (NLP). However, solving such problem is extremely computational intensive, which is difficult to meet the real-time requirement for online application.

In recent years, computational guidance has been proposed, which allows for onboard to generate guidance commands by the rapid trajectory optimization or

X. Yan (B) · C. Zhou

School of Astronautics, Northwestern Polytechnical University, 710072 Xi'an, People's Republic of China e-mail<EMAIL_ADDRESS>
C. Zhou e-mail<EMAIL_ADDRESS>
© The Author(s) 2023

Z. Song et al. (eds.), *Autonomous Trajectory Planning and Guidance Control for Launch Vehicles*, Springer Series in Astrophysics and Cosmology, https://doi.org/10.1007/978-981-99-0613-0\_3

other numerical computation [6]. Convex optimization method is the most popular one in this field [7]. The advantage of convex optimization is that the convex problem can be reliably solved by the interior point method to gain global optimality in polynomial time, regardless of the initial guess. Owing to the high efficiency, the convex optimization methods have been applied into various guidance designs, such as entry guidance [8], landing guidance [9], as well as ascent guidance [10]. A Newton-Kantorovich pseudospectral convex programming method was presented in [10] to solve the ascent trajectory planning problem. The combination of Newton-Kantorovich and pseudospectral discrete furtherly improves the computational efficiency. Similarly, Li et al. [11] presented a pseudospectral based convex optimization approach to solve the ascent guidance problem in the presence of thrust drop fault. Song et al. [12] investigated online joint optimization of the target orbit and flight trajectory for launch vehicles under a propulsion system fault. Most recently, the optimal abort orbit problem is studied in [13], which employed the SOCP approach to solve the circular abort orbit with the maximum of the radius. Through these literature, the computational guidance and online programming methods for ascent problems have been studied extensively. However, most of them are still limited in terms of the computational efficiency.

Another promising approach for online application is model predictive static programming (MPSP) [14]. Owing to featuring an explicit closed-form solution and avoiding numerical complexities of optimal control theory, this method exhibits a higher computational efficiency over convex optimization methods. In our previous work, a new developed generalized quasi-spectral MPSP (GS-MPSP) has been proposed [15], which furtherly improves the computational speed by using spectral discretization and collocation method. This new algorithm also offers the advantage of smooth control and better usability, and hence holds great promise for online application.

In this chapter, based on the GS-MPSP philosophy, an ascent guidance for the thrust drop fault of the launch vehicle is presented. For the ascent guidance problem, the thrust drop fault may lead to the flight time substantially changing compared to the nominal trajectory, and a proper final time or orbit injection point is not easy to give out. In this case, the ability to search the final time in a large range is required for the guidance algorithm. The original GS-MPSP [15] is able to address the free terminal time problem by formulating a sensitive relation between the final time and final outputs, and taking the final time as the additional variable to be determined. However, this operation is quite rough and feasible just for slightly adjusting the terminal time when a proper initial guess is provided. Therefore, an improved GS-MPSP (IGS-MPSP) method is first proposed, which furtherly enhances the convergence robustness for the free final time in the presence of the poor initial guess by introducing a scale factor of time interval. Then, the ascent guidance is systematically formulated using this improved algorithm. Consequently, a numerical simulation for the case of injecting into a circular orbit is conducted to verify the effectiveness of the proposed method. The comparison with the SOCP based method is also carried out.

### **3.2 Generic Theory of the IGS-MPSP Method**

The GS-MPSP method is proposed for efficiently solving a class of nonlinear terminal constraint problem. The considered nonlinear system dynamics is as follows:

$$X(t) = f(X, U, t), t \in [t\_0, t\_f] \tag{3.1}$$

$$Y(t) = h(X, t) \tag{3.2}$$

where *X* ∈ *Rn*, *U* ∈ *Rm* and *Y* ∈ *Rp* denote the state, control and output vectors, respectively. The purpose of this approach is to find suitable control history *U*(*t*) to ensure that the final system output *Yf*(*t f* ) approaches the desired value *Yd* with minimum control effort.

### *3.2.1 The Sensitivity Relation for Free-Terminal Time Continuous System*

In this subsection, a sensitivity relation for the continuous system with the free -terminal time is derived. It is based on the sensitivity relation developed in Ref. [15, 16]. The brief introduction is presented in the following.

In the proposed method, it is considered that the terminal time is adjusted by uniformly scaling the length of time interval [*t*0, *t f* ]. Thus, the updated terminal time can be written as:

$$t\_f^{l+1} = t\_f^l + \Delta \kappa \cdot (t\_f^l - t\_0) \tag{3.3}$$

where κ ∈ *R* is the scale factor, and the superscript *l* and *l* + 1 denotes the current step and the update step. The differential of Eq. (3.3) is given by

$$dt^{l+1} = dt^l + \Delta \kappa \cdot dt^l, t \in [t\_0, t\_f^l] \tag{3.4}$$

Note that in Eq. (3.4), the term κ · *dtl* denotes the changed time length for each infinitesimal time interval, *dt*, *t* ∈ [*t*0, *t l f* ]. Next, the final output vector of the system (3.1) can be expressed as follow by introducing a weighting matrix *W*(*t*) ∈ *Rp*×*n*

$$Y(X(t\_f)) = Y(X(t\_f)) + \int\_{t\_0}^{t\_f} \left[ W(t) \cdot (f(X, U, t) - \dot{X}) \right] dt \tag{3.5}$$

We then differentiate both sides of Eq. (3.5), it gives

$$\begin{split}dY(X(t\_f)) &= \frac{\partial Y(X(t))}{\partial X(t)} \cdot dX(t\_f) + \int\_{t\_0}^{t\_f} \left[ W(t) \cdot \frac{\partial f(X, U, t)}{\partial X(t)} \cdot \delta X(t) \right. \\ &\left. + W(t) \cdot \frac{\partial f(X, U, t)}{\partial U(t)} \cdot \delta U(t) - W(t) \cdot d\dot{X}(t) \right] dt \end{split} (3.6)$$

Note that in Eq. (3.6), *dX*(*t*) denotes the differential of *X*(*t*) taking into account the differential change of time interval, κ · *dt*, and δ*X* denotes the variation in *X* when the time interval is assumed to be fixed. They have the relationship as follow

$$dX(t) = \delta X(t) + X(t) \cdot \Delta \kappa \cdot dt \tag{3.7}$$

We conduct the first order differential of Eq.(3.7) respect to the time, it yields

$$d\dot{X}(t) = \delta \dot{X}(t) + \dot{X}(t) \cdot \Delta \kappa \tag{3.8}$$

Substituting the *dX***˙** (*t*) for Eq. (3.8) into the term of the integrand on right side of Eq. (3.6) leads to

$$\begin{split}dY(X(t\_f)) &= \frac{\partial Y(X(t))}{\partial X(t)} \cdot dX(t\_f) + \int\_0^{t\_f} \left[ W(t) \cdot \frac{\partial f(X, U, t)}{\partial X(t)} \cdot \delta X(t) \right. \\ &+ W(t) \cdot \frac{\partial f(X, U, t)}{\partial U(t)} \cdot \delta U(t) - W(t) \cdot \dot{X}(t) \cdot \Delta \kappa - W(t) \cdot \delta \dot{X}(t) \Big] dt \end{split} (3.9)$$

By integrating the last term on the right side of Eq. (3.6), we obtain

$$\begin{aligned} dY(X(t\_f)) &= \frac{\partial Y(X(t\_f))}{\partial X(t\_f)} \cdot \delta X(t\_f) - [W(t) \cdot \delta X(t)]\_{t=t\_f} + [W(t) \cdot \delta X(t)]\_{t=t\_0} \\ &+ \int\_0^{t\_f} \left[ W(t) \cdot \frac{\partial f(X, U, t)}{\partial X(t)} \cdot \delta X(t) + \dot{W}(t) \cdot \delta X(t) \right. \\ &+ W(t) \cdot \frac{\partial f(X, U, t)}{\partial U(t)} \cdot \delta U(t) - W(t) \cdot \dot{X}(t) \cdot \Delta \kappa \right] dt \end{aligned} (3.10)$$

Equation (3.10) can be further rearranged as

3 Ascent Predictive Guidance for Thrust Drop Fault of Launch Vehicles … 79

$$dY(X(t\_f)) = \left(\frac{\partial Y(X(t\_f))}{\partial X(t\_f)} - [W(t)]\_{t=t\_f}\right) \cdot \delta X(t\_f) + [W(t) \cdot \delta X(t)]\_{t=t\_0}$$

$$+ \int\_{t\_0}^{t\_f} \left[ \left(W(t) \cdot \frac{\partial f(X, U, t)}{\partial X(t)} + \dot{W}(t)\right) \cdot \delta X(t) + W(t) \cdot \frac{\partial f(X, U, t)}{\partial U(t)} \cdot \delta U(t) \right] dt$$

$$-W(t) \cdot \dot{X}(t) \cdot \Delta \kappa \Big] dt$$

Here, it is necessary to choose *W*(*t*) so that eliminates the coefficients of δ*X*(*t*) and δ*X*(*t f* ) in Eq. (3.11), which leads to the following weighting matrix dynamics with the associated boundary condition at the final time *t f* :

$$\dot{W}(t) = -W(t) \cdot \frac{\partial f(X, U, t)}{\partial X(t)}\tag{3.12}$$

$$\mathbf{W}(t\_f) = \frac{\partial \mathbf{Y}(\mathbf{X}(t\_f))}{\partial \mathbf{X}(t\_f)} \tag{3.13}$$

In Eq. (3.11), it is straight to obtain δ*X*(*t*0) = 0 with the fact that the initial conditions are specified. Then use this observation and the weighting matrix dynamic as presented in Eqs. (3.12) and (3.13), the Eq. (3.11) can be further simplified as

$$\begin{split}d\boldsymbol{Y}(\boldsymbol{X}(t\_{f})) &= \int\_{t\_{0}}^{t\_{f}} \left[\mathbf{W}(t) \cdot \frac{\partial \boldsymbol{f}(\boldsymbol{X}, \boldsymbol{U}, t)}{\partial \boldsymbol{U}(t)} \cdot \delta \boldsymbol{U}(t)\right] dt + \int\_{t\_{0}}^{t\_{f}} \left[-\mathbf{W}(t) \cdot \dot{\boldsymbol{X}}(t) \cdot \Delta \boldsymbol{\kappa}\right] dt \\ &= \int\_{t\_{0}}^{t\_{f}} \left[\mathbf{B}\_{\boldsymbol{\mathcal{S}}}(t) \cdot \delta \boldsymbol{U}(t)\right] dt + \mathbf{B}\_{\boldsymbol{\mathcal{K}}} \cdot \Delta \boldsymbol{\kappa} \end{split} \tag{3.14}$$

where

$$\mathcal{B}\_s(t) \stackrel{\triangle}{=} W(t) \cdot \frac{\partial f(X, U, t)}{\partial U(t)}\tag{3.15}$$

$$\mathcal{B}\_{\kappa} \stackrel{\triangle}{=} -\int\_{t\_0}^{t\_f} [\mathbf{W}(t) \cdot \dot{\mathbf{X}}(t)] dt \tag{3.16}$$

where *Bs*(*t*) is the sensitivity matrix that relates the error δ*U* to *dY* as per Ref. [16]. And *B*κ can be interpreted as the sensitivity matrix that relates the scale factor of time internal, κ to the final out error *dY*. Note that in here, the sensitive relation for terminal time is indirectly formulated by the scale factor of time internal, κ. Compared with the original way in GS-MPSP, this strategy is more accurate since it uniformly scales each infinitesimal time interval on [*t*0, *t f* ], rather than roughly adjust the final time.

### *3.2.2 The Mathematical of IGS-MPSP Method*

In IGS-MPSP, the control vector is represented by a weighted summation of some basic spectral functions to reduce the number of various to be optimized:

$$U(t) = \sum\_{i=1}^{N\_p} \mathcal{C}\_i P\_i(t) \tag{3.17}$$

where *C j* = [*c*1 *j*, *c*2 *j*,..., *cmj*] *T* is the coefficient vector corresponding to the *j*th spectral function. *Np* is the number of basic functions in the expression, and *P j*(*t*) is the basic spectral function. The spectral functions can be selected as Legendre series, Chebyshev series, etc.

Next, the new developed relationship as presented in Eq. (3.14) is applied to derive the GS-MPSP method for the free-time problem. Since the control history *U*(*t*), *t* ∈ [*t*0, *t f* ] has been represented by the spectral functions in Eq. (3.17), and a new scale factor of time internal, κ, is introduced to adjust the unspecified final time as in Eq. (3.3), both the coefficient vector [*C*1, *C*2,..., *CNp* ] and scale factor δκ are selected as variable to be determined.

First, the variation of the control history during the time *t* = [*t*0, *t f* ] can be obtained from Eq. (3.17):

$$d\boldsymbol{U}(t) = \sum\_{j=1}^{N\_p} d\mathbf{C}\_j \, \boldsymbol{P}\_j(t) \tag{3.18}$$

Substituting Eq. (3.18) into Eq. (3.14), it yields

$$\begin{split}d\boldsymbol{Y}\_{N} &= \int\_{0}^{t\_{f}} \left[\boldsymbol{\mathsf{B}}\_{s}(t) \cdot \sum\_{j=1}^{N\_{p}} \delta \mathbf{C}\_{j} \, \boldsymbol{P}\_{j}(t)\right] dt + \boldsymbol{\mathsf{B}}\_{\kappa} \cdot \boldsymbol{\Delta\kappa} \\ &= \sum\_{j=1}^{N\_{p}} \boldsymbol{A}\_{j} \cdot d\boldsymbol{\mathsf{C}}\_{j} + \boldsymbol{\mathsf{B}}\_{\kappa} \cdot \boldsymbol{\Delta\kappa} \end{split} \tag{3.19}$$

in which

$$A\_j = \int\_{t\_0}^{t\_f} \mathbf{B}\_s(t) \cdot P\_j(t) dt, j = 1, 2 \dots N\_p \tag{3.20}$$

In Eq. (3.19), *Aj* is the spectral sensitivity matrix as per Ref. [15], which relates the error of the coefficient of the *j*th spectral function, *dC j* , to the error of the output *dY N* . Thus, a linear formula for the error of the final output and the error of each coefficient vector as well as the scale factor is obtained. Then the desired coefficients and scale factor can be worked out by formulating a static programming problem:

The update equation for coefficient vectors can be written as

$$\mathbf{C}\_{j}^{l+1} = \mathbf{C}\_{j}^{l} - d\mathbf{C}\_{j}^{l} \tag{3.21}$$

where *Cl j* denotes as the *j*th coefficient at the current iteration (represented by superscript *l*), and *Cl*+1 *j* denotes the updated coefficient in the next iteration (represented by superscript *l* + 1). After substituting the expression of *dCl j* in Eq. (3.21) into Eq. (3.19), the error of final output can be written as

$$d\mathbf{Y}\_N = \sum\_{j=1}^{N\_p} \mathbf{F}\_j \delta \mathbf{C}\_j^l = \sum\_{j=1}^{N\_p} F\_j (\mathbf{C}\_j^l - \mathbf{C}\_j^{l+1}) = \mathbf{c}\_k - \sum\_{j=1}^{N\_p} A\_j \mathbf{C}\_j^{l+1} + \mathbf{B}\_k \cdot \Delta \boldsymbol{\kappa} \tag{3.22}$$

where

$$\mathbf{c}\_{\lambda} = \sum\_{j=1}^{N\_p} A\_j \mathbf{C}\_j^l \tag{3.23}$$

Equation (3.22) is apparently a linear equations set about *Cl*+1 *j* and κ, which contains *Np* × *m* + 1 unknowns and *p* equations where *Np* × *m* + 1 > *p*. Hence, Eq. (3.22) represents an under-constrained system and it is necessary to impose an appropriate performance index to facilitate a solution. In here, we consider to minimum the control effort, that is

$$J = \frac{1}{2} \int\_{t\_0}^{t\_f} \left[ U^{l+1}(t)^T \mathbf{R}(t) U^{l+1}(t) \right] \cdot dt + \mathcal{R}\_\kappa \Delta \kappa^2 \tag{3.24}$$

where *R*(*t*) and *R*κ are the weight matrix for control and scale factor. Substituting Eq. (3.17) into Eq. (3.24), it yields

$$J = \frac{1}{2} \int\_{t\_0}^{t\_f} \left[ \left( \sum\_{i=1}^{N\_p} \mathcal{C}\_i^{l+1} P\_i(t) \right)^T \mathcal{R}(t) \left( \sum\_{i=1}^{N\_p} \mathcal{C}\_i^{l+1} P\_i(t) \right) \right] \cdot dt + R\_\kappa \Delta \kappa^2 \quad (3.25)$$

Then, combined with the cost function (3.25) and the constraints given in Eq. (3.22), a static programming problem for the coefficient vector and scale vector can be formulated. The augmented cost function of this problem is given by

82 X. Yan and C. Zhou

$$\bar{J} = J + \lambda^T \left( dY\_N - \mathbf{c}\_{\lambda} + \sum\_{j=1}^{N\_p} A\_j \mathbf{C}\_j^{l+1} - \mathbf{B}\_{\kappa} \cdot \Delta \kappa \right) \tag{3.26}$$

where *λ* ∈ *Rp* is Lagrange multiplier. The first-order optimality conditions are

$$\frac{\partial \bar{J}}{\partial d\mathbf{C}\_j} = \sum\_{i=1}^{N\_p} \mathbf{R}\_{ij} \cdot \mathbf{C}\_j^{l+1} - A\_j \,^T \cdot \lambda = \mathbf{0},\tag{3.27}$$

$$j = 1, 2 \cdots N\_p$$
  $\text{a } \bar{I}$

∂ *J*¯ ∂κ = *R*κ · κ − *B*κ *T* · *λ* = **0** (3.28)

where

$$\mathcal{R}\_{ij} = \int\_{0}^{t\_{\parallel}'} \left[ P\_i(t) \mathcal{R}(t) P\_j(t) \right] \cdot dt \tag{3.29}$$

Thus, Eqs. (3.22), (3.27) and (3.28) make up a equation set about *Cl*+1 *j* , *λ* and δκ, which can be written as a compact form as follow

$$D \cdot X = E \tag{3.30}$$

where

$$\mathbf{D} = \begin{bmatrix} \mathbf{R}\_{11} & \cdots & \mathbf{R}\_{1N\_{p}} & -\mathbf{A}\_{1}^{T} & \mathbf{0} \\ \vdots & \ddots & \vdots & \vdots & 0 \\ \mathbf{R}\_{N\_{p}1} & \cdots & \mathbf{R}\_{N\_{p}N\_{p}} & -\mathbf{A}\_{N\_{p}}^{T} & \mathbf{0} \\ \mathbf{A}\_{1} & \cdots & \mathbf{A}\_{N\_{p}} & \mathbf{0} & -\mathbf{B}\_{\kappa} \\ \mathbf{0} & \cdots & \mathbf{0} & -\mathbf{B}\_{\kappa}^{T} & \mathbf{R}\_{\kappa} \end{bmatrix}, \tilde{\mathbf{X}} = \begin{bmatrix} \mathbf{C}\_{1}^{l+1} \\ \vdots \\ \mathbf{C}\_{N\_{p}}^{l+1} \\ \lambda \\ \Delta \kappa \end{bmatrix}, \mathbf{E} = \begin{bmatrix} \mathbf{0} \\ \vdots \\ \mathbf{0} \\ \mathbf{0} \\ \mathbf{c}\_{\lambda} - d\mathbf{Y}\_{N} \\ \mathbf{0} \end{bmatrix} (3.31)$$

The Eq. (3.30) contains *Np* × *m* + *p* + 1 unknowns (i.e., *Cl*+1 1 , *Cl*+1 2 ··· *Cl*+1 *Np* , *λ,*κ) and the same number of equations. Assuming that the matrix *D* is nonsingular, the unknown vector *X***˜** can be solved by

$$
\bar{X} = \mathcal{D}^{-1} \cdot \mathcal{E} \tag{3.32}
$$

Consequently, the updated coefficients *Cl*+1 *j* and scale factor κ are obtained from the solution of *X***˜** . Then substituting the κ into Eq. (3.3), the update terminal time can be obtained by

$$t\_f^{l+1} = t\_f^l + \Delta \kappa \cdot (t\_f^l - t\_0) \tag{3.33}$$

And substitute *Cl*+1 *j* (*j* = 1, 2 ··· *Np*) into Eq. (3.17), the updated control history at time *t* ∈ [*t*0, *t l*+1 *f* ] is eventually given by

$$\mathbf{U}^{l+1}(t) = \sum\_{i=1}^{N\_p} \mathbf{C}\_i^{l+1} P\_i(t), t \in [t\_0, t\_f^{l+1}] \tag{3.34}$$

*Remark 1* To implement the IGS-MPSP algorithm, the sensitivity matrix *Aj* , *B*κ and *Ri j* are necessary to be worked out in each iteration. The Gauss Quadrature Collocation method can be applied to efficiently compute such matrix and ensure the computational efficiency of this approach. The detailed procedure will be presented in the next subsection.

*Remark 2* Compared with the original algorithm [15], the improved method introduces a scale factor of time internal to adjust the terminal time, and accordingly a sensitive relation for this factor is derived. This way improves the accuracy of sensitive relation for terminal time and hence is able to search the final time in a wide range when a poor initial guess is provided. That is, the convergence robustness for initial guess of final time is improved.

### *3.2.3 The Computation of Sensitive Matrix by Gauss Quadrature Collocation*

In this subsection, the Gauss Quadrature Collocation method is applied to efficiently compute the sensitive matrix *Aj* , *B*κ and *Ri j* . The detailed procedure is presented as follow.

For the convenience of solving, the physical time *t* ∈ [*t*0, *t f* ] is converted to the scale time τ ∈ [−1, 1] by the following relation:

$$t \equiv t(\mathbf{r}, t\_0, t\_f) = \frac{t\_f - t\_0}{2}\mathbf{r} + \frac{t\_f + t\_0}{2} \tag{3.35}$$

Next, the collocation method is used to solve the weighting dynamic equation as presented in Eqs. (3.12) and (3.13). First, we rewrite the matrix equation (3.12) as the following vector equation with the independent variable τ :

$$\begin{aligned} \mathbf{W}\_k(\mathbf{r}) &= -\mathbf{W}\_k(\mathbf{r}) \cdot \mathbf{f}\_x(\mathbf{r}) \\ k &= 1, 2, \dots, p \end{aligned} \tag{3.36}$$

where *Wk* (τ ) denotes the *k*th row vector of matrix *W*(*t*), and *f x* (τ ) is defined by

$$f\_{\chi}(\mathbf{r}) \stackrel{\triangle}{=} \frac{\partial f(X, U, t)}{\partial X(t)} \cdot \frac{t\_f - t\_0}{2} \tag{3.37}$$

Then *N* Lagrange interpolating polynomials *Li*(τ )(*i* = 1, 2,..., *N*) are used to appropriate both sides of Eq. (3.36), by which the differential equation can be converted to a series of algeria equations at specified collocation points τ*i*(*i* = 1, 2 ··· *N*). In here, the Gauss-Lobatto type collocation is used, such as Legendre-Gauss-Lobatto (LGL), or Chebyshev-Gauss-Lobatto (CGL) series. Note that, in principle *Wk* (τ ) must satisfy Eq. (3.36) at all collocation points τ*i*(*i* = 1, 2,..., *N*). However, *Wk* is generally computed by integrating the matrix dynamics (3.36) backward from τ*N* to τ1 since the value at the final time *t f* (τ*N* ) is known. This means *Wk* (τ )|τ=τ1 is the last integral step as well as the integration result. Therefore, *Wk* (τ )|τ=τ1 is not necessary to strictly satisfy the differential equation (3.36), and we just consider the *N* − 1 collocation points τ*i*(*i* = 1, 2,..., *N*) for the according collocation equations. Consequently, the collocation equations are given in the compact form:

$$(\mathcal{D}\otimes I\_n)\cdot\Omega\_k = -f\cdot\Omega\_k\tag{3.38}$$

where *k* = [*Wk* (τ1), *Wk* (τ2), . . . *Wk* (τ*N* )] *T* ; *I n* is an *n* × *n* identity matrix and *D* ⊗ *I n* denotes the Kronecker product of *D* and *I n*; *D* ∈ R(*N*−1)×*N* is known as the differential matrix. The matrix *D* and *f* are given by

$$\mathbf{D} = \begin{vmatrix} \dot{L}\_1(\mathbf{r}\_2) \ \dot{L}\_2(\mathbf{r}\_2) \ \cdots \ \dot{L}\_N(\mathbf{r}\_2) \\ \dot{L}\_1(\mathbf{r}\_3) \ \dot{L}\_2(\mathbf{r}\_3) \ \cdots \ \dot{L}\_N(\mathbf{r}\_3) \\ \vdots \\ \dot{L}\_1(\mathbf{r}\_N) \ \dot{L}\_2(\mathbf{r}\_N) \ \cdots \ \dot{L}\_N(\mathbf{r}\_N) \end{vmatrix}, \mathbf{f} = \mathbf{d} \mathbf{iag} \begin{vmatrix} \mathbf{0}\_{n(N-1)\times n} \\ \mathbf{f}\_x^T(\mathbf{r}\_2) \\ \vdots \\ \mathbf{f}\_x^T(\mathbf{r}\_N) \end{vmatrix} . \tag{3.39}$$

Equation (3.38) can be further simplified by

$$A\Omega\_k = 0\tag{3.40}$$

where *A* = *f* + (*D* ⊗ *I n*) ∈ *R*(*N*−1)*n*×*N n*. Equation (3.40) contains (*N* − 1)*n* linear equations and the same number of unknowns (that is *Wk* (τ*i*)(*i* = 1, 2,..., *N* − 1)). Defining the unknown vector as *Xk* = *Wk* (τ1), *Wk* (τ2), . . . *Wk* (τ*N*−1) *T* , it is easy to obtain **k** = *X***k T** *W***k**(τ**N**) **T** . Next *A* is rearranged as *A* = [*AF* , *AN* ], where *AF* and *AN* denote the first (*N* − 1)*n* columns and the rest *n* columns of *A*, respectively. Using these relations, the linear equations can be further expressed as

$$A\Omega\_k = [A\_F, A\_N] \cdot \begin{bmatrix} X\_k \\ W\_k(\tau\_N)^T \end{bmatrix} = A\_F X\_k + A\_N W\_k(\tau\_N)^T = 0 \tag{3.41}$$

Assuming that the matrix *AF* is nonsingular, the *Xk* is eventually solved by

$$X\_k = -A\_F^{-1} A\_N \cdot W\_k(\tau\_N)^T \tag{3.42}$$

The solution of *Xk* gives the value of *k*th row of matrix *W*(τ ) at the collocation points τ*k* (τ1, τ2,...,τ*N*−1). By repeating the above calculation procedure for each row (*k* = 1, 2, ··· , *p*), the matrix *W*(τ ) at all the collocation points (τ1, τ2,...,τ*N* ) can be obtained.

Subsequently, the sensitivity matrix *Bs*(*t*) at collocation points can be calculated out according to Eq. (3.15):

$$\mathcal{B}\_{\boldsymbol{s}}(\tau\_{i}) = \left. \mathcal{W}(\tau\_{i}) \cdot \frac{\partial f(\boldsymbol{X}, \boldsymbol{U}, t)}{\partial \boldsymbol{U}(t)} \right|\_{t=\tau\_{i}}, i = 1, 2, \ldots, N \tag{3.43}$$

Lastly, the principle of Gaussian quadrature [16] is applied to compute the sensitive matrix *Aj* ,*B*κ and *Ri j* :

$$\begin{aligned} A\_j &= \int\_0^{t\_f} \mathbf{B}\_s(t) \cdot P\_j(t) dt = \frac{t\_f - t\_0}{2} \sum\_{i=1}^N \mathbf{B}\_s(\tau\_i) \cdot P\_j(\tau\_i) \cdot \eta\_i \\ j &= 1, 2, \dots, N\_p \end{aligned} \tag{3.44}$$

$$\mathcal{B}\_{\boldsymbol{\kappa}} \stackrel{\Delta}{=} -\int\_{t\_0}^{t\_f} \left[ \mathbf{W}(t) \cdot \dot{\mathbf{X}}(t) \right] = -\frac{t\_f - t\_0}{2} \sum\_{i=1}^{N} \mathbf{W}(\tau\_i) \cdot \dot{\mathbf{X}}(\tau\_i) \cdot \eta\_i \tag{3.45}$$

$$\begin{aligned} \mathcal{R}\_{ij} &= \int\_{t\_0}^{t\_f} \left[ P\_i(t) \mathcal{R}(t) P\_j(t) \right] \cdot dt = \frac{t\_f - t\_0}{2} \sum\_{k=1}^N P\_i(\tau\_k) \mathcal{R}(\tau\_k) P\_j(\tau\_k) \cdot \eta\_k \\\ i, j &= 1, 2, \dots, N\_p \end{aligned} \tag{3.46}$$

where η*i* is the weight coefficient of Gaussian quadrature corresponding to the collocation point τ*i* . In this way, such sensitive matrix is obtained by a set of algebraic operation at very few collocation points.

*Remark 3* In the calculation loop for each row of the matrix *W*, the matrix *A*remains unchanged and just need to be computed once, since which only upon to the given collocated points τ*i*(*i* = 1, 2,..., *N*) and *f x* (τ*i*). This feature effectively reduces the computational complexity.

*Remark 4* Since the spectral sensitivity matrix is directly worked out by the Gaussian quadrature method in Eqs. (3.42)–(3.46) and avoids heavy computational consumption produced by the numerical integration of a series of matrix differential equations, the computational efficiency is improved significantly.

### *3.2.4 The Implementation Step of IGS-MPSP*

The implementation procedure of this approach is provided in Algorithm 1. This method starts from the initial guess for spectral coefficients and terminal time. Then the final output errors and trajectory state are evaluated out. If the tolerance of the output errors is small enough, the desired control sequence is obtained. Otherwise, the corresponding sensitivity matrices are recalculated and the spectral coefficients as well as terminal time are updated. Then, the updated control history is generated, and the output errors are evaluated again. This iterative procedure is repeated until a specified criterion for the terminal output errors is met.

#### **Algorithm 1: IGS-MPSP algorithm**

	- (2.1) Compute the control history *Uk* (*t*), *t* ∈ [*t*0, *t f* ] by Eq.(3.17);

Integrate the system dynamic to obtain the trajectory state as well as the output error, *dY N* ∞;

(2.2)If *dY*∞ < δ*y*

Output the current control history *Uk* (*t*) and break the iteration;

Otherwise continue the iteration.

Since the spectral coefficients have no physical meaning, it's not straightforward to assign an initial guess with appropriate values. Therefore, the least-squares algorithm is applied to obtain the initial guess of the spectral coefficients when an initial guess of control sequence is provided.

Denoting *C***ˆ** 0 = [*C***ˆ** 0 1, *C***ˆ** 0 2,..., *C***ˆ** 0 *Np* ] as the initial guess of the spectral coefficients, as per Eq. (3.17), the control vector represented by the guess *C***ˆ** 0 at time step *tk* , *k* = 1, 2,..., *n* can be written as

$$\hat{\boldsymbol{U}}\_{k}^{0} = \sum\_{j=1}^{N\_{p}} P\_{j}(t\_{k}) \hat{\mathbf{C}}\_{j}^{0}, k = 1, 2, \dots, n \tag{3.47}$$

In matrix form, Eq. (3.47) can be written as

$$
\hat{\boldsymbol{U}}^{0} = \hat{\boldsymbol{C}}^{0}\boldsymbol{P} \tag{3.48}
$$

where *U***ˆ** 0 = [*U***ˆ** 0 1, *U***ˆ** 0 2,..., *U***ˆ** 0 *n*], and

#### 3 Ascent Predictive Guidance for Thrust Drop Fault of Launch Vehicles … 87

$$\mathbf{P} = \begin{vmatrix} P\_1(t\_1) & P\_1(t\_2) & \dots & P\_1(t\_n) \\ P\_2(t\_1) & P\_2(t\_2) & \dots & P\_2(t\_n) \\ \vdots & \vdots & \ddots & \vdots \\ P\_{N\_p}(t\_1) & P\_{N\_p}(t\_2) & \dots & P\_{N\_p}(t\_n) \end{vmatrix} \tag{3.49}$$

If the initial control guess sequence *U*0 = [*U*0 1, *U*0 2,..., *U*0 *n*] is given, the proper spectral function coefficients *C***ˆ** 0 are to be found so as to minimize |*U***ˆ** 0 − *U*0|. According to the principle of the least squares, the coefficients are estimated as

$$
\hat{\mathbf{C}}^{0} = U^{0} \mathbf{P}^{T} \cdot (\mathbf{P} \mathbf{P}^{T})^{-1} \tag{3.50}
$$

The initial guess of the spectral coefficients is obtained from Eq. (3.50).

### **3.3 The Ascent Predictive Guidance Under Thrust Drop Fault**

In this subsection, the proposed method is employed to solve the ascent guidance problem of launch vehicle under thrust drop fault. The problem formula is firstly introduced, then the detailed procedure to address this problem is presented.

### *3.3.1 Problem Formulation*

To be solved conveniently, a modified orbital inertial (MOI) coordinate system is firstly defined as follow. As shown in Fig. 3.1, *PF* is the position of launch vehicle when the fault occurs; *P F* is the projection of *PF* onto the injected orbital plane; *Pf* is the nominal injection point. Then, the origin of this modified orbital coordinate (MOC) is located at the center of the earth. The coordinate plane coincides with the injected orbital plane *O P F Pf* , in which the axis *OeY* directing to the midpoint of the arc *P F Pf* and the axis *OeX* perpendicular to the *OeY* . Lastly, the axis *OeZ* is determined by the right-hand-thread rule.

Note that the MOI coordinate can be determined by the position of launch vehicle at the time that the fault occurs and the injected orbit information (the inclination *i f* , longitude of ascending node *f* , and injection point for nominal trajectory). The relationship between the Modified orbital inertial (MOI) coordinate and the Earthcentered inertial (ECI) system is given by

$$\mathbf{X}^{MOI} = \mathbf{M}\_{ECI}^{MOI}(\boldsymbol{\Omega}\_f, \boldsymbol{i}\_f, \boldsymbol{P'}\_F) \cdot \mathbf{X}^{ECI} \tag{3.51}$$

**Fig. 3.1** Modified orbital inertial (MOI) coordinate system

where M*MOI ECI* denotes the transformation matrix from the ECI coordinate system to the MOI system.

It is considered the thrust drop fault occurs at the second stage of the launch vehicle. In this flight stage, the launch vehicle is assumed to fly out of the dense atmosphere and the aerodynamic forces can be ignored. Therefore, the three-dimensional point-mass dynamic equations of launch vehicles build in the (MOI) coordinate is given as follow:

$$\begin{cases} \dot{\mathbf{r}} = \mathbf{V} \\ \dot{\mathbf{V}} = T \cdot \mathbf{e}\_T / m \, -\mu \mathbf{r} / r^3 \\ \dot{m} = -m\_\varepsilon \end{cases} \tag{3.52}$$

where *r* = [*rx* ,*ry* ,*rz*] *T* is the position vector in the MOI coordinate system; *V* = [*Vx* , *Vy* , *Vz*] *T* is the inertial velocity vector; *T* is the thrust magnitude, which is considered to be constant; *m* is the mass of vehicle and *me* is the mass flow rate; *eT* denotes the thrust direction vector, which is generally aligned with the body longitudinal axis of the vehicle and can be given by

$$\mathbf{e}\_T = \begin{bmatrix} \cos \varphi \cos \psi , \sin \varphi \cos \psi , -\sin \psi \end{bmatrix}^T \tag{3.53}$$

where ϕ and ψ are the pitch angle and yaw angle relative to the MOI coordinate system, respectively. The dynamic equations as presented in Eqs. (3.52) and (3.53) can be written as the compact form

3 Ascent Predictive Guidance for Thrust Drop Fault of Launch Vehicles … 89

$$
\dot{\mathbf{x}} = f(\mathbf{x}, \mathbf{u}) \tag{3.54}
$$

where *x* = [*rx* ,*ry* ,*rz*, *Vx* , *Vy* , *Vz*] *T* is the state vector and *u* = [ϕ,ψ] *T* is the control vector of the system.

*Remark 5* Since the flight path angle and angle of attack of the launch vehicle is generally small in the second stage, the defined MOI and the according dynamic equation will ensure the pitch angle of the vehicle remain a small value. This can effectively reduce the nonlinearity of the controls as presented in Eq. (3.53) and improve the convergence of the algorithm. Moreover, such definition can simplify the terminal constraint to be introduced later.

### *3.3.2 Terminal Constraints*

It is assumed that the thrust fault of the launch vehicle takes place at the initial time *t*0, and the corresponding states are given by:

$$X(t\_0) = X\_0 \tag{3.55}$$

The final orbital injection time *t f* is constrained by

$$t\_f \le t\_{f,\text{max}}\tag{3.56}$$

where *t f*,max is the maximum burn time of the vehicle, which is determined by the remaining fuel and mass flow rate.

The terminal constraints of ascent guidance are determined by the orbital insertion conditions, which are generally provided by the semi-major axis *a f* , eccentricity *e f* , orbital inclination *i f* , and longitude of ascending node *f* . In here, we consider to entry into a circular orbit (*e f* = 0). Then, the first two conditions can be equivalently described by

$$\|\mathbf{r}(t\_f)\| = r\_f^\* \tag{3.57}$$

$$\left\|\mathbf{V}(t\_f)\right\| = V\_f^\* \tag{3.58}$$

$$(\mathbf{r}^T(t\_f)\mathbf{V}(t\_f) = 0\tag{3.59}$$

Moreover, in the modified orbital inertial (MOI) coordinate system, the final two orbital insertion conditions*i f* and *f* are equivalent to make the final position vector component *rz*(*t f* ) and velocity vector component *Vz*(*t f* ) to be zero. Therefore, the terminal constraints of this ascent guidance problem are defined by

$$h(\mathbf{x}(t\_f)) = \begin{bmatrix} \|\mathbf{r}(t\_f)\| = r\_f^\* \\ \|V(t\_f)\| = V\_f^\* \\ r^T(t\_f)V(t\_f) \\ r\_z(t\_f) \\ V\_z(t\_f) \end{bmatrix} = \mathbf{0} \tag{3.60}$$

Thus, the ascent guidance problem can be organized by

*P*0 : find *t f* , u(*t*), *t* ∈ [*t*0, *t f* ] subject to:

$$\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{u}(t), t) \tag{3.61}$$

$$\mathbf{x}(t\_0) = \mathbf{x}\_0 \tag{3.62}$$

$$h(\mathbf{x}(t\_f)) = 0\tag{3.63}$$

$$t\_f \le t\_{f,\text{max}}\tag{3.64}$$

### *3.3.3 Solved by the IGS-MPSP*

As introduced earlier, the proposed algorithm is able to solve the free-final time guidance problem. However, this algorithm cannot directly handle the inequality constraint as given in Eq. (3.64). Therefore, a numerical trick is additionally conducted to address this constraint.

First, the proposed method is employed to solve the problem *P*0 in which the constraint (3.64) is omitted. Then the obtained terminal time *t f* is checked. If this value is smaller than the maximum *t f*,max, it means the solution is feasible and the obtained control history can be used as the renewed commands of the launch vehicle. Otherwise, it implies the solution is infeasible for this problem. That is, the launch vehicle cannot directly entry into the required orbit. In this situation, a new guidance strategy is needed, such as entering into a new parking orbit or transfer orbit. This case is beyond the scope of this work.

The detailed implementation steps are summarized as follow. As presented in Fig. 3.2, the guidance strategy is triggered by a fault detection. Then the current state *t*0, x0 is employed as the initial state of the proposed algorithm, and the terminal time and control of nominal trajectory is used as the initial guess. Obviously, these controls and terminal time guess cannot steer the launch vehicle to well meet the required terminal states in the presence of thrust drop fault. Thus, the proposed algorithm is iteratively conducted to obtain the updated terminal time and control history. If the updated terminal time is smaller than the maximum value, then the corresponding control history is directly used as the new guidance for the launch vehicle. Otherwise, it implies that the launch vehicle cannot directly inject into the original orbit.

**Fig. 3.2** Ascent guidance strategy for launch vehicle under thrust drop fault

### **3.4 Numerical Results**

In this section, the numeric simulation is carried out to demonstrate the performance of the proposed method in term of accuracy and computational efficiency. It is considered that the fault of thrust drop occurs at the second stage flight phase of a launch vehicle. The nominal parameters of the launch vehicle are illustrated in Table 3.1. Both the initial conditions (at the fault occurring time) and target orbit's parameters are given in Table 3.2.

The proposed method is employed to solve such guidance problem when the thrust of the launch vehicle drops to 70% and 80% of the nominal value, respectively. In the algorithm implementation, the six order Legendre polynomials are used as the spectral function of control, and the Legendre-Gauss-Latto (LGL) points are selected as the collocation points in computing the spectral sensitivity matrix. In addition, the number of LGL nodes is taken as 15.

Additionally, a SOCP based method [11] is conducted in here as the comparison of the proposed method. This algorithm takes the minimum fuel consumption as the optimization object, and the classical Euler method is used to discrete the problems. The number of discretization nodes is set to be 50, which is determined by comprehensively considering the solving accuracy and efficiency.

**Table 3.1** Nominal parameters of the launch vehicle

**Table 3.2** Initial conditions and target orbit

In simulations, the terminal time and control history of nominal trajectory are used as the initial guess of the proposed method and SOCP based method. Moreover, all numerical simulations are implemented in the MATLAB 2021a environment on a personal desktop (Intel i7-8750H, 3.2 GHz). The CVX [17] optimization toolbox with SDPT3 4.0 [18] is employed as the solver of the SOCP based method.

### *3.4.1 The Results by the Proposed Method*

The proposed method reaches the required tolerance of terminal conditions by 8 iterations. Figure 3.3 depicts the altitude profiles of the trajectories with the thrust of 70% and 80% nominal value. Additionally, the nominal trajectory and the trajec-

**Fig. 3.3** Altitude profiles for nominal trajectory and the trajectories obtained by IGS-MPSP

**Fig. 3.4** Control profiles for the nominal trajectory

tories obtained by the nominal control in the case of thrust drop are also provided in Fig. 3.3. It can be seen that with the nominal control, the launch vehicle fails to enter into the target orbit under the thrust drop. The proposed method succeeds in regenerating the updated guidance commands (shown in Fig. 3.4) to steer the vehicle into the original target orbit. Specially, the trajectory states for each iteration of the proposed algorithm when the thrust is 80% of nominal value is presented in Fig. 3.5. It clearly reflects that the proposed method reaches the required orbit injection parameters by a few iterations, even if the relatively poor control and terminal time guess (nominal trajectory) are given. At the same time, the orbit injection time is considerably increased as compared to the nominal value, but still smaller than the maximum allowable value. These results demonstrate the effectiveness of the proposed method, which is able to re-plan the ascent trajectory under the thrust drop, and search the appropriate orbit injection time when a relatively accurate guess cannot not be given out.

### *3.4.2 Comparison with SOCP Method*

Furthermore, the comparison between the proposed method and SOCP based method are provided in Table 3.3, Figs. 3.6 and 3.7.

The Table 3.3 illustrates the terminal mass and terminal time achieved by the proposed algorithm and SOCP based method, respectively. It can be noted that the results obtained by the proposed algorithm are very close to that produced by SOCP method. The deviations are less than 0.025%. This means the proposed method achieves a near optimality for the fuel consumption compared to the SOCP based method. Moreover, the control histories and trajectories for the proposed method and SOCP are depicted in Figs. 3.6 and 3.7. As it can be seen, the control profiles as well as the trajectories by the proposed method are similar to those of SOCP, no matter in

**Fig. 3.5** The trajectory for each iteration in the 80% of nominal thrust

**Table 3.3** The results of final time and final mass

the case of 70% of the nominal thrust or 80% of the nominal thrust. This demonstrates the proposed method produces approximate effect of SOCP based method to solve the ascent guidance problem with the thrust drop fault.

Lastly, the computational efficiency of the proposed method and SOCP method are comparatively investigated by conducting the same simulation case. Figure 3.8 and Table 3.4 present the CPU time consumed by two methods, in which the time elapses for one iteration and the total are shown. It can be seen that, the CPU time consumed by the proposed method is almost one-sixtieth or one-seventieth of that

**Fig. 3.6** The Control profile for IGS-MPSP and SOCP method

by the SOCP method for one iteration and the total value. This result clearly demonstrates the superiority of the proposed method in the computational speed. Such highly computational efficiency is achieved by a series of careful design such as the spectral representation of control, sensitive matrix computation conducted by collocation method. Hence, the proposed method owns the great potential for online application.

### **3.5 Conclusion**

In this chapter, a predictive ascent guidance based on the IGS-MPSP for the thrust drop fault of the launch vehicle is presented. Firstly, an IGS-MPSP method is derived. Compared with the original GS-MPSP method, this approach introduces a scale factor for the time internal as the additional variable to adjust the terminal time. Then, a new sensitive relation for the final time is established. Since the accuracy of the sensitive relation is improved, the approach owns the better performance to search the appropriate final time in the presence of the poor initial guess. Hence, it is more suitable for the ascent guidance problem. Secondly, the application of the proposed method for the ascent guidance problem under the thrust drop fault is detailed introduced. The numerical simulation for a typical case and comparison with the SOCP based method are carried out. The results indicate the effectiveness of the proposed method, which generate approximate results with the SOCP method but with considerably higher computing efficiency.

**Fig. 3.7** Trajectories by various methods

**Table 3.4** The CPU time consumed by various methods

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 4 Birkhoff Pseudospectral Method and Convex Programming for Trajectory Optimization**

**Dangjun Zhao, Zhiwei Zhang, and Mingzhen Gui**

### **4.1 Introduction**

Trajectory optimization, an optimal control problem (OCP) in essence, is an important issue in many engineering applications including space missions, such as orbit insertion of launchers, orbit rescue, formation flying, etc. There exist two kinds of solving methods for OCP, i.e., indirect and direct methods. For some simple OCPs, using the indirect methods can result in analytic solutions, which are not easy to be obtained for complicated systems. Direct methods transcribe an OCPs into a finitedimensional nonlinear programming (NLP) problem via discretizing the states and the controls at a set of mesh points, which should be carefully designed via compromising the computational burden and the solution accuracy. In general, the larger number of mesh points, the more accurate solution as well as the larger computational cost including CPU time and memory [1]. There are many numerical methods have been developed for the transcription of OCPs, and the most common method is by using Pseudospectral (PS) collocation scheme [2], which is an optimal choice of mesh points in the reason of well-established rules of approximation theory [3]. Actually, there have several mature optimal control toolkits based PS methods, such as DIDO [4], GPOPS [5]. The resulting NLP problem can be solved by the well-known algorithm packages, such as IPOPT [6] or SNOPT [7]. However, these algorithms cannot obtain a solution in polynomial-time, and the resulting solution is locally

D. Zhao (B) · Z. Zhang · M. Gui

Z. Zhang e-mail<EMAIL_ADDRESS>
M. Gui e-mail<EMAIL_ADDRESS>
School of Automation in Central South University, Changsha, People's Republic of China e-mail<EMAIL_ADDRESS>
© The Author(s) 2023

Z. Song et al. (eds.), *Autonomous Trajectory Planning and Guidance Control for Launch Vehicles*, Springer Series in Astrophysics and Cosmology, https://doi.org/10.1007/978-981-99-0613-0\_4

optimal. Moreover, a good initial guess solution should be provided for complicated problems.

Convex optimization method provides a polynomial-time complexity for solving convex programming problem, and has been wildly used for solving OCP, especially for complicated space missions [8–11]. In convex optimization framework, the original OCP is converted into a second-order conic programming (SOCP) problem by using some convexification techniques and discretization on uniformly distributed points. To achieve sufficiently accurate solution, one way is to choose a large number of discretized points, which leads to extra computational burden, and an-other way is employing PS method [12–14]. In [15], the various discretization methods, including the zeroth order hold (ZOH), the first order hold (FOH), the classical Runge–Kutta fourth order integration method, and global PS methods (Chebyshev-Gauss-Lobatto, CGL, Legendre-Gauss-Radau, LGR) are considered in the convex optimization framework for solving the planetary powered landing trajectory optimization problem, and the authors con-clude that PS methods can provide more consistent solutions, which are less sensitive to the number of mesh points.

It is to note that only the cases of mesh points are considered in [15], while in some certain applications, especially aerospace applications, to achieve more accuracy, the more mesh points are required [16]. However, purely increasing will lead to ill-conditioned NLP problem [17], which is hard to solve or cannot to be solved. In order to overcome such ill-conditioned phenomenon, some well-conditioned PS methods for ordinary differential equations (ODEs) have been proposed in [17, 18], in which, Birkhoff PS method stemmed from Birkhoff interpolation [19] is introduced to solve higher-order ODEs. In [16, 19], the first-order Chebyshev Birkhoff PS (CBPS) method is proposed to transcribe general OCP into NLP, and the advantages of CBPS over other PS methods, especially for large mesh grids, is demonstrated by its application in solving an orbital transfer problem.

In this chapter, we apply PS methods using Bikrhoff polynomials to the convex optimization framework for OCPs, particularly, the first order and second-order Birkhoff PS (BPS) methods with LGL and CGL collocation schemes are applied to transcribe a class of cascaded second-order systems, which may be convex or not. The first-order BPS method in convex optimization framework is similar to that in [16], hence we focus on the PS transcription for convexified OCP by using secondorder Birkhoff polynomials. The main contributions of this paper lie in two aspects: (1) the unified matrix formulation of OCP by using BPS method within the convex optimization framework is proposed and validated; and (2) the computational performances and solution accuracies resulted from various PS schemes are extensively exploited thereby a useful conclusion that, using BPS method renders the remarkable drops in the condition number for the generated programming problem therefore lowering the computational cost.

The remainder of this chapter is organized as follows. The preliminaries of convex programming and PS method for general convexified optimal control problems are presented in Sect. 4.2, then the proposed well-conditioned PS method via Birkhoff polynomials within the convex optimization framework is detailed in Sect. 4.3. The demonstrated examples of a simple cart problem as well as a rescue orbit searching problem for validating the effective-ness and efficiency of the proposed method are presented in Sect. 4.4 followed by the conclusions in Sect. 4.5.

### **4.2 Preliminaries of Convex Programming and PS Method for Optimal Control**

In general, the free-final-time OCP can be converted into a fixed-final-time problem by normalizing the time into the domain of [0, 1], consequently, only the fixed-finaltime problem is considered in this chapter. A typical OCP is given as the following problem *G*.

$$\text{Problem}G\min \, J = f\_{t\_0}^{t\_f} \, L(\mathbf{x}, \mathbf{u})\text{d}t\tag{4.1}$$

$$\text{s.t. } \dot{\mathbf{x}} = \begin{bmatrix} \dot{\mathbf{x}}\_1 \\ \dot{\mathbf{x}}\_2 \end{bmatrix} = \begin{bmatrix} \mathbf{x}\_2 \\ \mathbf{f}(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{u}) \end{bmatrix} \stackrel{\Delta}{=} \mathbf{F}(\mathbf{x}, \mathbf{u}) \tag{4.2}$$

$$\mathbf{P}(\mathbf{x}, \mathbf{u}) \le \mathbf{0} \tag{4.3}$$

$$\mathbf{x}(t\_0) = \mathbf{x}\_0, \ \psi(\mathbf{x}(t\_f)) = \mathbf{0} \tag{4.4}$$

where the prescribed time instants*t*0 < *t f* < ∞, the state vector **x** = [**x***T* 1 , **x***T* 2 ] ∈ R2*n*, the control vector **u** ∈ R*m*, **x** ˙(**t0**) and **x** ˙(**tf**) are specified with the boundary constraints (4.4), the nonlinear path constraints in (4.3) may be various form in terms of **x** and **u**. The problem *G* is a kind of Lagrange problem, in which, the performance index (4.1) is the integral of Lagrange function [9] *L*(**x**, **u**) on time domain [*t*0, *t f* ]. According to the optimal control theory [20], other types including Mayer, Bollza, and Quadratic problems can be converted into a Lagrange problem. The nonlinear dynamics in (4.2) is a typical cascaded second-order system, which is widely used to formulate the linear motion of a spacecraft, such as the Mars landing [8], launch ascent [21], etc.

### *4.2.1 Convex Programming Method for OCP*

To solve the general problem **G** by convex programming method, the performance index (4.1) and the constraints (4.2)∼(4.4) should be the form of linear or secondorder cone (SOC), thus the convexification techniques are required for converting the nonlinear or concave constraints to the corresponding convex formulation. The most important issue lies in the conversion process is to guarantee the solution of the converted problem is still that of the original problem.

For some typical problems with the form of *G*, which are relatively simple, one can introduce some auxiliary state and (or) control variables as well as some unique relaxation to obtain the convex versions of the original problems. Such procedure is called a lossless convexification technique, which can be found in [8, 9, 22, 23]. For general OCPs, successive convex programming (SCP) method is proposed for handling more general nonconvex constraints [10, 24–26]. The main idea of SCP is to solve a nonlinear and nonconvex problems via iteratively solving a series of local convex approximate problem using linearization on the solution of last iteration. Denote the solution of the *k*th iteration to the original problem *G* as **x***k* , **u***k* , thereby the following sequential convex version

$$\begin{aligned} \text{Problem}: \text{min } J &= \boldsymbol{f}\_{\boldsymbol{t}\_{0}}^{\boldsymbol{t}\_{f}^{\boldsymbol{t}\_{f}}} \left[ L(\mathbf{z}^{k}) + \left[ L\_{\boldsymbol{z}}(\mathbf{z}^{k}) \right]^{T} (\mathbf{z} - \mathbf{z}^{k}) \right] dt \\ &\text{with } \quad \mathbf{z} = [\mathbf{x}^{T}, \mathbf{u}^{T}]^{T} \\ &\text{s.t. } \quad \dot{\mathbf{x}}\_{1} = \mathbf{x}\_{2} \\ &\dot{\mathbf{x}}\_{2} = \mathbf{a}\_{1}^{k} \mathbf{x}\_{1} + \mathbf{a}\_{2}^{k} \mathbf{x}\_{2} + \mathbf{g}^{k} \mathbf{u} + \mathbf{c}^{k} \\ &\mathbf{P}(\mathbf{z}^{k}) + \mathbf{P}\_{\boldsymbol{x}}^{k} (\mathbf{z} - \mathbf{z}^{k}) \le \mathbf{0} \\ &\qquad \qquad \qquad \mathbf{x}(t\_{0}) = \mathbf{x}\_{0} \\ &\boldsymbol{\Psi}(\mathbf{x}\_{t\_{f}}^{k}) + \boldsymbol{\Psi}\_{\boldsymbol{x}\_{t'}}^{k} (\mathbf{x}(t\_{f}) - \mathbf{x}\_{t\_{f}}^{k}) = 0 \\ &|\mathbf{z} - \mathbf{z}^{k}| \le \varepsilon \end{aligned} (4.5)$$

which is a first-order approximation of the problem *G*, and the symbols denote the following:

(1) The performance index *L*(**x**, **u**)is approximated by its Taylor first order expansion, and *Lz*(**z***k* ) = ∂*L*/∂**z**|**z**=**z***k* . Similarly the path constraints (4.3) and the boundary constraints (4.4) are first-orderly approximated, and

$$\mathbf{P}\_z^k = \partial \mathbf{P} / \partial \mathbf{z}|\_{\mathbf{z}=\mathbf{z}^\mathbf{k}}, \ \psi\_{\mathbf{x}\_{\mathbf{t}f}}^k = \partial \psi / \partial \mathbf{x}\_{\mathbf{t}\_\mathbf{f}}|\_{\mathbf{x}\_\mathbf{t}=\mathbf{x}\_\mathbf{t}^\mathbf{k}}.$$

(2) The dynamic constraints (4.2) are linearized on the trajectory **x***k* , **u***k* , and the linearizing coefficients are given by

$$\begin{aligned} \mathbf{a}\_1^k &= \partial \mathbf{f} / \partial \mathbf{x}\_1 |\_{\mathbf{x} = \mathbf{x}^k, \mathbf{u} = \mathbf{u}^k} , \mathbf{a}\_2^k = \partial \mathbf{f} / \partial \mathbf{x}\_2 |\_{\mathbf{x} = \mathbf{x}^k, \mathbf{u} = \mathbf{u}^k} \\ \mathbf{g}^k &= \partial \mathbf{f} / \partial \mathbf{u} |\_{\mathbf{x} = \mathbf{x}^k, \mathbf{u} = \mathbf{u}^k} \\ \mathbf{c}^k &= \mathbf{f} (\mathbf{x}^k, \mathbf{u}^k) - \mathbf{a}\_1^k \mathbf{x}\_1^k - \mathbf{a}\_2^k \mathbf{x}\_1^k - \mathbf{g}^k \mathbf{u}^k \end{aligned} \tag{4.6}$$

(3) The trust region constraints |**z** − **zk**| ≤ ε are introduced in problem SC to ensure its solution is still that of original problem G, and the reason lies in that the first-order approximation is only reasonable within a small neighborhood around {**xk**, **uk**}.

For convex problem SC, one can use some appropriate discretizing methods, and let the states and controls on the discrete points be decision variables henceforth a convex programming problem (or an SOCP problem), which can be efficiently solved by the classical conjugate gradient method [13] or the popular primal-dual interior method, such as ECOS [27], Mosek [28].

*Remark 1* Unlike the second-order approximations used in sequential quadratic programming (SQP) [29], the above procedure in the manner of first-order guarantees the approximations are convex and suitable be solved by efficient convex programming algorithms. In [10], the authors have presented a thoroughly theoretical analysis, and pointed out that, by introducing appropriate trust region constraints, the iterative solutions of the series problem SC will converge to that of the original problem *G*.

In this chapter, we concentrate on the PS method for problem SC, which contains various constraints and tedious bookkeeping, and may cause distractions. In fact, the emphasis of discretizing procedure is about the dynamic constraints, hence a distilled convex optimal control problem is given by

$$\text{Variables :} \mathbf{x} = [\mathbf{x}\_1^T, \mathbf{x}\_2^T]^T \in {}^{2n}, \mathbf{u} \in {}^m, \tau \in [-1, 1] \tag{4.7}$$

$$\text{ProblemC}: \min J[\mathbf{x}(\cdot), \mathbf{u}(\cdot)] = E(\mathbf{x}(-\mathbf{l}), \mathbf{x}(\mathbf{l})) \tag{4.8}$$

$$\begin{aligned} \text{s.t. } \dot{\mathbf{x}}\_1(\tau) &= s\mathbf{x}\_2(\tau) \\ \dot{\mathbf{x}}\_2(\tau) &= s \left[ \mathbf{a}\_1 \mathbf{x}\_1(\tau) + \mathbf{a}\_2 \mathbf{x}\_2(\tau) + \mathbf{g} \mathbf{u}(\tau) \right] \end{aligned} \tag{4.9}$$

$$\mathbf{e}(\mathbf{x}(-l), \mathbf{x}(l)) = \mathbf{0} \tag{4.10}$$

where the index function *E*(·) and the boundary constraints **e**(·) are reasonably supposed to be convex; the time mapping scale factor *s* = *dt*/*d*τ = (*t f* − *t*0)/2, which is resulted from that, the time dependent variable *t* ∈ [*t*0, *t f* ] is converted into the PS time τ ∈ [−1, 1] by using the mapping *t* = (*t f* − *t*0)τ/2 + (*t f* + *t*0)/2.

In problem C, the dynamic constraint (4.9) and boundary constraint (4.9) are considered, while the path constraints **P**(·) are not considered for the convenience of subsequent discussions. It is a reasonable simplification, since the PS results generated for problem C can be easily recovered to problem SC as shown in [16, 30].

### *4.2.2 PS Method for Convex Optimal Control Problem*

In PS optimal control techniques, any time signal *f* (·) can be expressed as a linear combination of a series of orthogonally Lagrange polynomials

$$f(\mathbf{r}) = \sum\_{i=0}^{N} f\_i \mathcal{L}\_i(\mathbf{r}) \tag{4.11}$$

where the spectral coefficient *fi* = *f* (τ*i*); π *N* := [τ0, τ1,...,τ*N* ]is the discrete point set such that −1 ≤ τ0 ≤ τ1 ≤···≤ τ*N*−1 ≤ τ*N* ≤ 1 and L*i* ∈ *PN* are the Lagrange orthogonally interpolating polynomials at the freedom of degree *N*, which satisfy the Kronecker relationship

$$\mathcal{L}\_j(\mathfrak{r}\_i) = \delta\_{ij} \tag{4.12}$$

Differentiating (4.11) and evaluating *f* (τ ) on the discrete points yields:

$$\dot{f}(\mathbf{r}\_i) = \sum\_{j=0}^{N} f\_j \dot{\mathcal{L}}\_j(\mathbf{r}\_i) \stackrel{\Delta}{=} \dot{\mathbf{f}} = \mathcal{D}\mathbf{f} \tag{4.13}$$

where the PS differentiation matrix (PSDM) D ∈ (*N*+1)×(*N*+1) is given by *di j* = L˙ *j*(τ*i*) for 0 ≤ *i*, *j* ≤ *N* ; and **f** = [ *f*0,..., *fN* ] *T* . It's to note that the higher-order differentiation matrix is given by [31]

$$\mathcal{D}^{(k)} = \mathcal{D}\mathcal{D} \dots \mathcal{D} = \mathcal{D}^k \tag{4.14}$$

Particularly for second differential of *f* (τ ), we have

$$
\ddot{\mathbf{f}} = \mathcal{D}^{(2)} \mathbf{f} \tag{4.15}
$$

The best choices for L*i* are the Legendre and Chebyshev polynomials [30] henceforth the Legendre and Chebyshev PS methods. When mesh points defined by τ*i* ∈ (−1, 1), the mesh points are Gauss points, while for Gauss-Labotto points, τ*i* ∈ [−1, 1] contain the two boundary points, and for Gauss-Radau points, τ*i* ∈ (−1, 1] or τ*i* ∈ [−1, 1). Since the boundary points must be constrained in many OCPs, Legend-Gauss-Labotto (LGL) and Chebyshev-Gauss-Labotto (CGL) collocation schemes are usually used in PS methods [32].

The PS transcription by using (4.13) for the general problem *G* can be found in [5, 33]. Here two types of PS transcription with matrix formulations for the convex problem *C* are presented as follows.

(1) Consider the dynamic constraints in problem *C* as a first-order system, and let

$$\begin{aligned} \mathbf{X} &= [\mathbf{x}\_0, \mathbf{x}\_1, \dots, \mathbf{x}\_N]^T = [\mathbf{x}(\tau\_0), \mathbf{x}(\tau\_1), \dots, \mathbf{x}(\tau\_N)]^T \in \prescript{(N+1)\times 2n}{}\\ \mathbf{U} &= [\mathbf{u}\_0, \mathbf{u}\_1, \dots, \mathbf{u}\_N]^T = [\mathbf{u}(\tau\_0), \mathbf{u}(\tau\_1), \dots, \mathbf{u}(\tau\_N)]^T \in \prescript{(N+1)\times m}{} \end{aligned} \tag{4.16}$$

Using (4.13) yields

$$\mathcal{D}\mathbf{X} = \dot{\mathbf{X}} \tag{4.17}$$

For the convenience of representation, rewrite the dynamic constraint (4.9) in problem *C* as

$$\dot{\mathbf{x}}(\tau) = \begin{bmatrix} \dot{\mathbf{x}}\_1(\tau) \\ \dot{\mathbf{x}}\_2(\tau) \end{bmatrix} = s \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{a}\_1 \ \mathbf{a}\_2 \end{bmatrix} \mathbf{x}(\tau) + s \begin{bmatrix} \mathbf{0} \\ \mathbf{g} \end{bmatrix} \mathbf{u}(\tau) \stackrel{\Delta}{=} \mathbf{F}(\mathbf{x}, \mathbf{u}) \tag{4.18}$$

Denote

$$\mathcal{F}(\mathbf{X}, \mathbf{U}) = \left[\mathbf{F}\_0, \mathbf{F}\_1, \dots, \mathbf{F}\_N\right]^T = \left[\mathbf{F}(\mathbf{x}(\tau\_0), \mathbf{u}(\tau\_0)), \mathbf{F}(\mathbf{x}(\tau\_1), \mathbf{u}(\tau\_1)), \dots, \mathbf{F}(\mathbf{x}(\tau\_N), \mathbf{u}(\tau\_N))\right]^T \tag{4.19}$$

Thus, problem *C* can be transcribed as

$$\begin{array}{ll}\text{Decision Variable}: \mathbf{X} = \in \prescript{(N+1)\times2n}{\end{array}}{\text{Problem } DPSC1: \min J^N[\mathbf{X}, \mathbf{U}]: = E(\mathbf{x}\_0, \mathbf{x}\_N) \\ \text{s.t.} \qquad \mathcal{D}\mathbf{X} = \mathcal{F}(\mathbf{X}, \mathbf{U}) \\ \mathbf{e}(\mathbf{x}\_0, \mathbf{x}\_N) = \mathbf{0} \end{array} \tag{4.20}$$

(2) An alternative PS method for problem *C* is taking its dynamic constraints as a second-order systems:

$$\ddot{\mathbf{x}}\_1(\tau) = s\dot{\mathbf{x}}\_2(\tau) = s^2 \left[ \mathbf{a}\_1 \mathbf{x}\_1(\tau) + \mathbf{a}\_2 \dot{\mathbf{x}}\_1(\tau) + \mathbf{g} \mathbf{u} \right] \stackrel{\Delta}{=} \mathbf{f}(\mathbf{x}\_1(\tau), \dot{\mathbf{x}}\_1(\tau), \mathbf{u}(\tau)) \quad (4.21)$$

and let

$$\mathbf{X}\_{1} = \begin{bmatrix} \mathbf{x}\_{10}, \mathbf{x}\_{11}, \dots, \mathbf{x}\_{1N} \end{bmatrix}^{T} = \begin{bmatrix} \mathbf{x}\_{1}(\tau\_{0}), \mathbf{x}\_{1}(\tau\_{1}), \dots, \mathbf{x}\_{1}(\tau\_{N}) \end{bmatrix}^{T} \in \prescript{(N+1)\times n}{} \tag{4.22}$$

$$\mathcal{F}(\mathbf{X}\_{\mathrm{l}}, \dot{\mathbf{X}}\_{\mathrm{l}}, \mathbf{U}) = [\mathbf{f}\_{0}, \dots, \mathbf{f}\_{N}]^{T} = \left[\mathbf{f}(\mathbf{x}\_{\mathrm{l}}(\tau\_{0}), \dot{\mathbf{x}}\_{\mathrm{l}}(\tau\_{0}), \mathbf{u}(\tau\_{0})) \dots, \mathbf{f}(\mathbf{x}\_{\mathrm{l}}(\tau\_{N}), \dot{\mathbf{x}}\_{\mathrm{l}}(\tau\_{N}), \mathbf{u}(\tau\_{N}))\right]^{T} \tag{4.23}$$

thus, the system (4.18) can be transcribed as

$$\mathcal{D}\mathbf{X}\_{1} = \dot{\mathbf{X}}\_{1} = \mathbf{s}\mathbf{X}\_{2} = \mathbf{s}[\mathbf{x}\_{20}, \mathbf{x}\_{21}, \dots, \mathbf{x}\_{2N}]^{T} = \mathbf{s}[\mathbf{x}\_{2}(\tau\_{0}), \mathbf{x}\_{2}(\tau\_{1}), \dots, \mathbf{x}\_{2}(\tau\_{N})]^{T} \tag{4.24}$$

and

$$\mathcal{D}^{(2)}\mathbf{X}\_{1} = \ddot{\mathbf{X}}\_{1} = \mathcal{F}(\mathbf{X}\_{1}, \dot{\mathbf{X}}\_{1}, \mathbf{U}) = \mathcal{F}(\mathbf{X}\_{1}, \mathbf{X}\_{2}, \mathbf{U}) \tag{4.25}$$

Combining (4.21), (4.25), problem *C* can be transcribed as the following secondorder form: