# Robust Linear Quadratic Regulator via Sliding Mode Guidance for Spacecraft Orbiting a Tumbling Asteroid

Authors: Peng Zhang; Tianhao Ma; Bo Zhao; Bo Dong; Yuanchun Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Aiming to ensure the stability of the spacecraft with multiuncertainties and mitigate the threat of the initial actuator saturation,
a Robust Linear Quadratic Regulator (RLQR) via sliding mode guidance (SMG) for orbiting a tumbling asteroid is proposed
in this paper. The orbital motion of the spacecraft near a tumbling asteroid is modelled in the body-fixed frame considering the
sun-relative effects, and the orbiting control problem is formulated as a stabilization of a nonlinear time-varying system. RLQR
based on the adaptive feedback linearization is proposed to stabilize the spacecraft orbiting with the uncertainties of the asteroid’s
rotation and gravitational field. In order to avoid the initial actuator saturation, SMG is applied to generate the transition process
trajectory of the closed-loop system. The effectiveness of the proposed control scheme is verified by the simulations of orbiting the
asteroid Toutatis 4179.


---

## Body

## 1. Introduction

The exploration of asteroids has become a hot topic of interests. Equilibria and periodic orbits near asteroids are usually highly unstable due to the irregular shapes and the complex rotation of the asteroids. Hence, the control of a spacecraft in closed proximity of an asteroid is among the most challenging technical problems in the exploration of asteroids [1].Many researches have paid great attention to the orbital control of the spacecrafts near asteroids. Sawai et al. [2] presented a control based on the one-dimensional altimetry measurements to stabilize a hovering spacecraft. Broschart and Scheeres [3] investigated the stability of realistic hovering control in the body-fixed and the inertial reference frames, respectively. Then, Broschart and Scheeres [4] proposed the sufficient conditions for a dead-band controller to bound spacecraft hovering motion in time-invariant Lagrangian dynamical systems. Wie [5] presented the dynamic modeling and control analysis of multiple gravity tractors in halo orbits. Furfaro et al. [6] investigated a novel closed-loop autonomous guidance law based on multiple sliding surfaces for the soft landing of the spacecraft on the designated point on the asteroid. Liu et al. [7] presented orbital control law for the spacecraft, which consists of PD controller and a nonsingular terminal sliding mode controller, to track the soft landing trajectory. Guelman [8] investigated a simple three-dimensional guidance law for the orbit transfer to a quasi-circular orbit about a rotating small celestial body using continuous thrust. In these previous contributions, asteroids are assumed to be a pure rotation configuration about the axis with its maximum inertia. However, the rotation of asteroids is very complex in practice and lots of them are time-varying [9]. The asteroid with nonprincipal axis rotational motion is so-called “tumbling asteroid” [10]. Nazari et al. [11] investigated the observer based body-frame hovering control over a tumbling asteroid, which is based on the time-varying LQR or the Lyapunov-Floquet transformation and time-invariant LQR, respectively. However, the rotation of the asteroid needs to be modelled in his control, which will cost a huge amount of telescope observing time to be determined [12], and the effect of the gravitational fields uncertainty is not analyzed. Hence, designing a robust control law for the spacecraft under the multiuncertainties of the asteroid’s gravitational field and rotation is urgent.The initial actuator saturation problem also needs to be considered. The output of the controller may have a high peak at the start, whose magnitude depends on the size of the initial error. As the continuous adjustable thrust is small, the control output may rise over the limitation of the thrust and threaten the stability of the system [13]. Aiming to avoid this problem, arranging the transition process is a good idea to acquire an acceptable controlled quantity. In many researches of the active disturbance rejection controllers, Tracking Differentiator (TD) is applied to arrange the transition process to avoid much larger control outputs [14–17]. However, the stability of the TD is difficult to be approved and the magnitude of the control output via TD is hardly adjustable. The sliding mode control is an important improvement of the control theory [18]. Based on it, the sliding mode guidance (SMG) is confirmed due to its good performance in the orbital control of the spacecraft near asteroids [19]. We apply the SMG to arrange the transition process of the controller. By doing this, the max magnitude of the control output becomes adjustable and the threat of the initial actuator saturation problem is reduced.In this paper, we proposed RLQR via the SMG for stabilizing the orbit of the spacecraft around a tumbling asteroid, which does not need to model the rotation of the asteroid. Firstly, the spacecraft orbital motion near a tumbling asteroid is modeled as a restricted three-body problem. RLQR based on the adaptive feedback linearization is proposed. The feedback linearization consists of a feedforward control, which is based on the spheric harmonic coefficient model of the asteroids gravitational field, and an adaptive compensator to ensure the robust stability against model uncertainties. Aiming to mitigate the threat of the initial actuator saturation, the SMG is applied to arrange the transition process of the RLQR. Simulations of orbiting the Toutatis 4769 are performed to verify the effectiveness of the proposed controller. The results of applying the RLQR with and without the SMG are compared to show the advantages.The rest of the paper is organized in the following form. In Section2, the orbital dynamic of spacecraft orbital motion is modeled and the problem formulation for control is proposed. In Section 3, the RLQR based on the adaptive feedback linearization is proposed. The SMG is applied to arrange the transition process of the proposed RLQR. In Section 4, simulations of orbiting the Toutatis 4769 are performed. Conclusions are drawn in Section 5.


## 2. Problem Formulation

For the relative orbital motion of a spacecraft near an asteroid, the dynamic model includes the nonspherical gravity field of the asteroid and the solar radiation pressure (SRP) [1]. In most previous researches, the orbital dynamic of the spacecraft was formulated into two regimes: the gravity dominated regime, in which the effects on the spacecraft from the sun were ignored, and T the solar dominated regime, in which the gravity of the asteroid was out of consideration. Few researches considered both of the effects from the sun and the asteroid [20]. In this section, the orbital dynamic of the spacecraft near a tumbling asteroid is modeled as a restricted three-body problem with the joint perturbations.


### 2.1. Orbital Dynamic of the Spacecraft near Tumbling Asteroids

Before the modeling, two relative frames need to be defined. As shown in Figure1, the inertial frame o-xIyIzI centered at the mass center of the asteroid. The xI-axis is parallel to the sun-line. The yI-axis is coinciding with the direction of the velocity vector of the asteroid. The xI-, yI-, and zI-axis compose the right-handed coordinate system. The body-fixed frame o-xyz fixes on asteroid with the origin coinciding with the mass center of the asteroid. The z, x, and y are coinciding with the axis of the asteroids maximum, minimum, and intermediate moment of inertia, respectively.Figure 1
Two relative frames of the orbiting.The orbital motion of spacecraft in the inertial frame can be expressed as(1)r¨I+Ω˙sI×rI+2ΩsI×r˙I+ΩsI×ΩsI×rI=gI+psI-μsrIraI3+μs1raI3-1rpI3raI+uI,where ΩsI is the orbital angular velocity vector of the asteroid, raI is the vector of the sun to the asteroid, rpI is the vector of the sun to the probe, rI is the vector of the asteroid to the probe, rI=rpI-raI, μs is the gravitational parameter of the sun, gI is the gravitational acceleration of the asteroid, psI is the solar radiation pressure on the probe, and uI is the vector of control.Assume that the asteroid orbits the sun in a circle; (1) can be simplified as(2)r¨I+2ΩsI×r˙I+ΩsI×ΩsI×rI=gI+psI-μsrIraI3+μsraI1raI3-1rpI3+uI.Transform (2) into the body-fixed frame; it can be expressed as(3)r¨+Ω˙a×r+2Ωa×r˙+Ωa×Ωa×r=fs+g+u,(4)fs=-2Ωs×r˙-2Ωs×Ωa×r-Ωs×Ωs×r-μsrra3+μsra1ra3-1rp3+ps,where Ωa is the rotation vector of the asteroid and fs is the sun-relative effect on the spacecraft.


#### 2.1.1. Sun-Relative Effects

DefineΩsI=[0,0,ωs]T, and the following equation can be established according to the assumption of asteroids circle orbit:(5)μsraI3=ωs2.Apply the Taylor expansion; an approximate linearization of (4) can be expressed as(6)fs=-2Ωs×r˙-2Ωs×Ωa×r+ps-Ωs×Ωs×r+ωs23raraTra2-1r.Considering that the main effect of the SRP arises from its first-order perturbations [21], psI can be defined as(7)psI=ξraIraI3,where ξ is a defined constant parameter of SRP.Define the transformation matrix from the inertial frame to the body-fixed frame asH(t)=hij3×3 and the rotation vector as Ωa=[ωx,ωy,ωz]T. Then, (6) can be rearranged as(8)fs=-ωsH0-10100000r˙-ωs2Mtr+ξra2H100,where(9)Mt=2ωsH-ωz0ωx0-ωzωy000-1-0-h33h23h330h13-h23h1302+3h13h23h33h13h23h33T.As usuallyωs≪ωi≪1,  i=x,y,z, and the elements of transformation matrix hij≤1,  i,j=1,2,3, hence (9) can be simplified as(10)Mt≈2ωsH-ωz0ωx0-ωzωy000.


#### 2.1.2. Asteroid’s Gravitational Field

The gravitational acceleration of the asteroid can be modelled by the gradient of the gravitational potential function of the asteroid as(11)g=∂Ur∂r.Different from the regular large objects, the bodies of asteroids are very irregular. Several methods have been investigated to approach the gravitational field of asteroids, such as the Spherical Harmonic Expansion Model (SHEM) [22], the polyhedron model [23], and the ellipsoidal harmonic expansion model [24]. The second-order SHEM, which models the asteroid’s perturbation potential with the most significant gravity coefficients and is widely used in the previous researches, is selected to formulate the nominal gravitational field in this paper. It can be expressed as(12)U0=μar1+r0r2C201-1.5cos2⁡δ+3C22cos2⁡δcos⁡2γ,where μa is the gravitational parameter of the asteroid, r0 is the normalizing radius, C20 and C22 are the spherical harmonic gravity coefficient, and r, δ, and γ are the radius, latitude, and longitude of the field point, respectively.Regard the asteroid as a three-axis ellipsoid with three main axesl1≥l2≥l3, and the spherical harmonic gravity coefficients C20 and C20 can be estimated by the three axes of the asteroid as(13)C20=152l32-l222l12-12,C22=1201-l22l12.Then, the real gravity of the asteroid can be formulated as(14)g=∂U0∂r+Δg,where Δg is the uncertainty of the asteroids gravitational model.Hypothesis 1.
Δ
g is bounded and the constant Γ exists, which satisfy the following conditions:(15)Γ=max⁡Δgr.


### 2.2. Problem Formulation for Controller Design

The stable hovering in the inertial frame is equivalent to a stable orbiting in the body-fixed frame. Hence, the stable orbiting condition in the body-fixed frame is expressed as(16)r˙=-Ωat×r.Define the error vector as(17)e=r˙+Ωat×r.Substitute (3) into (17); the time derivative of the error state is expressed as(18)e˙=-Ωat×e+fs+g+u.Then, the stable orbiting control problem can be formulated as the stabilization of the nonlinear system equation (18).The rotation of the asteroid is difficult to be modelled before the close proximity orbiting observation. Hence, the following assumptions need to be declared before the controller design:(i)
The rotation vectorΩa(t) is unmodeled, but its upper bound is known.(ii)
The current state of the rotation can be observed, so that the error vectore can be applied to the feedback.(iii)
The positionr and the velocity r˙ of the spacecraft in the body-fixed frame are provided by the navigation system.


## 2.1. Orbital Dynamic of the Spacecraft near Tumbling Asteroids

Before the modeling, two relative frames need to be defined. As shown in Figure1, the inertial frame o-xIyIzI centered at the mass center of the asteroid. The xI-axis is parallel to the sun-line. The yI-axis is coinciding with the direction of the velocity vector of the asteroid. The xI-, yI-, and zI-axis compose the right-handed coordinate system. The body-fixed frame o-xyz fixes on asteroid with the origin coinciding with the mass center of the asteroid. The z, x, and y are coinciding with the axis of the asteroids maximum, minimum, and intermediate moment of inertia, respectively.Figure 1
Two relative frames of the orbiting.The orbital motion of spacecraft in the inertial frame can be expressed as(1)r¨I+Ω˙sI×rI+2ΩsI×r˙I+ΩsI×ΩsI×rI=gI+psI-μsrIraI3+μs1raI3-1rpI3raI+uI,where ΩsI is the orbital angular velocity vector of the asteroid, raI is the vector of the sun to the asteroid, rpI is the vector of the sun to the probe, rI is the vector of the asteroid to the probe, rI=rpI-raI, μs is the gravitational parameter of the sun, gI is the gravitational acceleration of the asteroid, psI is the solar radiation pressure on the probe, and uI is the vector of control.Assume that the asteroid orbits the sun in a circle; (1) can be simplified as(2)r¨I+2ΩsI×r˙I+ΩsI×ΩsI×rI=gI+psI-μsrIraI3+μsraI1raI3-1rpI3+uI.Transform (2) into the body-fixed frame; it can be expressed as(3)r¨+Ω˙a×r+2Ωa×r˙+Ωa×Ωa×r=fs+g+u,(4)fs=-2Ωs×r˙-2Ωs×Ωa×r-Ωs×Ωs×r-μsrra3+μsra1ra3-1rp3+ps,where Ωa is the rotation vector of the asteroid and fs is the sun-relative effect on the spacecraft.


### 2.1.1. Sun-Relative Effects

DefineΩsI=[0,0,ωs]T, and the following equation can be established according to the assumption of asteroids circle orbit:(5)μsraI3=ωs2.Apply the Taylor expansion; an approximate linearization of (4) can be expressed as(6)fs=-2Ωs×r˙-2Ωs×Ωa×r+ps-Ωs×Ωs×r+ωs23raraTra2-1r.Considering that the main effect of the SRP arises from its first-order perturbations [21], psI can be defined as(7)psI=ξraIraI3,where ξ is a defined constant parameter of SRP.Define the transformation matrix from the inertial frame to the body-fixed frame asH(t)=hij3×3 and the rotation vector as Ωa=[ωx,ωy,ωz]T. Then, (6) can be rearranged as(8)fs=-ωsH0-10100000r˙-ωs2Mtr+ξra2H100,where(9)Mt=2ωsH-ωz0ωx0-ωzωy000-1-0-h33h23h330h13-h23h1302+3h13h23h33h13h23h33T.As usuallyωs≪ωi≪1,  i=x,y,z, and the elements of transformation matrix hij≤1,  i,j=1,2,3, hence (9) can be simplified as(10)Mt≈2ωsH-ωz0ωx0-ωzωy000.


### 2.1.2. Asteroid’s Gravitational Field

The gravitational acceleration of the asteroid can be modelled by the gradient of the gravitational potential function of the asteroid as(11)g=∂Ur∂r.Different from the regular large objects, the bodies of asteroids are very irregular. Several methods have been investigated to approach the gravitational field of asteroids, such as the Spherical Harmonic Expansion Model (SHEM) [22], the polyhedron model [23], and the ellipsoidal harmonic expansion model [24]. The second-order SHEM, which models the asteroid’s perturbation potential with the most significant gravity coefficients and is widely used in the previous researches, is selected to formulate the nominal gravitational field in this paper. It can be expressed as(12)U0=μar1+r0r2C201-1.5cos2⁡δ+3C22cos2⁡δcos⁡2γ,where μa is the gravitational parameter of the asteroid, r0 is the normalizing radius, C20 and C22 are the spherical harmonic gravity coefficient, and r, δ, and γ are the radius, latitude, and longitude of the field point, respectively.Regard the asteroid as a three-axis ellipsoid with three main axesl1≥l2≥l3, and the spherical harmonic gravity coefficients C20 and C20 can be estimated by the three axes of the asteroid as(13)C20=152l32-l222l12-12,C22=1201-l22l12.Then, the real gravity of the asteroid can be formulated as(14)g=∂U0∂r+Δg,where Δg is the uncertainty of the asteroids gravitational model.Hypothesis 1.
Δ
g is bounded and the constant Γ exists, which satisfy the following conditions:(15)Γ=max⁡Δgr.


## 2.1.1. Sun-Relative Effects

DefineΩsI=[0,0,ωs]T, and the following equation can be established according to the assumption of asteroids circle orbit:(5)μsraI3=ωs2.Apply the Taylor expansion; an approximate linearization of (4) can be expressed as(6)fs=-2Ωs×r˙-2Ωs×Ωa×r+ps-Ωs×Ωs×r+ωs23raraTra2-1r.Considering that the main effect of the SRP arises from its first-order perturbations [21], psI can be defined as(7)psI=ξraIraI3,where ξ is a defined constant parameter of SRP.Define the transformation matrix from the inertial frame to the body-fixed frame asH(t)=hij3×3 and the rotation vector as Ωa=[ωx,ωy,ωz]T. Then, (6) can be rearranged as(8)fs=-ωsH0-10100000r˙-ωs2Mtr+ξra2H100,where(9)Mt=2ωsH-ωz0ωx0-ωzωy000-1-0-h33h23h330h13-h23h1302+3h13h23h33h13h23h33T.As usuallyωs≪ωi≪1,  i=x,y,z, and the elements of transformation matrix hij≤1,  i,j=1,2,3, hence (9) can be simplified as(10)Mt≈2ωsH-ωz0ωx0-ωzωy000.


## 2.1.2. Asteroid’s Gravitational Field

The gravitational acceleration of the asteroid can be modelled by the gradient of the gravitational potential function of the asteroid as(11)g=∂Ur∂r.Different from the regular large objects, the bodies of asteroids are very irregular. Several methods have been investigated to approach the gravitational field of asteroids, such as the Spherical Harmonic Expansion Model (SHEM) [22], the polyhedron model [23], and the ellipsoidal harmonic expansion model [24]. The second-order SHEM, which models the asteroid’s perturbation potential with the most significant gravity coefficients and is widely used in the previous researches, is selected to formulate the nominal gravitational field in this paper. It can be expressed as(12)U0=μar1+r0r2C201-1.5cos2⁡δ+3C22cos2⁡δcos⁡2γ,where μa is the gravitational parameter of the asteroid, r0 is the normalizing radius, C20 and C22 are the spherical harmonic gravity coefficient, and r, δ, and γ are the radius, latitude, and longitude of the field point, respectively.Regard the asteroid as a three-axis ellipsoid with three main axesl1≥l2≥l3, and the spherical harmonic gravity coefficients C20 and C20 can be estimated by the three axes of the asteroid as(13)C20=152l32-l222l12-12,C22=1201-l22l12.Then, the real gravity of the asteroid can be formulated as(14)g=∂U0∂r+Δg,where Δg is the uncertainty of the asteroids gravitational model.Hypothesis 1.
Δ
g is bounded and the constant Γ exists, which satisfy the following conditions:(15)Γ=max⁡Δgr.


## 2.2. Problem Formulation for Controller Design

The stable hovering in the inertial frame is equivalent to a stable orbiting in the body-fixed frame. Hence, the stable orbiting condition in the body-fixed frame is expressed as(16)r˙=-Ωat×r.Define the error vector as(17)e=r˙+Ωat×r.Substitute (3) into (17); the time derivative of the error state is expressed as(18)e˙=-Ωat×e+fs+g+u.Then, the stable orbiting control problem can be formulated as the stabilization of the nonlinear system equation (18).The rotation of the asteroid is difficult to be modelled before the close proximity orbiting observation. Hence, the following assumptions need to be declared before the controller design:(i)
The rotation vectorΩa(t) is unmodeled, but its upper bound is known.(ii)
The current state of the rotation can be observed, so that the error vectore can be applied to the feedback.(iii)
The positionr and the velocity r˙ of the spacecraft in the body-fixed frame are provided by the navigation system.


## 3. Robust Linear Quadratic Regulator via Sliding Mode Guidance

### 3.1. Robust Linear Quadratic Regulator

The nonlinear system equation (18) can be transformed into (19). Consider(19)e˙=Ate+fs+g+u,(20)At=0ωzt-ωyt-ωzt0ωxtωyt-ωxt0=aij3×3.Because the rotation vectorΩa(t) is unmodeled, the time-varying matrix A(t) is uncertain. However, the upper bound of Ωa(t) is known. Hence, we can define the upper bound matrix of A(t) as(21)A+=aij+3×3,aij+=max⁡aijt.Design RLQR for the stabilization of the nonlinear system equation (19) as(22)u=uff+u1+ua.u
f
f is a feedforward control as(23)uff=-fs-∂U0∂r.u
1 is a linear quadratic regulator, which minimizes the cost function as (25). It can be expressed as(24)u1=-R-1Pe,(25)J=12∫t0∞eTtQeTt+u1TtRu1tdt,where Q and R are positive weight matrices and P is a positive symmetric matrix, which is the solution of the Riccati equation as(26)PA++A+TP-PR-1P+Q=0.u
a is an adaptive compensator, which can be expressed as(27)ua=-βPe,β˙=eTPPe,e∉δ0,e∈δβ0=0,where β is the adaptive parameter and δ is an adjustable field in the error space that contains the origin.Theorem 1.
Considering the spacecraft orbiting problem which is formulated as the stabilization of a nonlinear system as (19), define the error state as (18) and apply the robust LQR as (24). If Hypothesis 1 is satisfied, the closed-loop system is globally uniformly stable.Proof.
Define a Lyapunov function candidate as(28)V1=eTPe+12β∗-β2,where β∗ is a constant, which make all e∉δ satisfy the condition as(29)min⁡β∗κei≥Γ,i=x,y,z,κe=Pe=κex,κey,κezT.
The time derivative ofV1 can be expressed as(30)V˙1=Ate+fs+g+uTPe-β˙β∗-β+eTPAte+fs+g+u=Ate-R-1PeTPe+eTPAte-R-1Pe+ΔgTPe-βeTPPe-β˙β∗-β=eTAt-R-1PTPe+eTPAt-R-1Pe+ΔgTPe-β∗eTPPe-βeTPPe-β˙≤eTA+-R-1PTPe+eTPA+-R-1Pe+Γsign⁡e-β∗PePe-βeTPPe-β˙=Γsign⁡e-β∗κePe-eTPR-1P+Qe,e∉δ,where the function sign⁡e is defined as sign⁡ex,sign⁡ey,sign⁡(ez)T and the function κ is defined as κx,κy,κzT.
AsP, R, and Q are positive matrices, for all e∉δ,  V˙1<0. According to the Lyapunov stability theorem, the nonlinear system as (18) under the control as (22) is uniformly stable. The bound of the steady state error depends on the adjustable field δ.


### 3.2. RLQR via the Sliding Mode Guidance

The proposed RLQR could ensure the robust stability of the orbiting. However, the proposed controller has the potential threat from the actuator saturation. Firstly, depending on the error state, the control output may have a high peak at the start. Secondly, the gain of the adaptive compensator may be too conservative, which is much larger than necessary. As the adjustable continuous thrusters only could provide a small thrust, the saturation of the thrust may threaten the safety of the spacecraft. To avoid the initial actuator saturation problem of the RLQR, the sliding mode guidance is applied to arrange the transition process of the controller.


#### 3.2.1. The Sliding Mode Guidance

Define the outputs of the sliding mode guidance asɛ1 and ɛ2, and their initial values are(31)ɛ10=e0,ɛ20=0.The sliding mode guidance is expressed as(32)ɛ˙1=ɛ2,ɛ˙2=-k1ɛ2-k2sign⁡s,where s is the sliding mode function as s=ɛ2+k1ɛ1 and k1 and k2 are selected parameters of the SMG.Define a Lyapunov function candidate as(33)V2=12sTs.Then, the time derivative ofV2 is expressed as(34)V˙2=12s˙Ts+12sTs˙=ɛ˙2+k1ɛ˙1Ts=-k2sign⁡sTs≤0.Obviously,V˙2=0 only when s=0. According to the Lyapunov stability theorem, the SMG is globally asymptotically stable. When t→+∞, ɛ1→0 and ɛ2→0.


#### 3.2.2. The RLQR Based on the Virtual Error Vector

Define the virtual error vector as(35)υ=e-ɛ1.The time derivative of the virtual errorυ can be expressed as(36)υ˙=Atυ+g+fs+u~+Atɛ1-ɛ2.As the SMG is globally asymptotically stable and the matrixA(t) is bounded the constant Γ′ exists, which satisfy the condition as(37)Γ′=max⁡Atɛ1t.Then, design the RLQR based on the virtual error as(38)u~=-R-1+βIPυ-g0-fs+ɛ2,(39)β˙=υTPPυ,υ∉δ0,υ∈δ,β0=0.Theorem 2.
Considering the spacecraft orbiting problem which is formulated as the stabilization of a nonlinear system as (18), define the sliding mode guidance as (32) and the virtual error vector as (17) and (35), and apply the robust LQR based on the virtual error as (38). If Hypothesis 1 is satisfied, the closed-loop system is globally uniformly stable.Proof.
Define a Lyapunov function candidate as(40)V3=12υTPυ+12β-β′2+12sTs,where β′ is a constant, which make all υ∉δ satisfy the condition as(41)min⁡β′κυi≥Γ+Γ′,i=x,y,z,κυ=Pυ.
Then, the time derivative ofV3 is represented as(42)V˙3=υ˙TPυ+β˙β-β′+s˙Ts=-β′υTPPυ+Δg+Atɛ1TPυ-k2sign⁡sTs+υTAtP-2PR-1P+PAtυ-βυTPPυ-β˙≤Γ+Γ′sign⁡υ-β′PυPυ-k2sign⁡sTs-υTPR-1P+Qυ-βυTPPυ-β˙=-k2sign⁡sTs+Γ+Γ′sign⁡υ-β′PυPυ-υPR-1P+Qυ,υ∉δ.
AsP, Q, and R are positive matrices, for all υ∉δ, V˙3<0. Hence, when t→+∞, ɛ1→0, for all e∉δ, V˙3<0. According to the Lyapunov stability theorem, the RLQR via SMG could make the system as (18) uniformly stable. The bound of the steady state error is δ.


## 3.1. Robust Linear Quadratic Regulator

The nonlinear system equation (18) can be transformed into (19). Consider(19)e˙=Ate+fs+g+u,(20)At=0ωzt-ωyt-ωzt0ωxtωyt-ωxt0=aij3×3.Because the rotation vectorΩa(t) is unmodeled, the time-varying matrix A(t) is uncertain. However, the upper bound of Ωa(t) is known. Hence, we can define the upper bound matrix of A(t) as(21)A+=aij+3×3,aij+=max⁡aijt.Design RLQR for the stabilization of the nonlinear system equation (19) as(22)u=uff+u1+ua.u
f
f is a feedforward control as(23)uff=-fs-∂U0∂r.u
1 is a linear quadratic regulator, which minimizes the cost function as (25). It can be expressed as(24)u1=-R-1Pe,(25)J=12∫t0∞eTtQeTt+u1TtRu1tdt,where Q and R are positive weight matrices and P is a positive symmetric matrix, which is the solution of the Riccati equation as(26)PA++A+TP-PR-1P+Q=0.u
a is an adaptive compensator, which can be expressed as(27)ua=-βPe,β˙=eTPPe,e∉δ0,e∈δβ0=0,where β is the adaptive parameter and δ is an adjustable field in the error space that contains the origin.Theorem 1.
Considering the spacecraft orbiting problem which is formulated as the stabilization of a nonlinear system as (19), define the error state as (18) and apply the robust LQR as (24). If Hypothesis 1 is satisfied, the closed-loop system is globally uniformly stable.Proof.
Define a Lyapunov function candidate as(28)V1=eTPe+12β∗-β2,where β∗ is a constant, which make all e∉δ satisfy the condition as(29)min⁡β∗κei≥Γ,i=x,y,z,κe=Pe=κex,κey,κezT.
The time derivative ofV1 can be expressed as(30)V˙1=Ate+fs+g+uTPe-β˙β∗-β+eTPAte+fs+g+u=Ate-R-1PeTPe+eTPAte-R-1Pe+ΔgTPe-βeTPPe-β˙β∗-β=eTAt-R-1PTPe+eTPAt-R-1Pe+ΔgTPe-β∗eTPPe-βeTPPe-β˙≤eTA+-R-1PTPe+eTPA+-R-1Pe+Γsign⁡e-β∗PePe-βeTPPe-β˙=Γsign⁡e-β∗κePe-eTPR-1P+Qe,e∉δ,where the function sign⁡e is defined as sign⁡ex,sign⁡ey,sign⁡(ez)T and the function κ is defined as κx,κy,κzT.
AsP, R, and Q are positive matrices, for all e∉δ,  V˙1<0. According to the Lyapunov stability theorem, the nonlinear system as (18) under the control as (22) is uniformly stable. The bound of the steady state error depends on the adjustable field δ.


## 3.2. RLQR via the Sliding Mode Guidance

The proposed RLQR could ensure the robust stability of the orbiting. However, the proposed controller has the potential threat from the actuator saturation. Firstly, depending on the error state, the control output may have a high peak at the start. Secondly, the gain of the adaptive compensator may be too conservative, which is much larger than necessary. As the adjustable continuous thrusters only could provide a small thrust, the saturation of the thrust may threaten the safety of the spacecraft. To avoid the initial actuator saturation problem of the RLQR, the sliding mode guidance is applied to arrange the transition process of the controller.


### 3.2.1. The Sliding Mode Guidance

Define the outputs of the sliding mode guidance asɛ1 and ɛ2, and their initial values are(31)ɛ10=e0,ɛ20=0.The sliding mode guidance is expressed as(32)ɛ˙1=ɛ2,ɛ˙2=-k1ɛ2-k2sign⁡s,where s is the sliding mode function as s=ɛ2+k1ɛ1 and k1 and k2 are selected parameters of the SMG.Define a Lyapunov function candidate as(33)V2=12sTs.Then, the time derivative ofV2 is expressed as(34)V˙2=12s˙Ts+12sTs˙=ɛ˙2+k1ɛ˙1Ts=-k2sign⁡sTs≤0.Obviously,V˙2=0 only when s=0. According to the Lyapunov stability theorem, the SMG is globally asymptotically stable. When t→+∞, ɛ1→0 and ɛ2→0.


### 3.2.2. The RLQR Based on the Virtual Error Vector

Define the virtual error vector as(35)υ=e-ɛ1.The time derivative of the virtual errorυ can be expressed as(36)υ˙=Atυ+g+fs+u~+Atɛ1-ɛ2.As the SMG is globally asymptotically stable and the matrixA(t) is bounded the constant Γ′ exists, which satisfy the condition as(37)Γ′=max⁡Atɛ1t.Then, design the RLQR based on the virtual error as(38)u~=-R-1+βIPυ-g0-fs+ɛ2,(39)β˙=υTPPυ,υ∉δ0,υ∈δ,β0=0.Theorem 2.
Considering the spacecraft orbiting problem which is formulated as the stabilization of a nonlinear system as (18), define the sliding mode guidance as (32) and the virtual error vector as (17) and (35), and apply the robust LQR based on the virtual error as (38). If Hypothesis 1 is satisfied, the closed-loop system is globally uniformly stable.Proof.
Define a Lyapunov function candidate as(40)V3=12υTPυ+12β-β′2+12sTs,where β′ is a constant, which make all υ∉δ satisfy the condition as(41)min⁡β′κυi≥Γ+Γ′,i=x,y,z,κυ=Pυ.
Then, the time derivative ofV3 is represented as(42)V˙3=υ˙TPυ+β˙β-β′+s˙Ts=-β′υTPPυ+Δg+Atɛ1TPυ-k2sign⁡sTs+υTAtP-2PR-1P+PAtυ-βυTPPυ-β˙≤Γ+Γ′sign⁡υ-β′PυPυ-k2sign⁡sTs-υTPR-1P+Qυ-βυTPPυ-β˙=-k2sign⁡sTs+Γ+Γ′sign⁡υ-β′PυPυ-υPR-1P+Qυ,υ∉δ.
AsP, Q, and R are positive matrices, for all υ∉δ, V˙3<0. Hence, when t→+∞, ɛ1→0, for all e∉δ, V˙3<0. According to the Lyapunov stability theorem, the RLQR via SMG could make the system as (18) uniformly stable. The bound of the steady state error is δ.


## 3.2.1. The Sliding Mode Guidance

Define the outputs of the sliding mode guidance asɛ1 and ɛ2, and their initial values are(31)ɛ10=e0,ɛ20=0.The sliding mode guidance is expressed as(32)ɛ˙1=ɛ2,ɛ˙2=-k1ɛ2-k2sign⁡s,where s is the sliding mode function as s=ɛ2+k1ɛ1 and k1 and k2 are selected parameters of the SMG.Define a Lyapunov function candidate as(33)V2=12sTs.Then, the time derivative ofV2 is expressed as(34)V˙2=12s˙Ts+12sTs˙=ɛ˙2+k1ɛ˙1Ts=-k2sign⁡sTs≤0.Obviously,V˙2=0 only when s=0. According to the Lyapunov stability theorem, the SMG is globally asymptotically stable. When t→+∞, ɛ1→0 and ɛ2→0.


## 3.2.2. The RLQR Based on the Virtual Error Vector

Define the virtual error vector as(35)υ=e-ɛ1.The time derivative of the virtual errorυ can be expressed as(36)υ˙=Atυ+g+fs+u~+Atɛ1-ɛ2.As the SMG is globally asymptotically stable and the matrixA(t) is bounded the constant Γ′ exists, which satisfy the condition as(37)Γ′=max⁡Atɛ1t.Then, design the RLQR based on the virtual error as(38)u~=-R-1+βIPυ-g0-fs+ɛ2,(39)β˙=υTPPυ,υ∉δ0,υ∈δ,β0=0.Theorem 2.
Considering the spacecraft orbiting problem which is formulated as the stabilization of a nonlinear system as (18), define the sliding mode guidance as (32) and the virtual error vector as (17) and (35), and apply the robust LQR based on the virtual error as (38). If Hypothesis 1 is satisfied, the closed-loop system is globally uniformly stable.Proof.
Define a Lyapunov function candidate as(40)V3=12υTPυ+12β-β′2+12sTs,where β′ is a constant, which make all υ∉δ satisfy the condition as(41)min⁡β′κυi≥Γ+Γ′,i=x,y,z,κυ=Pυ.
Then, the time derivative ofV3 is represented as(42)V˙3=υ˙TPυ+β˙β-β′+s˙Ts=-β′υTPPυ+Δg+Atɛ1TPυ-k2sign⁡sTs+υTAtP-2PR-1P+PAtυ-βυTPPυ-β˙≤Γ+Γ′sign⁡υ-β′PυPυ-k2sign⁡sTs-υTPR-1P+Qυ-βυTPPυ-β˙=-k2sign⁡sTs+Γ+Γ′sign⁡υ-β′PυPυ-υPR-1P+Qυ,υ∉δ.
AsP, Q, and R are positive matrices, for all υ∉δ, V˙3<0. Hence, when t→+∞, ɛ1→0, for all e∉δ, V˙3<0. According to the Lyapunov stability theorem, the RLQR via SMG could make the system as (18) uniformly stable. The bound of the steady state error is δ.


## 4. Simulations

We specifically use the numerical data of the asteroid Toutatis 4179 in the simulation. To make the simulation more realistic, the gravitational field of the asteroid is modeled by the polyhedron method [23] based on Hudson’s 33996-face shape data [25], which can be expressed as (43). The second-order SHEM is applied as the nominal model in the controller in the same time:(43)U=12Gσ∑e∈edgesνe·Ee·νe·Le-12Gσ∑f∈facesνfFfνf·ψf,where G is the gravitational constant, σ is the constant density of the asteroid, νe is a vector from the field point to an arbitrary point on each edge, νf is a vector from the field point to an arbitrary point on each face, Ee is a dyad defined in terms of the face and edge normal vectors associated with each edge, Ff is the outer product of face normal vectors, Le is a logarithmic term expressing the potential of 1D straight wire, and ψf is the solid angle subtended by a face when viewed from the field point.According to the shape data, the maximum, intermediate, and minimum major axes of Toutatis 4179 arel1,l2,l3=4.5,2.4,1.9 km. The rotation of Toutatis in two periods is shown in Figure 2, which starts at Ωa(0)=(4.81,0,-12.53)×10-6 rad/s. The initial attitude of the asteroid expressed by 3-2-1 Euler angle is 0,π/2,π/3. The initial transformation matrix determined is by the initial attitude and the rotation of the asteroid. The gravitational parameters of Toutatis are shown in Table 1.Table 1
The gravitational parameters of the Toutatis.
Parameter
Value
r
0/m
4.500
e
2
μ
a/(km3/s2)
1.279
e
-
6
σ/(g/cm3)
2.100
e
0
ω
s/(rad/s)
4.937
e
-
8
ξ/(km3/s2)
3.750
e
6
r
a/m
3.806
e
11Figure 2
The rotation of the Toutatis.


### 4.1. The Orbiting under the Control of the RLQR

The weight matrixes of RLQR areR=10000I3×3 and Q=0.0001I3×3. The adjustable field of the error state is δ:eTe≤1×10-6. The initial position of the spacecraft in the body-fixed frame is r0=[8000,12000,11000]T m. The initial velocity of the spacecraft in the body-fixed frame is r˙0=[15.036,-15.606,3.8480]T m/s.As illustrated in Figure3, the error vector under RLQR control converges to zero quickly. Then, the errors keep near to the zero axes during the remaining times of the simulation and the system becomes stable. The outputs of the RLQR are shown in Figure 4. In the first several seconds, the control acceleration has a high peak, and the adaptive gain β grows quickly (as shown in Figure 5). Then, the control acceleration decreases and becomes very small that mainly restrains the effect of the sun and the gravity of the asteroid. The adaptive gain β is almost stable at 0.2213. In the simulation, the control cost of the steady-stable orbiting is about 0.0228 m/s per hour, which is acceptable for the deep space missions lasting months. Through these simulations, the effectiveness of the RLQR has been verified.Figure 3
The error under the control of RLQR.Figure 4
The control acceleration of RLQR.Figure 5
The adaptive gain of the RLQR.


### 4.2. The Orbiting under the RLQR via the SMG

The parameters of the SMG arek1=0.003 and k2=0.005. The adjustable field of the virtual error state is δ:υTυ≤1×10-6.The error states under the control of the RLQR via SMG are shown in Figure6. The convergence of the error states takes more time than the ones in Figure 3. The reduction of convergence speed is the price of the decrease of the control magnitude. As illustrated in Figure 7, the max value of the control acceleration is 0.046 m/s2, which is much less than the one in Figure 4. What is more, the max magnitude of the control acceleration can be adjusted by the parameters k1 and k2. Hence, the risk of the initial actuator saturation problem caused by the large initial error is mitigated. The gain of the adaptive compensator with the SMG is much smaller than the one without the SMG (see Figure 8). The control cost of the steady-stable orbiting is similar to the one without the SMG.Figure 6
The error under the control of the RLQR via SMG.Figure 7
The control acceleration of the RLQR with SMG.Figure 8
The comparison of the adaptive parameter in two kinds of controls.Through the simulations in Figures6–8, it was verified that the RLQR via the SMG could stabilize the orbiting and limit the magnitude of the control output. The threat of the initial actuator saturation has been mitigated.


### 4.3. The Proposed Controller under Different Magnitudes of Uncertainty

To analyze the effects of the gravity model’s uncertainties on the proposed algorithm, we did more simulations which apply no gravity model, the second-order SHEM, the fourth-order SHEM, and the polyhedron model (no uncertainty), respectively. Adding the order of the SHEM to the fourth, the gravitational potential of the asteroid can be expressed as(44)U0=μar1+r0r2C201-1.5cos2⁡δ+3C22cos2⁡δcos⁡2γ+r0r40.125C4035sin4⁡δ-30sin2⁡δ+3+7.5C42cos2⁡δ7sin2⁡δ-1cos⁡2γ+105C44cos2⁡δcos⁡4γ.The spheric harmonic coefficient can be determined as(45)C40=9l14+l24+24l34+6l12l22-24l12+l22l32140r04,C42=l12-l222l32-l12-l22280r04,C44=l12-l2222240r04.In the simulation, the initial position of the spacecraft is[8,10,11] km in the body-fixed frame. The initial velocity of the spacecraft is [21.3,-8.67,9.62] m/s in the body-fixed frame. The simulation results are shown as Figures 9–14.Figure 9
The error curve in thex-axis.Figure 10
The error curve in they-axis.Figure 11
The error curve in thez-axis.Figure 12
The control acceleration in thex-axis.Figure 13
The control acceleration in they-axis.Figure 14
The control acceleration in thez-axis.As shown in Figures9–11, the errors in three axes are converged near zero no matter which model applied. Also, the error curves controlled with different uncertainties have little difference. As shown in Figures 12–14, the control outputs are similar too.The reason of these results may be shown in Figures15–17. The asteroid’s gravitational acceleration is too small, which is in the magnitude of 10-5~10-6. Although we limited the control output, it is still much larger than the asteroid’s gravity.Figure 15
The gravitational acceleration in thex-axis.Figure 16
The gravitational acceleration in they-axis.Figure 17
The gravitational acceleration in thez-axis.As a conclusion, the uncertainty of the asteroid’s gravity has little effect on the proposed control algorithm. The proposed controller can be used in the orbiting around the asteroids which are not well mapped in terms of gravity field.


## 4.1. The Orbiting under the Control of the RLQR

The weight matrixes of RLQR areR=10000I3×3 and Q=0.0001I3×3. The adjustable field of the error state is δ:eTe≤1×10-6. The initial position of the spacecraft in the body-fixed frame is r0=[8000,12000,11000]T m. The initial velocity of the spacecraft in the body-fixed frame is r˙0=[15.036,-15.606,3.8480]T m/s.As illustrated in Figure3, the error vector under RLQR control converges to zero quickly. Then, the errors keep near to the zero axes during the remaining times of the simulation and the system becomes stable. The outputs of the RLQR are shown in Figure 4. In the first several seconds, the control acceleration has a high peak, and the adaptive gain β grows quickly (as shown in Figure 5). Then, the control acceleration decreases and becomes very small that mainly restrains the effect of the sun and the gravity of the asteroid. The adaptive gain β is almost stable at 0.2213. In the simulation, the control cost of the steady-stable orbiting is about 0.0228 m/s per hour, which is acceptable for the deep space missions lasting months. Through these simulations, the effectiveness of the RLQR has been verified.Figure 3
The error under the control of RLQR.Figure 4
The control acceleration of RLQR.Figure 5
The adaptive gain of the RLQR.


## 4.2. The Orbiting under the RLQR via the SMG

The parameters of the SMG arek1=0.003 and k2=0.005. The adjustable field of the virtual error state is δ:υTυ≤1×10-6.The error states under the control of the RLQR via SMG are shown in Figure6. The convergence of the error states takes more time than the ones in Figure 3. The reduction of convergence speed is the price of the decrease of the control magnitude. As illustrated in Figure 7, the max value of the control acceleration is 0.046 m/s2, which is much less than the one in Figure 4. What is more, the max magnitude of the control acceleration can be adjusted by the parameters k1 and k2. Hence, the risk of the initial actuator saturation problem caused by the large initial error is mitigated. The gain of the adaptive compensator with the SMG is much smaller than the one without the SMG (see Figure 8). The control cost of the steady-stable orbiting is similar to the one without the SMG.Figure 6
The error under the control of the RLQR via SMG.Figure 7
The control acceleration of the RLQR with SMG.Figure 8
The comparison of the adaptive parameter in two kinds of controls.Through the simulations in Figures6–8, it was verified that the RLQR via the SMG could stabilize the orbiting and limit the magnitude of the control output. The threat of the initial actuator saturation has been mitigated.


## 4.3. The Proposed Controller under Different Magnitudes of Uncertainty

To analyze the effects of the gravity model’s uncertainties on the proposed algorithm, we did more simulations which apply no gravity model, the second-order SHEM, the fourth-order SHEM, and the polyhedron model (no uncertainty), respectively. Adding the order of the SHEM to the fourth, the gravitational potential of the asteroid can be expressed as(44)U0=μar1+r0r2C201-1.5cos2⁡δ+3C22cos2⁡δcos⁡2γ+r0r40.125C4035sin4⁡δ-30sin2⁡δ+3+7.5C42cos2⁡δ7sin2⁡δ-1cos⁡2γ+105C44cos2⁡δcos⁡4γ.The spheric harmonic coefficient can be determined as(45)C40=9l14+l24+24l34+6l12l22-24l12+l22l32140r04,C42=l12-l222l32-l12-l22280r04,C44=l12-l2222240r04.In the simulation, the initial position of the spacecraft is[8,10,11] km in the body-fixed frame. The initial velocity of the spacecraft is [21.3,-8.67,9.62] m/s in the body-fixed frame. The simulation results are shown as Figures 9–14.Figure 9
The error curve in thex-axis.Figure 10
The error curve in they-axis.Figure 11
The error curve in thez-axis.Figure 12
The control acceleration in thex-axis.Figure 13
The control acceleration in they-axis.Figure 14
The control acceleration in thez-axis.As shown in Figures9–11, the errors in three axes are converged near zero no matter which model applied. Also, the error curves controlled with different uncertainties have little difference. As shown in Figures 12–14, the control outputs are similar too.The reason of these results may be shown in Figures15–17. The asteroid’s gravitational acceleration is too small, which is in the magnitude of 10-5~10-6. Although we limited the control output, it is still much larger than the asteroid’s gravity.Figure 15
The gravitational acceleration in thex-axis.Figure 16
The gravitational acceleration in they-axis.Figure 17
The gravitational acceleration in thez-axis.As a conclusion, the uncertainty of the asteroid’s gravity has little effect on the proposed control algorithm. The proposed controller can be used in the orbiting around the asteroids which are not well mapped in terms of gravity field.


## 5. Conclusion

In this paper, RLQR via SMG is proposed to ensure the robust stability of the spacecraft orbiting a tumbling asteroid. Based on the adaptive feedback linearization, the controller is able to suppress the effects of the multiuncertainties in the asteroids rotation and gravity model. Applying the SMG to arrange the transition process, the max magnitude of the control output is limited and threatening of the initial actuator saturation of the controller is mitigated. Also, the control cost of the steady-stable orbiting is acceptable for the deep space missions lasting months. As a consequence, the proposed control scheme can be properly selected to maintain spacecraft orbiting a tumbling asteroid and does not need to model the rotation of the asteroid.



---

*Source: 947843-2015-10-20.xml*