# Fault-Tolerant Guidance of Rocket Vertical Landing Phase Based on MPC Framework

Authors: Jingqi Li; Yaosong Long; Mao Su; Lei Liu; Bo Wang; Zhongtao Cheng
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

For the vertical landing process of reusable rockets, the landing accuracy is likely to be affected by disturbances and faults during flight. In this paper, a fault-tolerant guidance method based on the MPC framework is put forward. First, we propose a piecewise guidance algorithm that combines a trajectory optimization algorithm based on convex optimization with the MPC framework. With the fast trajectory optimization algorithm and the MPC framework that recursively introduces the real-time state, this algorithm forms a robust closed loop. Then, we design an integrated guidance, navigation, and control (GNC) system to enhance the fault tolerance and robustness of the guidance method. Simulation experiments verify that this method is fault-tolerant to various fault conditions including navigation system failures, control system failures, drag coefficient deviations, and atmospheric density deviations. This guidance method is robust enough to overcome disturbances and faults, and it has great potential for online use.


---

## Body

## 1. Introduction

Launching a rocket is a high-cost and high-risk investment, and there is growing concern about how to make it economical. Building reusable rockets is an effective way to reduce costs and create business value. At present, many countries have invested plenty of research in rocket recovery technology. Two US companies, Blue Origin and SpaceX, have successfully conducted rocket recovery experiments, proving the feasibility of vertical landing technology. The rocket recovery mission is generally divided into four flight phases: the attitude adjustment phase, the power deceleration phase, the aerodynamic deceleration phase, and the vertical landing phase [1], so as to ensure that the rocket lands vertically and softly at the predetermined location. The recovery process is confronted with difficulties including large spans of airspace and velocity domain, large changes in the flight environment, complex flight constraints, and strong disturbances and uncertainties, any of which may lead to errors. The accumulated errors generated by preceding phases must be eliminated in the vertical landing phase which is the last phase of the recovery mission [2]. The time is short but the precision is high. Therefore, extremely high requirements are placed on the guidance speed and accuracy of this phase. Research on fast and accurate guidance method of rocket vertical landing phase is the basis for successfully recovering a rocket.The study of the landing problem began with the Apollo project. Due to the limited computing ability, researchers could only apply analytical guidance law in that era [3]. Unlike the lunar landing mission, more complex path constraints and rigorous terminal conditions must be considered for rocket vertical landing problems. However, the derivation of the analytical method is cumbersome, making it unsuitable for most complex nonlinear problems [4]. Thanks to the research on numerical methods by mathematicians and the development of computers’ computing ability [2], the trajectory optimization of rocket vertical landing problems mainly relies on numerical algorithms these days.The vertical landing problem can be expressed as an optimal control problem with the objective of minimum fuel consumption, subject to the constraints of dynamic equations and the constraints on state and control variables. It can be transformed into a nonlinear programming problem (NLP) through time discretization. Early studies found that the pseudospectral methods had high accuracy in solving NLP. In recent years, these methods have been widely used for trajectory optimization problems. Chen and Xia [5] used a pseudospectral method to analyze the ascent trajectory characteristics of a solid rocket-powered launch vehicle. Mu et al. [6] planned the Mars landing trajectory by the Gauss pseudospectral method. However, the real-time performance and convergence of these methods cannot meet the requirements of online guidance [2]. Subsequent studies have shown that convex optimization is more advantageous for online trajectory optimization. As long as the problem is transformed into a convex optimization form, it can theoretically be solved in polynomial time by the interior point method [7]. Convex optimization was first applied to Mars landing [8, 9]. In order to convert the landing problem into a convex problem, Açıkmeşe et al. [10] proposed a lossless convexification method to deal with thrust constraints and proved that the lossless convexification problem is equivalent to the original problem, which greatly promoted the development of convex optimization methods in the field of trajectory optimization. However, not all constraints can be convexified by lossless convexification. A more general approach is successive convexification, which linearizes the original nonconvex problem into a series of convex subproblems and then iteratively finds the optimal solution of the subproblems. In 2013, Lu and Liu [11] proposed successive convexification for convexifing dynamic equations and other nonconvex constraints when they studied the rendezvous problem. They later used this method to solve the hypersonic vehicle reentry problem [12] and the rocket landing problem [13]. The advantages of convex optimization including computational efficiency and insensitivity to initial guesses make it an efficient method for solving trajectory optimization problems. It is widely used in the aerospace field, such as low-thrust transfers [14], rendezvous problem [15], re-entry problem [12], and the vertical landing problem studied in this paper.There are two commonly used coordinate systems for modeling rocket vertical landing problems, the velocity coordinate system [16–20] and the landing point coordinate system [21–25]. In this paper, the three-degree-of-freedom dynamics of the rocket are derived based on the landing point coordinate system. Since Szmuk et al. [21] applied successive convex optimization to solve the rocket vertical landing problem in 2016, research on trajectory optimization of landing problems has become more and more abundant and mature. The current research focus is to improve the convergence performance and solution efficiency of successive convex optimization, so that the approach can meet the needs of online guidance. In Ref. [24], a convex feasible set (CFS) method is proposed to convexify the angle of attack constraint which is a nonconvex-nonconcave inequality. Ref. [25] presents a two-stage successive convexification method. Simulation experiments show that the performance of the two-stage method is more stable and efficient than the single-stage method. Ref. [26] compares the computational performance and solution accuracy of six discretization methods. In Ref. [18], an online update strategy for trust regions is proposed to speed up the convergence of successive convex optimization.However, the rocket will be affected by wind interference and environmental uncertainties after entering the atmosphere, and modules such as the navigation system and the control system may malfunction during the landing process. Trajectory optimization alone cannot overcome these disturbances and faults. Therefore, it is necessary to design an online fault-tolerant guidance method to ensure landing accuracy. A fault-tolerant control method has been proposed to deal with parametric uncertainties and unknown actuator failures [27]. But there is a lack of research on guidance methods for vertical landing problems at present. In Ref. [22], a receding horizon guidance method based on convex optimization is proposed. Ref. [28] and Ref. [29] both propose to construct a guidance, navigation, and control (GNC) system for closed-loop guidance, but they do not give a specific algorithm.In recent years, some researchers put forward a guidance method based on trajectory optimization and model predictive control (MPC) framework. MPC is a control strategy that recursively solves an optimal control problem with updated system states at each sampling time. Ref. [30] presents the MPC algorithm for the optimal guidance and reconfiguration of swarms of spacecraft. In Ref. [31], MPC is used for asteroid landing. Ref. [32] reviews applications of MPC in the aerospace guidance field. For the online guidance problem of the rocket vertical landing phase, a successive convexification + MPC guidance algorithm is proposed by Ref. [33]. Ref. [19] embeds a pseudospectral-improved successive convexification (PISC) algorithm in the MPC framework to construct a parallel feasibility-guaranteed guidance algorithm. Ref. [34] designs an antidelay model predictive control (AD-MPC) scheme for carrier landing. Ref. [35] implements a successive convexification MPC-based guidance algorithm to solve the six-degree-of-freedom powered descent guidance problem.In this paper, we propose a piecewise guidance algorithm that embeds a convex optimization-based trajectory optimization algorithm in the MPC framework. An integrated GNC system is then designed to further improve the fault tolerance and robustness of the entire system. The rest of this paper is organized as follows. In Section2, a mathematical description of the rocket vertical landing problem is given. Section 3 elaborates the online trajectory optimization algorithm based on convex optimization, including convexification and discretization methods for transforming the original problem into a second-order cone programming (SOCP) problem. In Section 4, the trajectory optimization algorithm is embedded in the MPC framework, and a piecewise guidance method is proposed to cope with the bang-bang characteristic of the control. Section 5 designs an integrated GNC system. In Section 6, simulation experiments are carried out to verify the reliability, accuracy, robustness, and fault tolerance of the guidance method. Section 7 concludes the whole work.


## 2. Problem Formulation

In this section, we use the inertial reference frame to formulate the fuel-optimal rocket vertical landing problem as a standard optimal control problem.


### 2.1. Dynamics and Constraints

The rocket’s flight time during the vertical landing phase is so short that the Earth’s surface can be assumed to be a horizontal plane. Based on this assumption, we build an inertial reference frame with the origin O located at the predetermined landing point as shown in Figure1. The OX axis is perpendicular to the plane, and the upward direction is positive. The OY axis is parallel to the plane, and the direction pointing to the rocket’s launch point is positive. The OZ axis and the other two axes form a right-handed Cartesian coordinate system.Figure 1
Inertial reference frame.In this reference frame, the rocket’s three-degree-of-freedom dynamics are(1)ṙt=vt,v̇t=Tt+Dtmt+g,ṁt=−TtIspg0,where rt, vt, and mt represent the position vector, velocity vector, and mass of the rocket, respectively. Tt is the thrust vector of the rocket engine. · represents the 2-norm of the vector.During the landing phase, the constraints on thrust magnitude and direction are expressed as(2)Tmin≤Tt≤Tmax,(3)Ttcosηmax≤exTt,where Tmin and Tmax represent the minimum and maximum thrust magnitude that the engine can provide, respectively. ηmax is the maximum allowable value of the angle between the thrust direction and the OX direction. Dt is the aerodynamic drag calculated by the following formula:
(4)Dt=−12ρSDCDvtvt,where ρ is the air density, SD is the drag reference area, and CD is the drag coefficient.To prevent the rocket from colliding with the ground or being interfered with by the ground protrusion during flight, a glide-slope constraint is imposed to restrict the path of the rocket to lie within an upward-facing cone:(5)rtcosθmax≤exrt,where ex represents the unit vector in the OX direction. θmax is the maximum allowable half-cone angle.We also need to consider boundary conditions. The rocket’s initial position vector, velocity vector, and mass are specified as fixed parameters. Its final position and velocity vectors are fixed, too. The thrust vector at the terminal time must be in the OX direction. The fuel remaining of the rocket must be nonnegative, which means that the landing mass of the rocket must be greater than the dry mass which is denoted asmdry. The boundary conditions are summarized as
(6)r0=r0,v0=v0,m0=m0,(7)rtf=0,vtf=0,Ttf=Ttfex,mtf≥mdry.


### 2.2. Performance Index

The performance index of the fuel-optimal rocket vertical landing problem is selected as minimizing the fuel consumption, which is equivalent to maximizing the terminal mass of the rocket. Therefore, the objective function is expressed as(8)J=−mtf.To sum up, with thrust vectorT as the control variable and rT,vT,mT as the state variables, the rocket vertical landing problem can be formulated as an optimal control problem with free terminal time:
(9)Problem0:minTJ=−mtf,subjectto123567.


## 2.1. Dynamics and Constraints

The rocket’s flight time during the vertical landing phase is so short that the Earth’s surface can be assumed to be a horizontal plane. Based on this assumption, we build an inertial reference frame with the origin O located at the predetermined landing point as shown in Figure1. The OX axis is perpendicular to the plane, and the upward direction is positive. The OY axis is parallel to the plane, and the direction pointing to the rocket’s launch point is positive. The OZ axis and the other two axes form a right-handed Cartesian coordinate system.Figure 1
Inertial reference frame.In this reference frame, the rocket’s three-degree-of-freedom dynamics are(1)ṙt=vt,v̇t=Tt+Dtmt+g,ṁt=−TtIspg0,where rt, vt, and mt represent the position vector, velocity vector, and mass of the rocket, respectively. Tt is the thrust vector of the rocket engine. · represents the 2-norm of the vector.During the landing phase, the constraints on thrust magnitude and direction are expressed as(2)Tmin≤Tt≤Tmax,(3)Ttcosηmax≤exTt,where Tmin and Tmax represent the minimum and maximum thrust magnitude that the engine can provide, respectively. ηmax is the maximum allowable value of the angle between the thrust direction and the OX direction. Dt is the aerodynamic drag calculated by the following formula:
(4)Dt=−12ρSDCDvtvt,where ρ is the air density, SD is the drag reference area, and CD is the drag coefficient.To prevent the rocket from colliding with the ground or being interfered with by the ground protrusion during flight, a glide-slope constraint is imposed to restrict the path of the rocket to lie within an upward-facing cone:(5)rtcosθmax≤exrt,where ex represents the unit vector in the OX direction. θmax is the maximum allowable half-cone angle.We also need to consider boundary conditions. The rocket’s initial position vector, velocity vector, and mass are specified as fixed parameters. Its final position and velocity vectors are fixed, too. The thrust vector at the terminal time must be in the OX direction. The fuel remaining of the rocket must be nonnegative, which means that the landing mass of the rocket must be greater than the dry mass which is denoted asmdry. The boundary conditions are summarized as
(6)r0=r0,v0=v0,m0=m0,(7)rtf=0,vtf=0,Ttf=Ttfex,mtf≥mdry.


## 2.2. Performance Index

The performance index of the fuel-optimal rocket vertical landing problem is selected as minimizing the fuel consumption, which is equivalent to maximizing the terminal mass of the rocket. Therefore, the objective function is expressed as(8)J=−mtf.To sum up, with thrust vectorT as the control variable and rT,vT,mT as the state variables, the rocket vertical landing problem can be formulated as an optimal control problem with free terminal time:
(9)Problem0:minTJ=−mtf,subjectto123567.


## 3. Trajectory Optimization Algorithm Based on Convex Optimization

In this section, we will elaborate on the trajectory optimization algorithm based on convex optimization and embedded in the subsequent guidance system. This includes the use of convexification and discretization methods to convert Problem0 into a convex problem, as well as iteratively solving the converted problem afterward.


### 3.1. Convexification

The thrust magnitude constraint in Problem0 is nonconvex. Its nonconvex feasible region can be relaxed through lossless convexification to form a high-dimensional convex feasible region, transforming the nonconvex constraint into a relaxed convex constraint. Using the above method, we introduce a slack variableΓ to transform Eq. (2) into
(10)Tt≤Γt,(11)Tmin≤Γt≤Tmax.By replacing the nonconvex constraint (Eq. (2)) in Problem0 with the convex constraints (Eqs. (10) and (11)), the nonconvex problem is converted into a relaxed problem. Studies have shown that the optimal solution to the relaxed problem is also the optimal solution to Problem0 [8–10]. Taking T and Γ as control variables, that is, the control variables are redefined as u=Γ,TTT. The dynamics are transformed into
(12)fx,u=ṙt=vtv̇t=Tt+Dtmt+gṁt=−ΓIspg0.The nonlinearity of Eq. (12) is caused by the free terminal time, aerodynamic drag Dt, and the denominator mt. It can be linearized through successive convexification. In this paper, the terminal time is added to the optimization variables, which will be optimized together with the state and control variables in the subsequent optimization process. tf in Eq. (12) is a hidden variable. Define τ≜t/tf. Apply the chain rule to make tf explicit:
(13)dxdt=dxdτdτdt.Eq. (12) is transformed into
(14)dxdτ=fx,utf.We adopt successive convexification to iteratively solve Eq. (14). The first-order Taylor expansion is performed at the kth iteration to transform Eq. (14) into
(15)dxdτ=fxk,uktfk+Axk,ukx−xktfk+Bxk,uku−uktfk+fxk,uktf−tfk,where Axk,uk and Bxk,uk are the gradients of fx,u with respect to the state variables x=rT,vT,mT and control variables u=Γ,TTT, respectively.


### 3.2. Discretization

The problem is still infinite-dimensional after convexification and needs to be discretized into a finite-dimensional problem. We adopt the trapezoidal discretization to discretize the problem. This requires the time of flight to be evenly divided intoN discrete intervals, i.e., N+1 discrete points (the value of N in each guidance cycle will be discussed later). The subscript ii=1,⋯,N+1 represents the ith discrete point of the state or control variable, then, Eq. (15) is discretized as
(16)xi+1=xi+12NAikxi−xiktfk+Bikui−uiktfk+fiktf+12NAi+1kxi+1−xi+1ktfk+Bi+1kui+1−ui+1ktfk+fi+1ktf,where Aik=Axik,uik, Bik=Bxik,uik, and fik=fxik,uik.After discretization, constraints (3), (5), (6), (7), (10), and (11) are transformed into
(17)Ticosηmax≤exTi,(18)ricosθmax≤exri,(19)r1=r0,v1=v0,m1=m0,(20)rN+1=0,vN+1=0,TN+1=TN+1ex,mN+1≥mdry,(21)Ti≤Γi,(22)Tmin≤Γi≤Tmax.


### 3.3. Trajectory Optimization Algorithm

DefineHi=I+tfk/2NAik, Hi+1=−I+tfk/2NAi+1k, Gi=tfk/2NBik, Gi+1=tfk/2NBi+1k, Fi=fik+fi+1k/2N, Ci=Aikxiktfk+Bikuiktfk. Equation (16) can be expressed as
(23)Hixi+Hi+1xi+1+Giui+Gi+1ui+1+Fi=12NCi+Ci+1.The terminal timetf is added to the control variables for optimization. Combine the discretized forms of state variables x=rT,vT,mT, control variables u=Γ,TTT and tf into a joint optimization variable Z, that is, define Z=x1T,⋯,xN+1T,u1T,⋯,uN+1T,tfT, then, Eq. (23) can be transformed into the following matrix form:
(24)MZ=C,whereM=I0⋯0H1H2⋯0⋮⋮⋱⋮00⋯HN00⋮HN+100⋯0G1G2⋯0⋮⋮⋱⋮00⋯GN00⋮GN+10F1⋮FN,C=1/2N2Nx1C1+C2⋮CN+CN+1.Problem0 is transformed into(25)Problem1:minΓ,T,tfJ=−mN+1,subjectto24171819202122.Problem1 is an SOCP problem that can be solved iteratively using the interior-point algorithm. The solution steps are shown in Figure2.Figure 2
Solution steps for the trajectory optimization problem.The convergence condition in Figure2 is
(26)maxixik+1−xik≤εx,maxiuik+1−uik≤εu,which means that if the difference between the values of the optimization variables of two consecutive iterations is within the specified tolerance range, the iteration should be stopped and the optimal solution is obtained.


## 3.1. Convexification

The thrust magnitude constraint in Problem0 is nonconvex. Its nonconvex feasible region can be relaxed through lossless convexification to form a high-dimensional convex feasible region, transforming the nonconvex constraint into a relaxed convex constraint. Using the above method, we introduce a slack variableΓ to transform Eq. (2) into
(10)Tt≤Γt,(11)Tmin≤Γt≤Tmax.By replacing the nonconvex constraint (Eq. (2)) in Problem0 with the convex constraints (Eqs. (10) and (11)), the nonconvex problem is converted into a relaxed problem. Studies have shown that the optimal solution to the relaxed problem is also the optimal solution to Problem0 [8–10]. Taking T and Γ as control variables, that is, the control variables are redefined as u=Γ,TTT. The dynamics are transformed into
(12)fx,u=ṙt=vtv̇t=Tt+Dtmt+gṁt=−ΓIspg0.The nonlinearity of Eq. (12) is caused by the free terminal time, aerodynamic drag Dt, and the denominator mt. It can be linearized through successive convexification. In this paper, the terminal time is added to the optimization variables, which will be optimized together with the state and control variables in the subsequent optimization process. tf in Eq. (12) is a hidden variable. Define τ≜t/tf. Apply the chain rule to make tf explicit:
(13)dxdt=dxdτdτdt.Eq. (12) is transformed into
(14)dxdτ=fx,utf.We adopt successive convexification to iteratively solve Eq. (14). The first-order Taylor expansion is performed at the kth iteration to transform Eq. (14) into
(15)dxdτ=fxk,uktfk+Axk,ukx−xktfk+Bxk,uku−uktfk+fxk,uktf−tfk,where Axk,uk and Bxk,uk are the gradients of fx,u with respect to the state variables x=rT,vT,mT and control variables u=Γ,TTT, respectively.


## 3.2. Discretization

The problem is still infinite-dimensional after convexification and needs to be discretized into a finite-dimensional problem. We adopt the trapezoidal discretization to discretize the problem. This requires the time of flight to be evenly divided intoN discrete intervals, i.e., N+1 discrete points (the value of N in each guidance cycle will be discussed later). The subscript ii=1,⋯,N+1 represents the ith discrete point of the state or control variable, then, Eq. (15) is discretized as
(16)xi+1=xi+12NAikxi−xiktfk+Bikui−uiktfk+fiktf+12NAi+1kxi+1−xi+1ktfk+Bi+1kui+1−ui+1ktfk+fi+1ktf,where Aik=Axik,uik, Bik=Bxik,uik, and fik=fxik,uik.After discretization, constraints (3), (5), (6), (7), (10), and (11) are transformed into
(17)Ticosηmax≤exTi,(18)ricosθmax≤exri,(19)r1=r0,v1=v0,m1=m0,(20)rN+1=0,vN+1=0,TN+1=TN+1ex,mN+1≥mdry,(21)Ti≤Γi,(22)Tmin≤Γi≤Tmax.


## 3.3. Trajectory Optimization Algorithm

DefineHi=I+tfk/2NAik, Hi+1=−I+tfk/2NAi+1k, Gi=tfk/2NBik, Gi+1=tfk/2NBi+1k, Fi=fik+fi+1k/2N, Ci=Aikxiktfk+Bikuiktfk. Equation (16) can be expressed as
(23)Hixi+Hi+1xi+1+Giui+Gi+1ui+1+Fi=12NCi+Ci+1.The terminal timetf is added to the control variables for optimization. Combine the discretized forms of state variables x=rT,vT,mT, control variables u=Γ,TTT and tf into a joint optimization variable Z, that is, define Z=x1T,⋯,xN+1T,u1T,⋯,uN+1T,tfT, then, Eq. (23) can be transformed into the following matrix form:
(24)MZ=C,whereM=I0⋯0H1H2⋯0⋮⋮⋱⋮00⋯HN00⋮HN+100⋯0G1G2⋯0⋮⋮⋱⋮00⋯GN00⋮GN+10F1⋮FN,C=1/2N2Nx1C1+C2⋮CN+CN+1.Problem0 is transformed into(25)Problem1:minΓ,T,tfJ=−mN+1,subjectto24171819202122.Problem1 is an SOCP problem that can be solved iteratively using the interior-point algorithm. The solution steps are shown in Figure2.Figure 2
Solution steps for the trajectory optimization problem.The convergence condition in Figure2 is
(26)maxixik+1−xik≤εx,maxiuik+1−uik≤εu,which means that if the difference between the values of the optimization variables of two consecutive iterations is within the specified tolerance range, the iteration should be stopped and the optimal solution is obtained.


## 4. Piecewise Guidance Algorithm Based on MPC Framework

In this section, on the basis of the trajectory optimization algorithm based on convex optimization detailed in the previous section, a piecewise guidance method combining the trajectory optimization algorithm with the MPC framework is proposed. At each sampling time, the real-time state of the rocket fed back by the navigation system is employed as the initial state to start the trajectory optimization algorithm. Then, the optimal values of a series of control variables from the current time to the landing point are obtained. But only the values in the current guidance cycle are applied to control the rocket until the next sampling time. Repeat these steps until the rocket lands at the predetermined point. Since the trajectory optimization algorithm needs to be provided with an initial guess, the optimization result obtained at the previous sampling time can be used as the initial trajectory at the current sampling time. The computational efficiency of the trajectory optimization algorithm based on convex optimization makes it possible for the guidance algorithm to form an effective closed loop. In addition, the practice of recursively introducing the real-time state of the rocket into the trajectory optimization ensures the stability, robustness, and fault tolerance of the closed-loop.The piecewise method is to cope with the thrust magnitude’s bang-bang characteristic. Bang-bang control means the control amount changes dramatically in a short period near the switching point. If the duration of each guidance cycle is too long and a fault occurs near the switching point, the trajectory will deviate from the optimal solution.The piecewise guidance algorithm based on MPC framework is shown in algorithm1, where ti represents the sampling time of the ith trajectory optimization, Ti represents the terminal time obtained by the ith trajectory optimization, T represents the thrust series from the current time to the landing point, while Ti contains only the portion within the ith guidance cycle, Δt is the duration of each guidance cycle which is set initially, and Δti is the actual duration of the ith guidance cycle, Δti=ti+1−ti.Algorithm 1: Piecewise guidance algorithm based on MPC framework.
Piecewise guidance algorithm based on MPC frameworkInitialization: Save the initial trajectory to the database, set relevant parameters required by the trajectory optimization algorithm, and set the update index i=1.1whileTi−1−ti−1>02  generate guidance commands according to the guidance cycle clock;3  sample the current statexti(i.e. rti,vti,mti) of the rocket;4  employ the optimal trajectoryJi−1 which was obtained in the previous guidance cycle as the initial guess. Start the trajectory optimization using the current state xti. Obtain a new optimal trajectory Ji and an optimal thrust series T. Save Ji to the database, implement Ti, and set i=i+1;5if a fault occurs and ti approaches the switching point6   shorten the duration of the guidance cycle:Δti=α⋅Δt0<α<1 ;7else8Δti=Δt ;9end if10end whileOver time, the time horizon of the optimization becomes shorter, and the number of discrete points for trajectory optimization within the guidance cycle should also be reduced accordingly. The number of discrete points taken by theith trajectory optimization algorithm is determined by
(27)Ni=Ni−1Ti−1−tiTi−1−ti−1,where the ceiling function is defined as x=minn∈ℤn≥x. Ti−1−ti−1is the time horizon of the i−1th optimization, and Ti−1−ti is the estimate of the time horizon of the ith optimization.


## 5. Integrated GNC System Design

In actual flight, a key technology to realize fault-tolerant guidance is designing a guidance, navigation, and control (GNC) system. The whole GNC system is closed-loop, where the navigation system consists of a variety of sensors for real-time measurement and evaluation of states including the rocket’s position, velocity, and mass, the guidance system consists of the guidance algorithm proposed in the previous section to solve the landing trajectory, and the control system consists of the actuators to control the rocket to track the trajectory. The guidance system needs to provide the real-time state of the rocket to the guidance system before the trajectory optimization. And the control commands derived from the guidance algorithm need to be executed by the control system. Thanks to the division of labor and rapid cooperation of each system, the integrated GNC system has strong fault tolerance and robustness against various disturbances and faults including navigation system failures, control system failures, drag coefficient deviations, and atmospheric density deviations.In this section, an integrated GNC system is built, as shown in Figure3. In a guidance cycle, based on the current state rti,vti,mti fed back by the navigation system, the guidance system adopts the trajectory optimization algorithm based on convex optimization to generate an optimal thrust series T (where Ji−1 represents the initial guess which is the optimal trajectory of the last guidance). The control system outputs the thrust command Ti of the current guidance cycle and transmits it to the rocket, and the rocket system executes the command accordingly. In the simulation experiments, we use numerical integration to integrate the state of the rocket after one guidance cycle, which is used to simulate the real-time state of the rocket. Yet in actual flight, the state is measured by the navigation system and transmitted to the guidance system. Here comes the next guidance cycle, repeat the above until the rocket lands. The piecewise guidance algorithm proposed in the previous section and the integrated GNC system designed in this section ensure the fault tolerance of the guidance method.Figure 3
Integrated GNC system.


## 6. Simulation Experiments

In this section, numerical simulations are provided to verify the reliability, accuracy, fault tolerance, and robustness of the proposed guidance method. All simulation experiments are carried out on MATLAB with the use of the modeling tool CVX [36, 37] to establish the guidance problem and the solver MOSEK [38] to solve it. We employ the fourth-order Runge-Kutta integration with a 0.01 s time step to obtain updated states at each sampling time. The parameters used in the simulations are shown in Table 1.Table 1
Parameter values.
ParameterValueUnitsr03500,700,0Tmv0−200,−90,0Tm/sm027000kgmdry20000kgTmin100kNTmax300kNθmax80°ηmax15°g09.8m/s2Isp270sρ1.225kg/m3SD10m2CD2.2—Δt1sThe number of discrete points for the first trajectory optimization is:N=40. Assuming ri0, vi0, and mi0 vary linearly from the initial value to the terminal value, then, the initial guesses for all optimization parameters are
(28)ri0=rf−r0i−1N+r0,vi0=vf−v0i−1N+v0,mi0=mdry−m0i−1N+m0,Γi0=Tmin,Ti0=−Γi0vi0vi0,tf0=40s.The convergence condition is(29)εr=10−3m,εv=10−3m/s,εm=10−3kg,εΓ=10−3kN.


### 6.1. Reliability and Accuracy

In this subsection, we conduct the simulations in the absence of any faults and deviations. Figures4 and 5 show the trajectories and control commands of the rocket vertical landing phase obtained by the trajectory optimization algorithm and the guidance method presented in this paper, respectively. The red curves represent the results of the trajectory optimization algorithm, and the black is the guidance method. The differences between the results are small, demonstrating the reliability of the guidance method. We then perform numerical integration using the control commands to obtain the actual terminal states. Tables 2 and 3 show the terminal states and landing errors, respectively. It can be seen that the landing position and velocity errors of the guidance method are smaller than those of the trajectory optimization algorithm, which proves that the guidance method proposed in this paper has higher accuracy.Figure 4
The trajectories of rocket vertical landing phase obtained by the guidance method and the trajectory optimization algorithm.Figure 5
The control commands of rocket vertical landing phase obtained by the guidance method and the trajectory optimization algorithm.Table 2
Terminal states.
tfsrxtfmrytfmvxtfm/svytfm/smtfkgTrajectory optimization43.18-2.0085-1.1821-0.0646-0.033223137.09Guidance43.18-0.0276-0.0241-0.0210-0.004323135.79Table 3
Landing errors.
Landing position error/mLanding velocity error/(m/s)Trajectory optimization2.33050.0726Guidance0.03660.0214


### 6.2. Fault Tolerance and Robustness

In this subsection, we assume three scenarios: faults occur in the navigation system, faults occur in the control system, and estimates of parameters (the drag coefficient and the atmospheric density) have deviations. Considering these scenarios, we apply the trajectory optimization algorithm and the guidance method proposed in this paper to conduct simulation experiments, respectively, and compare their results.


#### 6.2.1. Faults Occur in the Navigation System

It is assumed that the navigation system has faults. Rather, the measurement of the real-time state of the rocket has the following errors:(30)rmeasure=rtrue±min10,rtrue×ε%×random0,1,vmeasure=vtrue±min1,vtrue×ε%×random0,1,where 10 and 1 represent the absolute errors of the rocket’s position and velocity, ε represents the maximum range of the relative error, random0,1 represents a random number generated from [0,1], rtrue,vtrue represents the actual state of the rocket, and rmeasure,vmeasure represents the state measured by the navigation system and fed back to the guidance system. We take the value of ε as 5, 10, and 15 for experiments, and the landing errors under these fault conditions are shown in Table 4.Table 4
Landing errors and masses under the condition that the navigation system has faults.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgε=52.46655.777923331.19ε=103.57505.830123329.48ε=154.06246.551423355.94It can be seen from Table4 that as the relative error of the measured state increases, the landing errors of the rocket also increase, but they remain within a small range. This experiment proves that the guidance method proposed in this paper is fault-tolerant and robust to navigation system faults.


#### 6.2.2. Faults Occur in the Control System

It is assumed that the control system has two fault conditions: thrust cannot change continuously (that is, thrust is constant during each guidance cycle), and thrust magnitude has deviations. The trajectory optimization algorithm and the guidance method proposed in this paper are, respectively, applied to conduct simulations under the above two fault conditions, and their results are compared.(1) Thrust Cannot Change Continuously. For trajectory optimization, we assume that thrust is constant per second, and for guidance, thrust is constant during each guidance cycle. The landing errors and landing masses are shown in Table 5.Table 5
Landing errors and masses under the condition that thrust cannot change continuously.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgTrajectory optimization407.64014.806223081.57Guidance0.09140.314523136.54(2) Thrust Magnitude Has Deviations. It is assumed that there are deviations within ±5% of the thrust magnitude. In the experiment, random deviations within ±5% are applied to the thrust magnitude obtained by the trajectory optimization algorithm and the guidance method. Then, the deviated thrust is used to control the rocket. The landing errors and landing masses are shown in Table 6.Table 6
Landing errors and masses under the condition that thrust magnitude has deviations.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgTrajectory optimization359.13676.393923137.13Guidance0.01850.020923133.57It can be seen from Tables5 and 6 that when the thrust cannot be continuously changed or the thrust magnitude has deviations, the landing position error and velocity error obtained by trajectory optimization are significantly larger than those obtained by guidance. The landing position errors of the guidance method under the two fault conditions remain at the same order of magnitude as that under the no-fault condition in Subsection 6.1, which demonstrates that the two fault conditions have little effect on the landing position. As for the landing velocity error, the guidance method can make it at the same order of magnitude as the no-fault condition when the thrust magnitude has deviations within ±5%. The error is one order of magnitude greater when the thrust cannot be continuously changed, but it still does not exceed 1 m/s. The simulation experiments demonstrate that the guidance method proposed in this paper is fault-tolerant and robust to control system faults.


#### 6.2.3. Estimates of Parameters Have Deviations

It is assumed that the actual drag coefficient and atmospheric density differ by±15% from the estimated values. We apply the trajectory optimization algorithm and the guidance method proposed in this paper to conduct simulation experiments using the estimated parameters. Then, we perform numerical integration using the actual parameters to obtain real terminal states. The landing errors and landing masses of the two algorithms under the fault condition are shown in Tables 7 and 8, respectively, (where CD+15% indicates that the actual drag coefficient has a value of 1.15×CD).Table 7
Landing errors and masses of the trajectory optimization algorithm.
DeviationLanding position error/mLanding velocity error/(m/s)Landing mass/kgCD+15%289.59334.979023135.95CD−15%356.56896.259823135.95ρ+15%289.95844.979023135.95ρ−15%356.56896.259823135.95Table 8
Landing errors and masses of the guidance method.
DeviationLanding position error/mLanding velocity error/(m/s)Landing mass/kgCD+15%0.00350.023923230.63CD−15%4.162729.245623957.32ρ+15%0.00350.023923230.63ρ−15%4.162729.245623957.32It can be seen from Tables7 and 8 that the guidance method proposed in this paper can significantly improve the landing position accuracy in the case of the drag coefficient or the atmospheric density having deviations. When the actual values of the two parameters are larger than the estimated values, the accuracy of the landing velocity is also significantly improved. This experiment demonstrates the fault tolerance and robustness of the guidance method to the landing position in the presence of parameter deviations. However, when the actual values of the two parameters are smaller than the estimated values, the landing velocity errors obtained by the guidance method will become larger. We refer to the analysis of Ref. [19] and learn that when the actual values of the two parameters are decreased, part of the mechanical energy of the rocket cannot be dissipated as expected; and when the thrust magnitude saturates, there is no additional energy to compensate for the undissipated mechanical energy, leading to the error of the optimal trajectory becomes larger. Subsequent research on guidance methods should attempt to address this issue.


#### 6.2.4. Combination of Various Faults

Finally, we consider navigation system failures, control system failures (thrust magnitude has deviations), drag coefficient deviations, and atmospheric density deviations at the same time and carry out Monte Carlo simulations (200 cases). For the navigation system and control system, the failures are added in the same way as the previous experiments. Drag coefficient deviations and atmospheric density deviations are considered to be normally distributed, and their means are set to zero. The3σ values are, respectively, taken as CD~15%, ρ~15%. Results of terminal position and terminal velocity are shown in Figure 6. The maximum errors are (7.2551 m, 31.6108 m/s). This experiment proves that the guidance method can still keep the terminal states within acceptable limits under the combination of various faults.Figure 6
Monte Carlo results of terminal position and terminal velocity.


## 6.1. Reliability and Accuracy

In this subsection, we conduct the simulations in the absence of any faults and deviations. Figures4 and 5 show the trajectories and control commands of the rocket vertical landing phase obtained by the trajectory optimization algorithm and the guidance method presented in this paper, respectively. The red curves represent the results of the trajectory optimization algorithm, and the black is the guidance method. The differences between the results are small, demonstrating the reliability of the guidance method. We then perform numerical integration using the control commands to obtain the actual terminal states. Tables 2 and 3 show the terminal states and landing errors, respectively. It can be seen that the landing position and velocity errors of the guidance method are smaller than those of the trajectory optimization algorithm, which proves that the guidance method proposed in this paper has higher accuracy.Figure 4
The trajectories of rocket vertical landing phase obtained by the guidance method and the trajectory optimization algorithm.Figure 5
The control commands of rocket vertical landing phase obtained by the guidance method and the trajectory optimization algorithm.Table 2
Terminal states.
tfsrxtfmrytfmvxtfm/svytfm/smtfkgTrajectory optimization43.18-2.0085-1.1821-0.0646-0.033223137.09Guidance43.18-0.0276-0.0241-0.0210-0.004323135.79Table 3
Landing errors.
Landing position error/mLanding velocity error/(m/s)Trajectory optimization2.33050.0726Guidance0.03660.0214


## 6.2. Fault Tolerance and Robustness

In this subsection, we assume three scenarios: faults occur in the navigation system, faults occur in the control system, and estimates of parameters (the drag coefficient and the atmospheric density) have deviations. Considering these scenarios, we apply the trajectory optimization algorithm and the guidance method proposed in this paper to conduct simulation experiments, respectively, and compare their results.


### 6.2.1. Faults Occur in the Navigation System

It is assumed that the navigation system has faults. Rather, the measurement of the real-time state of the rocket has the following errors:(30)rmeasure=rtrue±min10,rtrue×ε%×random0,1,vmeasure=vtrue±min1,vtrue×ε%×random0,1,where 10 and 1 represent the absolute errors of the rocket’s position and velocity, ε represents the maximum range of the relative error, random0,1 represents a random number generated from [0,1], rtrue,vtrue represents the actual state of the rocket, and rmeasure,vmeasure represents the state measured by the navigation system and fed back to the guidance system. We take the value of ε as 5, 10, and 15 for experiments, and the landing errors under these fault conditions are shown in Table 4.Table 4
Landing errors and masses under the condition that the navigation system has faults.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgε=52.46655.777923331.19ε=103.57505.830123329.48ε=154.06246.551423355.94It can be seen from Table4 that as the relative error of the measured state increases, the landing errors of the rocket also increase, but they remain within a small range. This experiment proves that the guidance method proposed in this paper is fault-tolerant and robust to navigation system faults.


### 6.2.2. Faults Occur in the Control System

It is assumed that the control system has two fault conditions: thrust cannot change continuously (that is, thrust is constant during each guidance cycle), and thrust magnitude has deviations. The trajectory optimization algorithm and the guidance method proposed in this paper are, respectively, applied to conduct simulations under the above two fault conditions, and their results are compared.(1) Thrust Cannot Change Continuously. For trajectory optimization, we assume that thrust is constant per second, and for guidance, thrust is constant during each guidance cycle. The landing errors and landing masses are shown in Table 5.Table 5
Landing errors and masses under the condition that thrust cannot change continuously.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgTrajectory optimization407.64014.806223081.57Guidance0.09140.314523136.54(2) Thrust Magnitude Has Deviations. It is assumed that there are deviations within ±5% of the thrust magnitude. In the experiment, random deviations within ±5% are applied to the thrust magnitude obtained by the trajectory optimization algorithm and the guidance method. Then, the deviated thrust is used to control the rocket. The landing errors and landing masses are shown in Table 6.Table 6
Landing errors and masses under the condition that thrust magnitude has deviations.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgTrajectory optimization359.13676.393923137.13Guidance0.01850.020923133.57It can be seen from Tables5 and 6 that when the thrust cannot be continuously changed or the thrust magnitude has deviations, the landing position error and velocity error obtained by trajectory optimization are significantly larger than those obtained by guidance. The landing position errors of the guidance method under the two fault conditions remain at the same order of magnitude as that under the no-fault condition in Subsection 6.1, which demonstrates that the two fault conditions have little effect on the landing position. As for the landing velocity error, the guidance method can make it at the same order of magnitude as the no-fault condition when the thrust magnitude has deviations within ±5%. The error is one order of magnitude greater when the thrust cannot be continuously changed, but it still does not exceed 1 m/s. The simulation experiments demonstrate that the guidance method proposed in this paper is fault-tolerant and robust to control system faults.


### 6.2.3. Estimates of Parameters Have Deviations

It is assumed that the actual drag coefficient and atmospheric density differ by±15% from the estimated values. We apply the trajectory optimization algorithm and the guidance method proposed in this paper to conduct simulation experiments using the estimated parameters. Then, we perform numerical integration using the actual parameters to obtain real terminal states. The landing errors and landing masses of the two algorithms under the fault condition are shown in Tables 7 and 8, respectively, (where CD+15% indicates that the actual drag coefficient has a value of 1.15×CD).Table 7
Landing errors and masses of the trajectory optimization algorithm.
DeviationLanding position error/mLanding velocity error/(m/s)Landing mass/kgCD+15%289.59334.979023135.95CD−15%356.56896.259823135.95ρ+15%289.95844.979023135.95ρ−15%356.56896.259823135.95Table 8
Landing errors and masses of the guidance method.
DeviationLanding position error/mLanding velocity error/(m/s)Landing mass/kgCD+15%0.00350.023923230.63CD−15%4.162729.245623957.32ρ+15%0.00350.023923230.63ρ−15%4.162729.245623957.32It can be seen from Tables7 and 8 that the guidance method proposed in this paper can significantly improve the landing position accuracy in the case of the drag coefficient or the atmospheric density having deviations. When the actual values of the two parameters are larger than the estimated values, the accuracy of the landing velocity is also significantly improved. This experiment demonstrates the fault tolerance and robustness of the guidance method to the landing position in the presence of parameter deviations. However, when the actual values of the two parameters are smaller than the estimated values, the landing velocity errors obtained by the guidance method will become larger. We refer to the analysis of Ref. [19] and learn that when the actual values of the two parameters are decreased, part of the mechanical energy of the rocket cannot be dissipated as expected; and when the thrust magnitude saturates, there is no additional energy to compensate for the undissipated mechanical energy, leading to the error of the optimal trajectory becomes larger. Subsequent research on guidance methods should attempt to address this issue.


### 6.2.4. Combination of Various Faults

Finally, we consider navigation system failures, control system failures (thrust magnitude has deviations), drag coefficient deviations, and atmospheric density deviations at the same time and carry out Monte Carlo simulations (200 cases). For the navigation system and control system, the failures are added in the same way as the previous experiments. Drag coefficient deviations and atmospheric density deviations are considered to be normally distributed, and their means are set to zero. The3σ values are, respectively, taken as CD~15%, ρ~15%. Results of terminal position and terminal velocity are shown in Figure 6. The maximum errors are (7.2551 m, 31.6108 m/s). This experiment proves that the guidance method can still keep the terminal states within acceptable limits under the combination of various faults.Figure 6
Monte Carlo results of terminal position and terminal velocity.


## 6.2.1. Faults Occur in the Navigation System

It is assumed that the navigation system has faults. Rather, the measurement of the real-time state of the rocket has the following errors:(30)rmeasure=rtrue±min10,rtrue×ε%×random0,1,vmeasure=vtrue±min1,vtrue×ε%×random0,1,where 10 and 1 represent the absolute errors of the rocket’s position and velocity, ε represents the maximum range of the relative error, random0,1 represents a random number generated from [0,1], rtrue,vtrue represents the actual state of the rocket, and rmeasure,vmeasure represents the state measured by the navigation system and fed back to the guidance system. We take the value of ε as 5, 10, and 15 for experiments, and the landing errors under these fault conditions are shown in Table 4.Table 4
Landing errors and masses under the condition that the navigation system has faults.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgε=52.46655.777923331.19ε=103.57505.830123329.48ε=154.06246.551423355.94It can be seen from Table4 that as the relative error of the measured state increases, the landing errors of the rocket also increase, but they remain within a small range. This experiment proves that the guidance method proposed in this paper is fault-tolerant and robust to navigation system faults.


## 6.2.2. Faults Occur in the Control System

It is assumed that the control system has two fault conditions: thrust cannot change continuously (that is, thrust is constant during each guidance cycle), and thrust magnitude has deviations. The trajectory optimization algorithm and the guidance method proposed in this paper are, respectively, applied to conduct simulations under the above two fault conditions, and their results are compared.(1) Thrust Cannot Change Continuously. For trajectory optimization, we assume that thrust is constant per second, and for guidance, thrust is constant during each guidance cycle. The landing errors and landing masses are shown in Table 5.Table 5
Landing errors and masses under the condition that thrust cannot change continuously.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgTrajectory optimization407.64014.806223081.57Guidance0.09140.314523136.54(2) Thrust Magnitude Has Deviations. It is assumed that there are deviations within ±5% of the thrust magnitude. In the experiment, random deviations within ±5% are applied to the thrust magnitude obtained by the trajectory optimization algorithm and the guidance method. Then, the deviated thrust is used to control the rocket. The landing errors and landing masses are shown in Table 6.Table 6
Landing errors and masses under the condition that thrust magnitude has deviations.
Landing position error/mLanding velocity error/(m/s)Landing mass/kgTrajectory optimization359.13676.393923137.13Guidance0.01850.020923133.57It can be seen from Tables5 and 6 that when the thrust cannot be continuously changed or the thrust magnitude has deviations, the landing position error and velocity error obtained by trajectory optimization are significantly larger than those obtained by guidance. The landing position errors of the guidance method under the two fault conditions remain at the same order of magnitude as that under the no-fault condition in Subsection 6.1, which demonstrates that the two fault conditions have little effect on the landing position. As for the landing velocity error, the guidance method can make it at the same order of magnitude as the no-fault condition when the thrust magnitude has deviations within ±5%. The error is one order of magnitude greater when the thrust cannot be continuously changed, but it still does not exceed 1 m/s. The simulation experiments demonstrate that the guidance method proposed in this paper is fault-tolerant and robust to control system faults.


## 6.2.3. Estimates of Parameters Have Deviations

It is assumed that the actual drag coefficient and atmospheric density differ by±15% from the estimated values. We apply the trajectory optimization algorithm and the guidance method proposed in this paper to conduct simulation experiments using the estimated parameters. Then, we perform numerical integration using the actual parameters to obtain real terminal states. The landing errors and landing masses of the two algorithms under the fault condition are shown in Tables 7 and 8, respectively, (where CD+15% indicates that the actual drag coefficient has a value of 1.15×CD).Table 7
Landing errors and masses of the trajectory optimization algorithm.
DeviationLanding position error/mLanding velocity error/(m/s)Landing mass/kgCD+15%289.59334.979023135.95CD−15%356.56896.259823135.95ρ+15%289.95844.979023135.95ρ−15%356.56896.259823135.95Table 8
Landing errors and masses of the guidance method.
DeviationLanding position error/mLanding velocity error/(m/s)Landing mass/kgCD+15%0.00350.023923230.63CD−15%4.162729.245623957.32ρ+15%0.00350.023923230.63ρ−15%4.162729.245623957.32It can be seen from Tables7 and 8 that the guidance method proposed in this paper can significantly improve the landing position accuracy in the case of the drag coefficient or the atmospheric density having deviations. When the actual values of the two parameters are larger than the estimated values, the accuracy of the landing velocity is also significantly improved. This experiment demonstrates the fault tolerance and robustness of the guidance method to the landing position in the presence of parameter deviations. However, when the actual values of the two parameters are smaller than the estimated values, the landing velocity errors obtained by the guidance method will become larger. We refer to the analysis of Ref. [19] and learn that when the actual values of the two parameters are decreased, part of the mechanical energy of the rocket cannot be dissipated as expected; and when the thrust magnitude saturates, there is no additional energy to compensate for the undissipated mechanical energy, leading to the error of the optimal trajectory becomes larger. Subsequent research on guidance methods should attempt to address this issue.


## 6.2.4. Combination of Various Faults

Finally, we consider navigation system failures, control system failures (thrust magnitude has deviations), drag coefficient deviations, and atmospheric density deviations at the same time and carry out Monte Carlo simulations (200 cases). For the navigation system and control system, the failures are added in the same way as the previous experiments. Drag coefficient deviations and atmospheric density deviations are considered to be normally distributed, and their means are set to zero. The3σ values are, respectively, taken as CD~15%, ρ~15%. Results of terminal position and terminal velocity are shown in Figure 6. The maximum errors are (7.2551 m, 31.6108 m/s). This experiment proves that the guidance method can still keep the terminal states within acceptable limits under the combination of various faults.Figure 6
Monte Carlo results of terminal position and terminal velocity.


## 7. Conclusions

In this paper, a fault-tolerant guidance method is proposed to realize online guidance of rocket vertical landing. The main contribution of this paper is that we propose a piecewise guidance algorithm. We first embed a trajectory optimization algorithm based on convex optimization in the MPC framework and then put forward a piecewise method to cope with the bang-bang characteristic of the thrust magnitude. An integrated GNC system is designed to enhance the fault-tolerance and robustness of the guidance method, which constitutes another contribution of this paper. Simulation experiments are conducted under conditions of no faults and deviations, navigation system failures, control system failures, drag coefficient deviations, and atmospheric density deviations, respectively, proving the reliability and fault-tolerance of the guidance method. The proposed method shows great potential for online use.



---

*Source: 9104823-2022-08-12.xml*