Powered Landing Control of Reusable Rockets Based on Softmax Double DDPG

Abstract

Multi-stage launch vehicles are currently the primary tool for humans to reach extraterrestrial space. The technology of recovering and reusing rockets can effectively shorten rocket launch cycles and reduce space launch costs. With the development of deep representation learning, reinforcement learning (RL) has become a robust learning framework capable of learning complex policies in high-dimensional environments. In this paper, a deep reinforcement learning-based reusable rocket landing control method is proposed. The mathematical process of reusable rocket landing is modelled by considering the aerodynamic drag, thrust, gravitational force, and Earth’s rotation during the landing process. A reward function is designed according to the rewards and penalties derived from mission accomplishment, terminal constraints, and landing performance. Based on this, the Softmax double deep deterministic policy gradient (SD3) deep RL method is applied to build a robust reusable rocket landing control method. In the constructed simulation environment, the proposed method can achieve convergent and robust control results, proving the effectiveness of the proposed method.

1. Introduction

With the expanding demand for launching satellites and spacecraft, expendable launch vehicles are becoming one of the main factors restricting the exploration of extraterrestrial space. Expendable launch vehicles discard critical components such as engine and navigation equipment directly after completing their missions, raising the cost of space access [1]. To tackle the above problems, reusable launch vehicles have gained widespread attention regarding their economy and flexibility [2,3,4,5].

Research in reusable launch vehicles has made continuous progress in recent years. As an extension of parachute landing technology, the parachute recovery method was widely studied and applied in various rockets [6,7,8,9]. Although it is a simple, reliable, and economical method, the unpowered descent process makes it difficult to predict the final landing point of the rocket. As a comparison, the powered landing control method has gradually become the preferred option considering the safety, controllability, and rapidity of the landing process. Academic research on powered landing generally apply offline or online trajectory optimization followed by various tracking control methods for pinpoint landing. Convex optimization is one of the typical methods to transform the optimal landing problem into a convex optimization problem where an optimal landing trajectory can be obtained [10,11,12]. Açıkmeşe et al. [10] introduced a convexification of the control constraints that was proven to be lossless and used interior methods of convex optimization to obtain optimal solutions of the original non-convex optimal control problem. Blackmore [11] combined convex optimization by CVXGEN [13] with powered divert guidance to optimize trajectory and eliminate error. While convex-optimization-based methods can achieve optimal solutions, they require the problem to be convex and often become computationally inefficient as the complexity of the problem increases. Besides convex optimization, many other numerical optimization methods can also be applied to solve the optimal landing problem. Sagliano et al. [14] applied a multivariate pseudospectral interpolation scheme to subspaces of a database of pre-computed trajectories to generate real-time capable entry guidance solutions. Seelbinder et al. [15] introduced parametric sensitivity analysis of non-linear programs for near-optimal control sequence generation. To take advantage of both convex optimization and pseudospectral methods, Wang et al. [16] adopted a pseudospectral-improved successive convexification algorithm to solve the trajectory optimization problem. The above-mentioned methods greatly enriched the solution diversity of the optimal landing problem, but they usually require extensive calculations and accurate modelling of the environment, which are difficult to satisfy in practice.

Leveraging the continuous development and breakthrough of deep learning technology, the idea of using deep neural networks (DNN) in the guidance, navigation, and control has been widely explored in the literature [17,18,19,20]. Among all deep learning technologies, supervised learning has been extensively discussed. DNNs have been applied to asteroid landing trajectory optimization [17] and first-stage recovery flight of reusable launch vehicles [18]. Furfaro et al. [19,20] designed a set of convolutional neural networks (CNN) and long-short-term memory (LSTM) architectures to handle the problem of the autonomous lunar landing. Although supervised learning exhibits strong fitting abilities for optimal landing controls, acquiring sufficient data required for the training process is very time-consuming [17,19]. An optimization method (e.g., GPOPS [19]) is usually employed before supervised learning to generate optimal state-control pairs as training samples (146,949 samples in [17] and 2,970,000 samples in [19]). Depending on the complexity of the optimization problem, a single sample can take anywhere from a few seconds to several hours to compute, and generating such a large number of samples consumes considerable time.

Deep reinforcement learning (DRL), on the other hand, uses trial-and-error interaction with the environment to optimize the policy instead of learning from enormous labelled data, thus making it useful in many fields such as chess [21], game [22,23], investment [24], and autonomous driving [25,26]. DRL has also drawn extensive attention to reusable launch vehicles. Ciabatti et al. [27] combined a deep deterministic policy gradient (DDPG) [28] algorithm with transfer learning to tackle the problem of planetary landing on unknown extra-terrestrial environments. Cheng et al. [29] derived an interactive reinforcement learning (RL) algorithm for the Moon fuel-optimal landing problem. Furfaro et al. [30] presented an adaptive guidance algorithm based on zero-effort-miss/zero-effort-velocity (ZEM/ZEV) and deterministic policy gradient [31]. Gaudet et al. [32] applied the principles of the deep Q network [33] into a six degrees-of-freedom Mars landing simulation environment, showing strong robustness.

In this paper, we present the final landing stage of a two-stage reusable rocket controlled by using DRL algorithms. In contrast with existing applied RL algorithms, we apply the Softmax double deep deterministic policy gradient (SD3) to keep the error between the value function estimated by the neural network and the true value function within a reasonable interval. This advantage is particularly useful when facing the coupling nature of the 3-DOF translational and 3-DOF rotational motion of the rocket, as the value function is difficult to accurately estimate in such scenarios. In addition, we designed a novel reward function by introducing a state-dependent Gaussian distribution. The landing result was evaluated by considering the landing position error, landing velocity, and landing impact related to the final acceleration. To verify the effectiveness of our proposed method, we built a three-dimensional simulation environment considering gravity, thrust, and aerodynamic drag to test the trained landing control policy. Our contribution are listed as follows:

1.

We introduced RL to obtain the optimal control commands for the reusable rocket landing problem. We apply the SD3 algorithm to simultaneously predict both the optimal thrust and its nozzle angle considering the three-axis translational and three-axis rotational dynamics of the rocket as well as environmental disturbances and uncertainties.

2.

We designed a novel reward function specific for the reusable rocket landing problem during the training of the RL models. We aim to give sufficient guidance to the rocket in terms of position, velocity, thrust, and attitude angle. We introduced Gaussian distribution rewards for position and attitude angle to give higher rewards under the condition of a higher probability of a successful landing.

2. Materials and Methods

2.1. Mathematical Modelling of the Reusable Rocket Landing

This paper focuses on the power landing stage of the reusable rocket landing process, where the position and attitude of the rocket are controlled by adjusting the thrust magnitude and nozzle angle. Four coordinate systems are used in this paper: the geocentric equatorial coordinate system, the launch point coordinate system, the body coordinate system of the rocket, and the velocity coordinate system of the rocket. Their definitions and relationships are derived in Appendix A.

The translational movement of the rocket is governed by the following translational dynamics equations expressed in the launch point coordinate system:

$$\left\{\begin{array}{cc}\hfill \dot{\mathit{r}}& =\mathit{v}\hfill \\ \hfill \dot{\mathit{v}}& =\frac{\mathit{T}+\mathit{D}}{m}+\mathit{g}-{\mathit{a}}_e-{\mathit{a}}_k\hfill \\ \hfill \dot{m}& =-\frac{1}{I_{sp}}\left|\right|\mathit{T}\left|\right|\hfill \end{array}\right.$$

(1)

where $\mathit{r}={\left[xyz\right]}^{\top }\in {\mathbb{R}}^3$ is the position vector, $\mathit{v}={\left[v_x{v}_y{v}_z\right]}^{\top }\in {\mathbb{R}}^3$ is the velocity vector, m is the mass of the rocket, $\mathit{T}$ is the thrust vector, $\mathit{D}$ is the aerodynamic drag vector, $\mathit{g}$ is the acceleration of gravity, and $I_{sp}$ is the specific impulse of the rocket. We define ${\omega }_{\mathit{e}}$ as the angular velocity of the Earth’s rotation, then ${\mathit{a}}_e={\mathit{\omega }}_e\times ({\mathit{\omega }}_e\times \mathit{r})$ is the convected acceleration, and ${\mathit{a}}_k=2{\mathit{\omega }}_e\times \frac{\delta \mathit{r}}{\delta t}$ is the Coriolis acceleration.

The acceleration of gravity $\mathit{g}$ can be decomposed into two components considering the oblateness of the Earth, i.e., the vertical component $g_r$ pointing from the Earth’s centre to the rocket’s centre of mass, and the horizontal component $g_m$ pointing from the rocket’s mass centre to North

$$\left\{\begin{array}{cc}\hfill g_r& =-\frac{\mu }{{\left|\right|\mathit{r}\left|\right|}^2}[1+\frac{3}{2}{J}_2{\left(\frac{a_e}{\left|\right|r\left|\right|}\right)}^2(1-5{sin}^2\mathsf{\Phi })]\hfill \\ \hfill g_m& =-\frac{3\mu }{{\left|\right|\mathit{r}\left|\right|}^2}{J}_2{\left(\frac{a_e}{\left|\right|r\left|\right|}\right)}^{2}sin\mathsf{\Phi }\hfill \end{array}\right.$$

where $\mu =3.986004\times {10}^{14}$ is the geocentric gravitational constant, $J_2=1.08263\times {10}^{-3}$ is the oblateness of the Earth, and $\mathsf{\Phi }$ is the latitude of the rocket.

The aerodynamic drag $\mathit{D}$ is described as

$$\mathit{D}=\frac{1}{2}\rho \left|\right|v\left|\right|S_{ref}{C}_{D}Mav$$

(2)

where $\rho$ is the atmospheric density determined by the height, v is the velocity of the rocket relative to the atmosphere, $S_{ref}$ is the reference area of the rocket, and $Ma$ is the Mach number. $C_D$ is the drag coefficient, which is a non-linear function of the velocity v. Detailed calculation of aerodynamic drag D is discussed in Appendix B.

The rotational movement of the rocket is governed by the attitude dynamics equation expressed as

$$\mathit{I}\frac{d{\mathit{\omega }}_b}{dt}+{\mathit{\omega }}_b\times \mathit{I}·{\mathit{\omega }}_b=\mathit{M}$$

(3)

where ${\mathit{\omega }}_b={\left[{\omega }_x{\omega }_y{\omega }_z\right]}^{\top }$ is the angular velocity of the rocket, $\mathit{M}$ is the external torque on the rocket, and $\mathit{I}$ is the inertia matrix of the rocket. The torque caused by aerodynamic drag is relatively small compared to the torque generated by the engine, it is neglected in the following discussion. The shape of the rocket is simplified and defined as a cylinder with radius r and height h. Therefore, the inertia matrix can be expressed as

$$\mathit{I}=\left[\begin{array}{ccc}\frac{1}{2}m{r}^2& 0& 0\\ 0& \frac{1}{12}m(3r^2+h^2)& 0\\ 0& 0& \frac{1}{12}m(3r^2+h^2)\end{array}\right].$$

(4)

2.2. SD3 Reinforcement Learning Algorithm

The RL problem can be formulated by a Markov decision process (MDP). We consider a standard RL setup consisting of an agent interacting with an environment in discrete time steps. In the powered landing problem, the agent is considered as the rocket itself. At each time step t, the agent receives an observation $s_t$, takes an action $a_t$ and receives a scalar reward $r_t(s_t,a_t)$. An agent’s behaviour is defined by a policy $\pi \left(s_t\right)$, which outputs the action $a_t$ based on the state $s_t$. The environment changes from one stage to another according to the transition probability $p\left(s_{t+1}\right|s_t,a_t)$. The return of an episode R is defined as the sum of discounted future reward with a discounting factor $\gamma \in [0,1]$

$$R_t=\sum _{i=0}^T{\gamma }^{i}r(s_i,a_i)$$

(5)

RL aims to learn a policy that maximizes the expected return from the starting distribution $J={\mathbb{E}}_{r_i,s_i\sim E,a_i\sim \pi }\left[R_1\right]$. Specifically, in the reusable rocket landing problem, policy $\pi$ refers to an actor network that can output thrust control commands based on states such as velocity and position, and return refers to the reward given by the environment. During the learning process, a rocket is randomly initialized repeatedly to land in a fixed landing zone under the control of policy $\pi$, and the policy is updated according to the environment-based rewards. The action-value function describes the expected return after taking an action $a_t$ in state $s_t$ and thereafter following policy $\pi$:

$$Q^{\pi }(s_t,a_t)={\mathbb{E}}_{r_{i\ge t},s_{i>t}\sim E,a_{i>t}\sim \pi }\left[r_t(s_t,a_t)\right]$$

(6)

According to the Bellman equation, the action-value function can be expressed as

$$Q^{\pi }(s_t,a_t)={\mathbb{E}}_{r_t,s_{t+1}\sim E}[r_t(s_t,a_t)+\gamma \underset{a_{t+1}}{max}{Q}^{\pi }(s_{t+1},a_{t+1})]$$

(7)

As the first method to extend DQN (deep Q network) to a continuous control space, DDPG is one of the standard algorithms for training an off-policy model-free RL network. DDPG learns both a policy network and a value function by decomposing them into an actor network $\pi \left(s\right|{\theta }^{\pi })$, target actor network ${\pi }^{\prime }\left(s\right|{\theta }^{{\pi }^{\prime }})$, Q network $Q(s,a|{\theta }^Q)$ and target Q network $Q^{\prime }(s,a|{\theta }^{Q^{\prime }})$. The parameterized actor network learns to determine the best action with the highest value estimates according to the critic network by policy gradient descent. The parameters of the target network are periodically replicated from the ontology to ensure the stability and convergence of training while alleviating the overestimation of Q. As an off-policy algorithm, it takes advantage of a replay buffer D to store past experiences and learn from them by randomly sampling from the replay buffer. The Q networks are updated by minimizing the loss L:

$$L=\underset{(s,a,r,s^{\prime })\sim D}{\mathbb{E}}{(Q(s,a|{\theta }^Q)-(r+\gamma \underset{a^{\prime }}{max}{Q}^{\prime }(s^{\prime },a^{\prime }|{\theta }^{Q^{\prime }})))}^2$$

(8)

The actor network is updated by policy gradient:

$${\nabla }_{{\theta }^{\pi }}J\approx \underset{(s,a,r,s^{\prime })\sim D}{\mathbb{E}}{\nabla }_{a}Q(s,a|{\theta }^Q)$$

(9)

TD3 (Twin-delayed deep deterministic policy gradient) is developed based on DDPG with three main improvements: clipped double Q learning, delayed policy updates and target policy smoothing.

• Clipped double Q
  As an extension of DQN in continuous control, DDPG inherits the over-estimation problem: the estimated Q value function is larger than the actual Q value function. To alleviate this problem, TD3 uses two Q value function networks at the same time, while in the update, the smaller value of the two networks is selected to calculate the updating target for the value function

  $$y(r,s^{\prime },d)=r+\gamma (1-d)\underset{i=1,2}{min}{Q}_{{\mathsf{\Phi }}_{i,arg}}(s^{\prime },a^{\prime }\left(s^{\prime }\right))$$
• Delayed policy updates
  Since an unstable value network can lead to a blind update of the policy network, it should be updated at a lower frequency than the value network. In this case, the error of the target networks is reduced to avoid updating the policy updates on high-error states which can cause divergent behaviour.
• Target policy smoothing
  Due to the deterministic policy, the update target of the Q network is more susceptible to the estimation error. TD3 introduces a regularization strategy for deep value learning in calculating the target Q, adding a small amount of random noise to the target policy and averaging over mini-batches.

  $$y=r+\gamma Q_{{\theta }^{\prime }}(s^{\prime },{\pi }_{{\varphi }^{\prime }}\left(s^{\prime }\right)+ϵ),ϵ\sim clip(\mathcal{N}(0,\sigma ),-c,c)$$

  (10)

This updated style can also lead to improvement in the stochastic domains with failure cases, making the estimation of Q more accurate and robust.

In this paper, we use the SD3 algorithm to alleviate the underestimation generated by directly applying the double Q-learning algorithm. A key component in TD3 is the clipped double Q-learning algorithm, which significantly reduces the overestimation but also leads to a large underestimation bias. To handle this problem, SD3 investigates the use of a Softmax operator in updating value functions in actor–critic methods for continuous control. It proposes to estimate the target value for the critic $Q_i$ by

$$y_i=r+\gamma {\mathcal{T}}_{SD3}\left(s^{\prime }\right)$$

(11)

where

$${\mathcal{T}}_{SD3}\left(s^{\prime }\right)=softma{x}_{\beta }\left({\widehat{Q}}_i(s^{\prime },·)\right)$$

(12)

$${\widehat{Q}}_i(s^{\prime },a^{\prime })=min\left(Q_i(s^{\prime },a^{\prime };{\theta }_i^{-}),Q_{-i}(s^{\prime },a^{\prime };{\theta }_{-i}^{-})\right)$$

(13)

By applying the Boltzmann Softmax operator for value function estimation, SD3 helps to smooth the optimization landscape and improve the underestimation bias of TD3. The SD3-based rocket recycling algorithm is shown in Algorithm 1.

Algorithm 1 SD3-based rocket recycling.

1: Initialize critic networks $Q_1$, $Q_2$, and actor networks ${\pi }_1$, ${\pi }_2$ with random parameters ${\theta }_1,{\theta }_2,{\varphi }_1,{\varphi }_2$ 2: Initialize target networks ${\theta }_1^{-}\leftarrow {\theta }_1,{\theta }_2^{-}\leftarrow {\theta }_2,{\varphi }_1^{-}\leftarrow {\varphi }_1,{\varphi }_2^{-}\leftarrow {\varphi }_2$ 3: Initialize replay buffer $\mathcal{B}$ 4: for episode $t=1$ to T do 5:    Observe a rocket state s 6:    Select a thruster action a with exploration noise $ϵ$ 7:    Execute action a, observe reward r, new state $s^i$ and done d 8:    Store transition tuple $(s,a,r,s^{\prime },d)$ in $\mathcal{B}$ 9:    while it’s time for update do 10:      Sample a mini-batch of N transitions $(s,a,r,s^{\prime },d)$ from $\mathcal{B}$ 11:      Sample K noises $ϵ\sim \mathcal{N}(0,\overline{\sigma })$ 12:      ${\widehat{a}}^{\prime }\leftarrow {\pi }_i(s^{\prime };{\varphi }_i^{-})+clip(ϵ,-c,c)$ 13:      $\widehat{Q}(s^{\prime },{\widehat{a}}^{\prime })\leftarrow {min}_{j=1,2}\left(Q_j(s^{\prime },{\widehat{a}}^{\prime };{\theta }_j^{-})\right)$ 14:      $softma{x}_{\beta }\left({\widehat{Q}}_i(s^{\prime },·)\right)\leftarrow {\mathbb{E}}_{{\widehat{a}}^{\prime }\sim p}\left[\frac{exp\left(\beta \widehat{Q}(s^{\prime },{\widehat{a}}^{\prime })\right)\widehat{Q}(s^{\prime },{\widehat{a}}^{\prime })}{p\left({\widehat{a}}^{\prime }\right)}\right]/{\mathbb{E}}_{{\widehat{a}}^{\prime }\sim p}\left[\frac{exp\left(\beta \widehat{Q}(s^{\prime },{\widehat{a}}^{\prime })\right)}{p\left({\widehat{a}}^{\prime }\right)}\right]$ 15:      $y_i\leftarrow r+\gamma (1-d)softma{x}_{\beta }\left({\widehat{Q}}_i(s^{\prime },·)\right)$ 16:      Update the critic ${\theta }_i$ according to Bellman loss: $\frac{1}{N}{\sum }_s{\left(Q_i(s,a;{\theta }_i)-y_i\right)}^2$ 17:      Update actor ${\varphi }_i$ by policy gradient: $\frac{1}{N}{\sum }_s\left[{\nabla }_{{\varphi }_i}\left(\pi (s;{\varphi }_i)\right){\nabla }_a{Q}_i(s,a;{\theta }_i){|}_{a=\pi (s;{\varphi }_i)}\right]$ 18:      Update target networks: ${\theta }_i^{-}\leftarrow \tau {\theta }_i+(1-\tau ){\theta }_i^{-},{\varphi }_i^{-}\leftarrow \tau {\varphi }_i+(1-\tau ){\varphi }_i^{-}$ 19:    end while 20: end for

2.3. Reinforcement Learning for Rocket Landing

In this section, we set up RL for reusable rocket landing, including action space, observation space, reward function, and network structure.

2.3.1. Action Space

The rocket is equipped with an engine with a deflectable nozzle, as is shown in Figure 1. The angle between $\mathit{T}$ and the plane $O_b{x}_b{z}_b$ is defined as ${\delta }_y$, and the angle between the projection of $\mathit{T}$ on plane $O_b{x}_b{z}_b$ and the opposite direction of $O_b{x}_b$ is defined as ${\delta }_z$. Therefore, the control variables of the rocket are $\left|\right|\mathit{T}\left|\right|,{\delta }_y$, and ${\delta }_z$.

To simplify the outputs of the agent for better numerical stability, the control variables are constrained to $[-1,1]$. The actual actions taken by the rocket can then be defined as $\left[a_1{a}_2{a}_3\right]$, and their relations with the original control variables are

$$\left\{\begin{array}{cc}\hfill \left|\right|\mathit{T}\left|\right|& =\frac{T_{max}-T_{min}}{2}{a}_1+\frac{T_{max}+T_{min}}{2}\hfill \\ \hfill {\delta }_y& ={\delta }_{max}{a}_2\hfill \\ \hfill {\delta }_z& ={\delta }_{max}{a}_3\hfill \end{array}\right.$$

(14)

The thrust vector in the body coordinate system $T_b$ can be expressed as

$${\mathit{T}}_b={\left[\begin{array}{ccc}\left|\right|\mathit{T}\left|\right|cos{\delta }_{y}cos{\delta }_z& \left|\right|\mathit{T}\left|\right|sin{\delta }_y& \left|\right|\mathit{T}\left|\right|cos{\delta }_{y}sin{\delta }_z\end{array}\right]}^{\top }$$

(15)

Since the torque caused by aerodynamic drag is relatively small, only the torque generated by thrust is considered in this paper. The shape of the rocket is defined as a cylinder with radius r and height h, then $\mathit{M}$ can be expressed as

$$\mathit{M}={\left[\begin{array}{ccc}0& \frac{1}{2}h\left|T\right||sin{\delta }_y& \frac{1}{2}h\left|\right|T\left|\right|cos{\delta }_{y}sin{\delta }_z\end{array}\right]}^{\top }$$

(16)

2.3.2. Observation Space

To make sure the agent receives enough information about the current states to make decisions during the reusable rocket landing process, we carefully design the information the agent can obtain from the environment, including normalized position information, velocity information, attitude angle information, and control information. Note that all the actions are obtained from the previous time step. Besides traditional velocity information, we introduce extra information for observation, i.e., the angle ${\theta }_v$ between the rocket’s velocity and a reference direction, to better stabilize the landing control process. The reference direction is defined to point from the rocket’s centre of mass to the ideal landing position. All angles except the roll angle are normalized by sine functions, whereas the roll angle requires both sine and cosine functions to cover $[0^{\circ },{360}^{\circ }]$. The observation space at time t can be expressed as

$$s_t={[\frac{r_x}{\left|\right|\mathit{r}\left|\right|}\frac{r_y}{\left|\right|\mathit{r}\left|\right|}\frac{r_z}{\left|\right|\mathit{r}\left|\right|}\frac{v_x}{\left|\right|\mathit{v}\left|\right|}\frac{v_y}{\left|\right|\mathit{v}\left|\right|}\frac{v_z}{\left|\right|\mathit{v}\left|\right|}sin{\theta }_{v}sin(\varphi +{90}^{\circ })sin\psi cos\psi sin\gamma a_1{a}_2{a}_3]}_t$$

(17)

2.3.3. Reward Function

The goal of the agent is to maximize the received reward, and a well-designed reward function can guide the agent to learn more robustly. We design a reward function including a successful landing reward, a failed landing penalty, a location-based reward, a velocity-based reward, and an angle-based reward.

• Reward for successful landing: The agent receives a reward when the rocket successfully lands at the target location. The conditions for successful landing can be expressed as $C_1$, which includes $\left|x\right|\le x_{lim}$, $0\le \left|y\right|\le y_{lim}$, $\left|z\right|\le z_{lim}$, $\parallel v_v\parallel <v_{vlim}$, $\parallel v_h\parallel <v_{hlim}$, $\left|\varphi \right|\le {\varphi }_{lim}$, $\left|\psi \right|\le {\psi }_{lim}$, and $\left|\gamma \right|\le {\gamma }_{lim}$. The reward can be defined as follows:

  $$R_{success}=\left\{\begin{array}{c}{\lambda }_{success},C_{1}satisfied\hfill \\ 0,C_{1}notsatisfied\hfill \end{array}\right.$$

  (18)

  where ${\lambda }_{success}$ is the positive constant hyperparameter, $x_{lim}=z_{lim}=10$ m are the length and width of the landing zone, $y_{lim}=0.5$ m is the upper bound of the landing height, $v_{vlim}=2$ m/s and $v_{hlm}=1$ m/s are the upper bounds of the vertical and horizontal velocity, respectively, ${\varphi }_{lim}={\psi }_{lim}={\gamma }_{lim}={10}^{\circ }$ are the upper bounds of the rocket’s landing Euler angle. Although there is no need to limit $\psi$, it is limited in the simulation.
• Penalty for failed landing: The agent will receive a penalty if the landing fails, including an excessive attitude angle or deviation from the landing zone, and the simulation will be terminated. A failed landing is that exact opposite of a successful landing. The penalty can be expressed as:

  $$R_{fail}=\left\{\begin{array}{c}{\lambda }_{fail},C_{1}notsatisfied\hfill \\ 0,C_{1}satisfied\hfill \end{array}\right.$$

  (19)
• Penalty for fuel consumption: The landing performance is evaluated by fuel consumption. The penalty aims to minimize fuel consumption and is defined by

  $$R_{fuel}={\lambda }_{fuel}\frac{1}{I_{sp}}\left|\right|\mathit{T}\left|\right|$$

  (20)
• Reward for accurate landing: The ideal landing location is set to be $\mathit{r}={[0,0,0]}^{\top }$. In order to encourage the rocket to land as close to the ideal landing location as possible, we construct the landing position reward function in the form of a Gaussian distribution,

  $$R_i=\frac{1}{\sqrt{2\pi }{\sigma }_i}exp\left(-\frac{r_i^2}{2{\sigma }_i^2}\right),i=x,y,z$$

  (21)

  where ${\sigma }_i$ determines the scope of the reward given, and the closer the rocket is to the ideal landing location, the bigger reward it will obtain. The total reward is:

  $$R_{pos}=\frac{1}{3}{\lambda }_{pos}\sum _{i=x,y,z}{R}_i$$

  (22)

  where ${\lambda }_{pos}$ is a hyperparameter to adjust the size of the reward.
• Reward for velocity direction: We construct a reward for velocity direction to stabilize the translational flight of the rocket. The reference direction is defined to point from the rocket’s centre of mass to the ideal location, and the angle ${\theta }_v$ between the current velocity and the reference direction is used as a guidance for learning. The smaller ${\theta }_v$ is, the bigger reward the rocket will obtain. The reward can be described as:

  $$R_{vel}={\lambda }_{vel}cos{\theta }_v$$

  (23)
• Reward for Euler angle: Similar to the reward for velocity direction, we want to stabilize the attitude motion of the rocket during landing and prefer to control the Euler angle of the rocket to arrive as close to our predefined reference attitude $\varphi =-{90}^{\circ }$, $\psi =0^{\circ }$ and $\gamma =0^{\circ }$. The reward function can be expressed as:

  $$R_{\varphi }=\frac{1}{\sqrt{2\pi }{\sigma }_{\varphi }}exp\left(-\frac{{(\varphi +{90}^{\circ })}^2}{2{\sigma }_{\varphi }^2}\right)$$

  (24)

  $$R_{\psi }=\frac{1}{\sqrt{2\pi }{\sigma }_{\psi }}exp\left(-\frac{{\psi }^2}{2{\sigma }_{\psi }^2}\right)$$

  (25)

  $$R_{\gamma }=\frac{1}{\sqrt{2\pi }{\sigma }_{\gamma }}exp\left(-\frac{{\gamma }^2}{2{\sigma }_{\gamma }^2}\right)$$

  (26)

  where ${\sigma }_{\varphi }$, ${\sigma }_{\psi }$ and ${\sigma }_{\gamma }$ limit the range of rewards given, and the total reward is:

  $$R_{angle}=\frac{1}{3}{\lambda }_{angle}\left(R_{\varphi }+R_{\psi }+R_{\gamma }\right)$$

  (27)

  where ${\lambda }_{angle}$ is a hyperparameter to adjust the size of the reward.
  Finally, the overall reward function is given in the form of a summation

  $$R=R_{success}+R_{fail}+R_{fuel}+R_{pos}+R_{vel}+R_{angle}$$

  (28)

2.3.4. Network Architecture

We apply the widely used actor–critic framework to construct the RL algorithm. The agent contains two types of networks: the actor network and the critic network. The actor network takes state $s_t$ as the input and outputs the action $a_t$. The critic network takes the concatenation of state $s_t$ and action $a_t$, and outputs the value $Q_{(s_t,a_t)}$. The network architecture is shown in Figure 2.

3. Results

In this section, we will demonstrate the simulation used to verify our proposed method and present the results.

3.1. Environmental Settings

The values of the simulation parameters are listed in

Table 1. Parameters of the simulation environment.

Parameters	Values
Launch azimuth A	${190}^{\circ }$
Latitude of launch point B	${41}^{\circ }$ N
Longitude of launch point $\lambda$	${100}^{\circ }$ E
Mass of rocket m	35,000 kg
Radius of rocket r	3.7 m
Height of rocket h	40 m
Minimum thrust $T_{min}$	250,000 N
Maximum thrust $T_{max}$	750,000 N
Reference area of rocket $S_{ref}$	10.75 m${}^2$
Specific impulse $I_{sp}$	345 s
Landing zone area	10 m × 10 m
Acceptable vertical landing velocity	$[-2.5,0]$ m/s
Acceptable translational landing velocity	$[0,1]$ m/s

. The overall training process is cut into episodes. In each episode, the position $r_x,r_y,r_z$, velocity $v_x,v_y,v_z$, and attitude angle $\varphi ,\psi ,\gamma$ of the rocket are initialized randomly according to

Table 2. Value ranges of the initial status of the rocket.

Initial States of the Rocket	Value Ranges
X-axis component of the initial position $r_x$	$[-500,500]$ m
Y-axis component of the initial position $r_y$	2500 m
Z-axis component of the initial position $r_z$	$[-500,500]$ m
X-axis component of the initial velocity $v_x$	$[-5,5]$ m/s
Y-axis component of the initial velocity $v_y$	$-230$ m/s
Z-axis component of the initial velocity $v_z$	$[-5,5]$ m/s
Initial pitch angle $\varphi$	${[-95,-85]}^{\circ }$
Initial yaw angle $\psi$	${[-5,5]}^{\circ }$
Initial roll angle $\gamma$	$0^{\circ }$

. The hyperparameters of the SD3 RL are shown in

Table 3. Hyperparameters of the SD3 RL algorithm.

Parameters	Values
Size of replay buffer $\left|\mathcal{B}\right|$	${10}^6$
Variance of policy noise $\overline{\sigma }$	$0.05$
Maximum noise c	$0.05$
Discount factor $\gamma$	$0.99$
Learning rate of actor networks	$5\times {10}^{-5}$
Learning rate of critic networks	$5\times {10}^{-5}$

.

3.2. Training Results

We apply DDPG, TD3, and SD3 separately to train our agents and compare their landing control performance. Their reward curves are shown in Figure 3.

We can see from the figure that rewards received by all three algorithms start to rise at around 3000 epochs, indicating a significant improvement in the agent’s control ability. All three rewards converge after 6000 episodes of training, meaning all three algorithms demonstrate a sufficient control ability. Among the three algorithms, SD3 gained the highest reward, reaching around 175.

To validate the effectiveness of the trained agent, we further select 100 initial simulation cases randomly according to Table 1 and Table 2 to conduct 100 landing tests. The landing trajectories are shown in Figure 4. According to the figure, all three algorithms are capable of guiding the rocket to approach the landing zone. However, when controlled by DDPG, the rocket failed to reach the ground successfully. When controlled by SD3, the rocket performs abnormal turning and lateral movement during the final landing stage. In comparison, the rocket controlled by SD3 managed to reach the landing zone vertically, smoothly and accurately.

The comparison of landing points is shown in Figure 5. The mean and standard deviation of the three sets of data are shown in the following

Table 4. Comparison of the landing points.

Landing Point Coordinates	DDPG	TD3	SD3
Average	$(5.534,-16.955)$ m	$(-5.795,-9.316)$ m	$(1.074,-5.133)$ m
Standard deviation	$(11.563,17.219)$ m	$(8.746,12.064)$ m	$(4.173,2.987)$ m
Best	$(-1.322,-1.549)$ m	$(-2.316,-0.248)$ m	$(-0.217,-0.534)$ m
Worst	$(-13.712,-46.027)$ m	$(-23.674,-25.946)$ m	$(-2.387,-9.976)$ m

. The DDPG agent failed to guide the rocket to land on the ground, the points in the graph are the vertical projection of the final position of the rocket on the ground surface. Although rockets controlled by TD3 agents successfully landed on the ground, all their landing points are out of the landing zone. The landing points controlled by SD3 achieved the best results compared to DDPG and TD3, with no landing points located outside the landing zone. The Comparison of landing points of the three methods are listed in

Table 5. Comparison of the reward function components for the different algorithms.

Reward Function Components	DDPG	TD3	SD3
$R_{fuel}$	$-0.147$	$-0.142$	$-0.118$
$R_{pos}$	$0.318$	$0.348$	$0.368$
$R_{vel}$	$0.441$	$0.454$	$0.469$
$R_{angle}$	$0.131$	$0.135$	$0.138$

. The results validate the effectiveness of the proposed SD3-based method.

Figure 6, Figure 7 and Figure 8 show the detailed landing status and control commands of the 26th trajectory controlled by the three algorithms, respectively.

4. Discussion

In RL, an optimal policy is linked to the optimal value function. When the reward function is fixed, the obtained landing policy is directly related to the estimation of the value function.

The DDPG method uses the Bellman equation that tends to overestimate the value function. Under this circumstance, the agents are more likely to explore areas with low rewards, thus making the training process unstable.

The TD3 method uses two value networks with the same structure but different learned parameters to estimate the value function. It takes the smaller value as the estimation result. This operation can avoid overestimation of the value function, but introduces the problem of underestimation. Therefore, although the training process of the TD3 algorithm is more stable compared to DDPG, the convergence speed of the algorithm is slower and the algorithm converges easily to areas with lower rewards, making the agents land poorly.

In comparison, the SD3 method adopts Softmax operation on the state-value function, making the error between the estimation result and true value within a controllable range. A better estimation of the state-value function leads to a better landing performance.

5. Conclusions

In this paper, an SD3-based RL method is proposed with a novel designed reward function to solve the reusable rocket landing control problem. Benefiting from the accurate estimation of the state-value function, the proposed method outperforms widely used DDPG and TD3 algorithms and has been proven to be effective through numerical simulation conducted on 100 randomly selected cases.