# Model Predictive Control to Autonomously Approach a Failed Spacecraft

Authors: Xun Wang; Zhaokui Wang; Yulin Zhang
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

In this study, a model predictive control (MPC) method is developed for a servicer spacecraft autonomously approaching a tumbling failed spacecraft at an ultraclose range. Flight safety and collision avoidance are basic requirements during the approach. Two types of a failed spacecraft with complex configurations are considered, and a double-ellipsoid composite envelope strategy is designed to model their keep-out zones. Given the keep-out zone of the servicer, two expanded ellipsoids are subsequently introduced to determine the collision and sufficient conditions for collision avoidance are derived by using the form of concave constraint. The tumbling motion of the target is considered, and a CW-based translational dynamics and derived attitude dynamics of the target are formulated to predict the motion of the docking point and keep-out zone. The MPC is formulated to drive the servicer tracking the docking point with collision avoidance and handle constraints including control input saturation and relative velocity bound. Convexification of the collision avoidance constraint and sequential convex programming are adopted for the implementation of MPC. Scenarios on the servicer with different initial positions approaching the target with different angular velocities are simulated, and the simulation results indicate that the proposed MPC method is effective.


---

## Body

## 1. Introduction

On-orbit services such as the repair of a failed spacecraft, spacecraft refueling, and spacecraft reorbiting are extremely important. A number of missions were conducted, such as the orbital express (OE) demonstration mission [1], the engineering test satellite VII (ETS-VII) project [2], and the spacecraft for the universal modification of orbits (SUMO) project [3]. Recently, on-orbit services for geosynchronous orbit (GEO) targets have attracted significant attention because of their specific importance. The Phoenix Program [4] and the Robotic Servicing of Geosynchronous Satellites Program (RSGS) [5] have thus been proposed in the last decade. An important stage, namely, “approaching targets at an ultraclose range,” is significant for the implementation of on-orbit services. However, it is challenging to ensure an effective and safe approach when the target exhibits a complex configuration and one that especially features uncontrolled tumbling motion.The study focuses on the control problem for a servicer spacecraft approaching a failed spacecraft with tumbling motion at an ultraclose range. Two types of the on-orbit failed spacecraft with complex configurations are considered. The first is a spacecraft with two large solar panels symmetrically mounted on its body, i.e., most GEO failed spacecraft [6, 7]. The second is a spacecraft with solar panels deployed on only one side where deployment of panels on the other side fails, i.e., SinoSat-II and TV-Sat-I in GEO [8]. The objective of the approach involves tracking a docking point (DP) that is fixed to the body frame of the target and is a few meters from its surface such as the docking mechanism on an apogee motor [5]. The implementation of the approach is conducive to the next step in on-orbit operation, i.e., catching the target with a space manipulator [9]. Flight safety is the major consideration in the approach. The keep-out zone that covers the outer surface of the target is always used for collision avoidance [10, 11]. Most related studies adopted a sphere envelope to model target’s keep-out zone [12, 13] although this is unsuitable for the approach problem discussed in the study. Consequently, a double-ellipsoid composite envelope strategy is designed to model the keep-out zone of the above two types of the failed spacecraft, and sufficient conditions for the collision are derived.Numerous guidance and control methods are used for spacecraft approaching at an ultraclose range [14–16]. However, most of the aforementioned methods do not consider collision avoidance with the target. Guidance methods based on the optimal control principle are widely used. Given the complex state or control constraints, it is difficult to solve analytic solutions to optimal control. In reference [17], the Gauss pseudospectral method was employed to solve the optimal rendezvous trajectory with certain constraints. In references [18–20], Pontryagin’s minimum principle was combined with numerical methods to formulate and solve the optimal rendezvous trajectory by considering collision avoidance. Another widely investigated method is the artificial potential function (APF) method. Dong et al. [21] introduced a potential function and two constrained zones to plan a safe rendezvous path. An adaptive control law based on a time-varying sliding manifold was subsequently used to track the desired path. Ge et al. [22] addressed the problem of docking with a tumbling target by using the APF method along with a sliding control. Tumbling motion was classified into three cases, and different safe boundaries were discussed in detail.Model predictive control is another method that should be considered. Given its working principle, it is extremely suitable for dealing with complex constraints such as multivariable state constraints and control input constraints. Several studies focused on investigating how to conduct close-range rendezvous and docking using the MPC method [23–25]. In reference [24], strategies for handling collision avoidance constraints, control input constraints, velocity constraints, and line-of-sight constraints were discussed in detail. Several studies also investigated the application of MPC to docking with a tumbling target [13, 26–28]. Under the condition, the docking point and the collision avoidance constraints are time variant. Collision avoidance is always modeled as a nonlinear and nonconvex constraint, and thus, it is difficult to ensure computational convergence and efficiency while solving the underlying optimization problem in MPC. Typically, the linearization of the collision avoidance constraint at the desired docking point is adopted [27], which is conducive to solving the optimization problem. However, the derived new collision avoidance zone is subsequently converted into a half-space, which is overly conservative, and the object in safe zones may be misjudged as corresponding to prohibited zones. Other effective methods for handling the collision avoidance constraints such as sequential convex programming and mixed integer programming have been investigated in spacecraft swarm mission and spacecraft rendezvous [29–31]. Based on the aforementioned studies, an MPC is developed for a servicer approaching a tumbling failed spacecraft. A novel collision discrimination method is proposed with the aim of handling the collision avoidance constraint with the modeled double-ellipsoid envelope of the target. Furthermore, a convexification method for the constraints and a sequential convex programming are adopted to implement the MPC.The remainder of the study is organized as follows: Section2 presents the mathematical formulation of relative translational dynamics and the attitude dynamics of the target. The control object and constraints on modeling for the approach are stated. Section 3 presents the design of the MPC controller and its solution in detail, and Section 4 provides the results of a numerical simulation to verify the performance of the proposed algorithms. The concluding remarks are discussed in Section 5.


## 2. Mathematical Formulation

In this section, the relative translational dynamics of the servicer and target are set along with the rotational dynamics of the target. Both the aforementioned dynamics form the basis of the state prediction and construction of time-variant constraints used in the design of the MPC.The relative geometrical relationship and coordinate systems are shown in Figure1.Figure 1
Geometrical relationship and coordinate systems.As shown in Figure1, r∈ℜ3 denotes the vector directed from the center of mass (c.m.) of the target to the c.m. of the servicer. Additionally, Rt and Rs denote the distances of the target and servicer, respectively, relative to the center of the Earth, and Re denotes the radius of the Earth.


### 2.1. Inertial Coordinate System∑I

The origin of the inertial coordinate systemoI−xIyIzI is centered on Earth, the zI-axis is along the rotational axis, and the xI-axis points toward the vernal equinox. The yI-axis completes the right-handed orthogonal.


### 2.2. Target Local Vertical-Local Horizontal (LVLH) Coordinate System∑t

The origin of the LVLH coordinate systemo−xtytzt lies on the c.m. of the target. The xt-axis is in the opposite direction to Earth’s center, the yt-axis is along the direction of flight, and the zt-axis completes the right-handed orthogonal.


### 2.3. Target Body-Fixed Coordinate System∑b

The origin of the target body-fixed coordinate systemo−xbybzb is in the c.m. of the target. Without loss of generality, it is assumed that xb, yb, and zb are aligned with its principal axes.


#### 2.3.1. Relative Translational Dynamics

It is assumed that the target is on a circular orbit and the distance between the servicer and target is extremely low. Given a short-period approach, the disturbances caused by solar pressure, atmospheric drag, and the effects of nonspherical gravity perturbation are neglected in both orbital and attitude motions. The Clohessy–Wiltshire (CW) equation is subsequently used to describe the relative translational dynamics [13]:
(1)x¨−2nẏ−3n2x=ux,y¨+2nẋ=uy,z¨+n2ż=uz,where ux,uy,uz denote the components of the control acceleration and are resolved into ∑t, n=μ/Rt3 denotes the orbital rate of the target with respect to ∑I, and μ denotes the gravitational constant.Eq. (1) is rewritten in the form of a state space model as follows:
(2)Ẋ=AX+BU,where X∈ℜ6 denotes the state vector, U∈ℜ3 denotes the control vector, and
(3)X=xyzẋẏż,A=0001000000100000013n20002n0000−2n0000−n2000,B=000000000100010001,U=uxuyuz.With respect to the implementation of the MPC, the discrete-time model of the relative translational dynamics is derived with a sampling periodTs as follows:
(4)Xk+1=AdXk+BdUk,where Xk and Uk denote the state and control vectors, respectively, at sampling instant k∈Z+. The matrices Ad∈ℜ6×6 and Bd∈ℜ6×3 are defined as follows [28]:
(5)Ad=eATs,Bd=∫0TseAτdτB.


#### 2.3.2. Attitude Dynamics of the Target

The attitude dynamics of the target are modeled to describe the relationship between∑t and ∑b. The attitude of the target is parameterized by using the rotation quaternion q=q0,q¯TT, where the first component denotes the scalar part and the other is a three-dimensional vector. The rotation matrix associated with q, denoted by Rtbq, is given as follows [32]:
(6)Rtbq=q02−q¯Tq¯Ι3+2q¯q¯T−2q0q¯×,where Rtbq transforms a vector from ∑t into ∑b, Ι3 denotes a three-dimensional identity matrix, and ∗× denotes a skew-symmetric matrix that represents the cross product operator.The angular velocityw of ∑b relative to ∑t is defined as follows:
(7)w=wb−wt,where wb and wt denote the angular velocities of ∑b and ∑t, respectively, relative to ∑I.The first time derivative of Eq. (7) in the ∑I is given as follows:
(8)dwdtI=dwbdtI−dwtdtI.Based on Coriolis’s theorem, the following expression is obtained:(9)dwdtI=dwdtt+wt×w.Given thatdwb/dt∣I=dwb/dt∣b and dwt/dt∣I=dwt/dt∣t, a combination of Eqs. (8) and (9) yields the following expression:
(10)dwbdtb=dwdtt+dwtdtt+wt×w.The attitude dynamics of the target are given by the following expression:(11)Itdwbdtb+wb×Itwb=Nt,where It and Nt denote the inertial tensors of the target and external torques, respectively. It is assumed that the target is on a circular orbit and environmental perturbations are neglected, and dwt/dt∣t=0 and Nt=0 are subsequently obtained, respectively.Eq. (10) is substituted into Eq. (11) to yield the following expression:
(12)Itẇ+wt×w+wt+w×It⋅wt+w=0,where wt, resolved into ∑t, is given as wt=0,0,nT.Subsequently, the change in target’s attitude is given as follows [32]:
(13)q̇=12q0Ι3+q¯×−q¯×w.According to Eq. (6), the rotation matrix Rtbk+1 at the sampling instant k+1 is a function of rotation quaternion qk+1. According to Eqs. (12) and (13), qk+1 can be obtained through solving the numerical integration with the given wk and qk. Thereby, Rtbk+1 can be represented as a nonlinear function expression:
(14)Rtbk+1=fwk,qk.The fourth-order Runge–Kutta algorithm is utilized to solve the numerical integration in Eqs. (12) and (13). It is assumed that there is no attitude maneuvering of the failed target and that there is no external torque applied to the target from the servicer. The offline calculation can be adopted to predict target’s attitude, which will be used in MPC to predict the relative translational states and servicer’s control input.


#### 2.3.3. Control Objective

The study focuses on two types of the on-orbit failed spacecraft. The first includes two large solar panels that are symmetrically mounted on the body as shown in Figure2. The second includes solar panels on only one side where deployment fails on the other side as shown in Figure 3.Figure 2
Failed spacecraft with two mounted solar panels.Figure 3
Failed spacecraft with solar panels deployed on only one side.Thexb-axis is along the center line of the solar panels, the yb-axis is along the symmetric axis of target’s body and is orthogonal to the xb-axis, and the zb-axis completes the dextral triad. In the study, the target is assumed as tumbling freely. It is also assumed that target’s DP is located along the yb-axis and the position vector of the DP is rDPb=0,−l,0T resolved in ∑b. The tracking error is defined as δrt=Rtbtrt−rDPb, where rt denotes servicer’s trajectories resolved in ∑t. Therefore, the objective of the approach procedure involves driving the servicer while tracking with the target DP, and this is expressed as follows:
(15)limt→∞δrt≤ε,where ε is a constant denoting the upper bound of tracking error.


#### 2.3.4. Constraint Modeling

With respect to the problem of the ultraclose approach, collision avoidance is the basic requirement. Given the aforementioned two types of the failed spacecraft, a double-ellipsoid composite envelope strategy is introduced to model target’s keep-out zone in which the spacecraft body and the solar panels are modeled as two different ellipsoids based on their sizes. Furthermore, a sphere envelope with radiusre is simplified as servicer’s keep-out zone.Intuitively, collision avoidance is attained when the envelopes of the target and the servicer do not intersect. In order to determine collision avoidance, expanded ellipsoids considering the sizes of envelopes of the target and the servicer are introduced, as shown in Figure4. The original ellipsoids denote target’s double-ellipsoid envelope. It is assumed that centers of the original ellipsoids for the first type of the failed spacecraft are located at its mass center. With respect to the original ellipsoid enveloping solar panels of the second type of the failed spacecraft, its center is located at the center of the solar panels.Figure 4
The double-ellipsoid composite envelope strategy for (a) the first type of a failed spacecraft and (b) the second type of a failed spacecraft.
(a)
(b)Furthermore, the scheme for discriminating collision avoidance between the original ellipsoid and the sphere is shown in Figure5.Figure 5
Scheme for discriminating collision avoidance between the original ellipsoid and the sphere.Without loss of generality, the original ellipsoid is as follows:(16)x−d2a2+y2b2+z2c2=1,where a, b, and c are constants designed to envelop the target body or the solar panels based on their sizes and a>b>c is satisfied. If the center of ellipsoid is located at target’s mass center, d=0 is satisfied. The parametric equation of the original ellipsoid is given as follows [17]:
(17)x=d+asinθcosφ,y=bsinθsinφ,z=ccosθ,where θ∈0,π and φ∈0,2π.The expanded ellipsoid envelope is as follows:(18)x−d2a+are/c2+y2b+bre/c2+z2c+re2=1.It should be noted that the value of the minimum half-axis of the expanded ellipsoid is the summation of the radius of the sphere and the value of the minimum half-axis of the original ellipsoid. The other two half-axes of expanded ellipsoid are augmented based on the product of the radius and the ratio between the corresponding half-axis and minimum half-axis.The parametric equation of the expanded ellipsoid is given as follows:(19)x=d+a+arecsinθcosφ,y=b+brecsinθsinφ,z=c+recosθ.Proposition 1.
Given the original ellipsoid in Eq. (16) and the sphere with radius re, the sufficient condition for collision avoidance between the original ellipsoid and the sphere is defined as follows:
(20)x−d2a+are/c2+y2b+bre/c2+z2c+re2−1>0,where x,y,z denote the coordinates of the sphere center. Evidently, Eq. (20) indicates that the center of the sphere is outside the expanded ellipsoid shown in Eq. (18).Proof.
Given a pointp0 using spherical coordinates θ,φ describing its direction, there are two other corresponding points located at the original ellipsoidal surface and the expanded ellipsoidal surface (points p1 and p2, as shown in Figure 5), and their coordinates are given as follows:
(21)p1:x1=d+asinθcosφ,y1=bsinθsinφ,z1=ccosθ,p2:x2=d+a+arecsinθcosφ,y2=b+brecsinθsinφ,z2=c+recosθ.
Subsequently, the tangent planes of the two points (planesL1 and L2, as shown in Figure 5) are given as follows:
(22)L1:sinθcosφax+sinθsinφby+cosθcz−sinθcosφad−1=0,L2:sinθcosφax+sinθsinφby+cosθcz−sinθcosφad−1−rec=0.
According to Eq. (22), it is noted that the two tangent planes are parallel. Subsequently, the minimum distance from the point p2 to the tangent plane L1 is equal to the distance between tangent planes L2 and L1, and this is given as follows:
(23)rc=rec2/a2sin2θcos2φ+c2/b2sin2θsin2φ+cos2θ.
Given thata>b>c is defined, rc≥re is subsequently obtained and it denotes that the minimum distance between the original ellipsoid and the expanded ellipsoid exceeds or is equal to the radius of the sphere. Thus, when the center of the sphere lies outside the expanded ellipsoid, the collision avoidance of the sphere and the original ellipsoid is satisfied. The proof is completed.Remark 1.
Given the double-ellipsoid envelope modeling target’s keep-out zone, the collision avoidance condition for the target and servicer is the collision avoidance between the sphere and the two original ellipsoids. The modeled ellipsoid envelope is fixed in target’s body frame. Therefore, the mathematical description of collision avoidance between the servicer and the target can be conducted in target’s body frame. Although this is a sufficient but not a necessary condition to determine collision avoidance, it is derived by using the form of concave constraint, which is conducive to the implementation of MPC as subsequently discussed.In practice, all thrusters are designed with limited control capability. The infinite norm of vector is applied to describe the thrust saturation constraint, and this is given as follows [24]:
(24)U∞≤umax,where umax denotes the maximum thrust output.The relative velocity constraint is also considered to ensure that the servicer exhibits the capability to adjust the thrust output and avoid collision in case of emergency, and this is extremely important in actual on-orbit approaches. The velocity constraint is expressed as follows [27]:
(25)ẋ≤ẋmax,ẏ≤ẏmax,ż≤żmax,where ẋmax, ẏmax, and żmax denote the components of the maximum approach velocity in each direction.


## 2.1. Inertial Coordinate System∑I

The origin of the inertial coordinate systemoI−xIyIzI is centered on Earth, the zI-axis is along the rotational axis, and the xI-axis points toward the vernal equinox. The yI-axis completes the right-handed orthogonal.


## 2.2. Target Local Vertical-Local Horizontal (LVLH) Coordinate System∑t

The origin of the LVLH coordinate systemo−xtytzt lies on the c.m. of the target. The xt-axis is in the opposite direction to Earth’s center, the yt-axis is along the direction of flight, and the zt-axis completes the right-handed orthogonal.


## 2.3. Target Body-Fixed Coordinate System∑b

The origin of the target body-fixed coordinate systemo−xbybzb is in the c.m. of the target. Without loss of generality, it is assumed that xb, yb, and zb are aligned with its principal axes.


### 2.3.1. Relative Translational Dynamics

It is assumed that the target is on a circular orbit and the distance between the servicer and target is extremely low. Given a short-period approach, the disturbances caused by solar pressure, atmospheric drag, and the effects of nonspherical gravity perturbation are neglected in both orbital and attitude motions. The Clohessy–Wiltshire (CW) equation is subsequently used to describe the relative translational dynamics [13]:
(1)x¨−2nẏ−3n2x=ux,y¨+2nẋ=uy,z¨+n2ż=uz,where ux,uy,uz denote the components of the control acceleration and are resolved into ∑t, n=μ/Rt3 denotes the orbital rate of the target with respect to ∑I, and μ denotes the gravitational constant.Eq. (1) is rewritten in the form of a state space model as follows:
(2)Ẋ=AX+BU,where X∈ℜ6 denotes the state vector, U∈ℜ3 denotes the control vector, and
(3)X=xyzẋẏż,A=0001000000100000013n20002n0000−2n0000−n2000,B=000000000100010001,U=uxuyuz.With respect to the implementation of the MPC, the discrete-time model of the relative translational dynamics is derived with a sampling periodTs as follows:
(4)Xk+1=AdXk+BdUk,where Xk and Uk denote the state and control vectors, respectively, at sampling instant k∈Z+. The matrices Ad∈ℜ6×6 and Bd∈ℜ6×3 are defined as follows [28]:
(5)Ad=eATs,Bd=∫0TseAτdτB.


### 2.3.2. Attitude Dynamics of the Target

The attitude dynamics of the target are modeled to describe the relationship between∑t and ∑b. The attitude of the target is parameterized by using the rotation quaternion q=q0,q¯TT, where the first component denotes the scalar part and the other is a three-dimensional vector. The rotation matrix associated with q, denoted by Rtbq, is given as follows [32]:
(6)Rtbq=q02−q¯Tq¯Ι3+2q¯q¯T−2q0q¯×,where Rtbq transforms a vector from ∑t into ∑b, Ι3 denotes a three-dimensional identity matrix, and ∗× denotes a skew-symmetric matrix that represents the cross product operator.The angular velocityw of ∑b relative to ∑t is defined as follows:
(7)w=wb−wt,where wb and wt denote the angular velocities of ∑b and ∑t, respectively, relative to ∑I.The first time derivative of Eq. (7) in the ∑I is given as follows:
(8)dwdtI=dwbdtI−dwtdtI.Based on Coriolis’s theorem, the following expression is obtained:(9)dwdtI=dwdtt+wt×w.Given thatdwb/dt∣I=dwb/dt∣b and dwt/dt∣I=dwt/dt∣t, a combination of Eqs. (8) and (9) yields the following expression:
(10)dwbdtb=dwdtt+dwtdtt+wt×w.The attitude dynamics of the target are given by the following expression:(11)Itdwbdtb+wb×Itwb=Nt,where It and Nt denote the inertial tensors of the target and external torques, respectively. It is assumed that the target is on a circular orbit and environmental perturbations are neglected, and dwt/dt∣t=0 and Nt=0 are subsequently obtained, respectively.Eq. (10) is substituted into Eq. (11) to yield the following expression:
(12)Itẇ+wt×w+wt+w×It⋅wt+w=0,where wt, resolved into ∑t, is given as wt=0,0,nT.Subsequently, the change in target’s attitude is given as follows [32]:
(13)q̇=12q0Ι3+q¯×−q¯×w.According to Eq. (6), the rotation matrix Rtbk+1 at the sampling instant k+1 is a function of rotation quaternion qk+1. According to Eqs. (12) and (13), qk+1 can be obtained through solving the numerical integration with the given wk and qk. Thereby, Rtbk+1 can be represented as a nonlinear function expression:
(14)Rtbk+1=fwk,qk.The fourth-order Runge–Kutta algorithm is utilized to solve the numerical integration in Eqs. (12) and (13). It is assumed that there is no attitude maneuvering of the failed target and that there is no external torque applied to the target from the servicer. The offline calculation can be adopted to predict target’s attitude, which will be used in MPC to predict the relative translational states and servicer’s control input.


### 2.3.3. Control Objective

The study focuses on two types of the on-orbit failed spacecraft. The first includes two large solar panels that are symmetrically mounted on the body as shown in Figure2. The second includes solar panels on only one side where deployment fails on the other side as shown in Figure 3.Figure 2
Failed spacecraft with two mounted solar panels.Figure 3
Failed spacecraft with solar panels deployed on only one side.Thexb-axis is along the center line of the solar panels, the yb-axis is along the symmetric axis of target’s body and is orthogonal to the xb-axis, and the zb-axis completes the dextral triad. In the study, the target is assumed as tumbling freely. It is also assumed that target’s DP is located along the yb-axis and the position vector of the DP is rDPb=0,−l,0T resolved in ∑b. The tracking error is defined as δrt=Rtbtrt−rDPb, where rt denotes servicer’s trajectories resolved in ∑t. Therefore, the objective of the approach procedure involves driving the servicer while tracking with the target DP, and this is expressed as follows:
(15)limt→∞δrt≤ε,where ε is a constant denoting the upper bound of tracking error.


### 2.3.4. Constraint Modeling

With respect to the problem of the ultraclose approach, collision avoidance is the basic requirement. Given the aforementioned two types of the failed spacecraft, a double-ellipsoid composite envelope strategy is introduced to model target’s keep-out zone in which the spacecraft body and the solar panels are modeled as two different ellipsoids based on their sizes. Furthermore, a sphere envelope with radiusre is simplified as servicer’s keep-out zone.Intuitively, collision avoidance is attained when the envelopes of the target and the servicer do not intersect. In order to determine collision avoidance, expanded ellipsoids considering the sizes of envelopes of the target and the servicer are introduced, as shown in Figure4. The original ellipsoids denote target’s double-ellipsoid envelope. It is assumed that centers of the original ellipsoids for the first type of the failed spacecraft are located at its mass center. With respect to the original ellipsoid enveloping solar panels of the second type of the failed spacecraft, its center is located at the center of the solar panels.Figure 4
The double-ellipsoid composite envelope strategy for (a) the first type of a failed spacecraft and (b) the second type of a failed spacecraft.
(a)
(b)Furthermore, the scheme for discriminating collision avoidance between the original ellipsoid and the sphere is shown in Figure5.Figure 5
Scheme for discriminating collision avoidance between the original ellipsoid and the sphere.Without loss of generality, the original ellipsoid is as follows:(16)x−d2a2+y2b2+z2c2=1,where a, b, and c are constants designed to envelop the target body or the solar panels based on their sizes and a>b>c is satisfied. If the center of ellipsoid is located at target’s mass center, d=0 is satisfied. The parametric equation of the original ellipsoid is given as follows [17]:
(17)x=d+asinθcosφ,y=bsinθsinφ,z=ccosθ,where θ∈0,π and φ∈0,2π.The expanded ellipsoid envelope is as follows:(18)x−d2a+are/c2+y2b+bre/c2+z2c+re2=1.It should be noted that the value of the minimum half-axis of the expanded ellipsoid is the summation of the radius of the sphere and the value of the minimum half-axis of the original ellipsoid. The other two half-axes of expanded ellipsoid are augmented based on the product of the radius and the ratio between the corresponding half-axis and minimum half-axis.The parametric equation of the expanded ellipsoid is given as follows:(19)x=d+a+arecsinθcosφ,y=b+brecsinθsinφ,z=c+recosθ.Proposition 1.
Given the original ellipsoid in Eq. (16) and the sphere with radius re, the sufficient condition for collision avoidance between the original ellipsoid and the sphere is defined as follows:
(20)x−d2a+are/c2+y2b+bre/c2+z2c+re2−1>0,where x,y,z denote the coordinates of the sphere center. Evidently, Eq. (20) indicates that the center of the sphere is outside the expanded ellipsoid shown in Eq. (18).Proof.
Given a pointp0 using spherical coordinates θ,φ describing its direction, there are two other corresponding points located at the original ellipsoidal surface and the expanded ellipsoidal surface (points p1 and p2, as shown in Figure 5), and their coordinates are given as follows:
(21)p1:x1=d+asinθcosφ,y1=bsinθsinφ,z1=ccosθ,p2:x2=d+a+arecsinθcosφ,y2=b+brecsinθsinφ,z2=c+recosθ.
Subsequently, the tangent planes of the two points (planesL1 and L2, as shown in Figure 5) are given as follows:
(22)L1:sinθcosφax+sinθsinφby+cosθcz−sinθcosφad−1=0,L2:sinθcosφax+sinθsinφby+cosθcz−sinθcosφad−1−rec=0.
According to Eq. (22), it is noted that the two tangent planes are parallel. Subsequently, the minimum distance from the point p2 to the tangent plane L1 is equal to the distance between tangent planes L2 and L1, and this is given as follows:
(23)rc=rec2/a2sin2θcos2φ+c2/b2sin2θsin2φ+cos2θ.
Given thata>b>c is defined, rc≥re is subsequently obtained and it denotes that the minimum distance between the original ellipsoid and the expanded ellipsoid exceeds or is equal to the radius of the sphere. Thus, when the center of the sphere lies outside the expanded ellipsoid, the collision avoidance of the sphere and the original ellipsoid is satisfied. The proof is completed.Remark 1.
Given the double-ellipsoid envelope modeling target’s keep-out zone, the collision avoidance condition for the target and servicer is the collision avoidance between the sphere and the two original ellipsoids. The modeled ellipsoid envelope is fixed in target’s body frame. Therefore, the mathematical description of collision avoidance between the servicer and the target can be conducted in target’s body frame. Although this is a sufficient but not a necessary condition to determine collision avoidance, it is derived by using the form of concave constraint, which is conducive to the implementation of MPC as subsequently discussed.In practice, all thrusters are designed with limited control capability. The infinite norm of vector is applied to describe the thrust saturation constraint, and this is given as follows [24]:
(24)U∞≤umax,where umax denotes the maximum thrust output.The relative velocity constraint is also considered to ensure that the servicer exhibits the capability to adjust the thrust output and avoid collision in case of emergency, and this is extremely important in actual on-orbit approaches. The velocity constraint is expressed as follows [27]:
(25)ẋ≤ẋmax,ẏ≤ẏmax,ż≤żmax,where ẋmax, ẏmax, and żmax denote the components of the maximum approach velocity in each direction.


## 2.3.1. Relative Translational Dynamics

It is assumed that the target is on a circular orbit and the distance between the servicer and target is extremely low. Given a short-period approach, the disturbances caused by solar pressure, atmospheric drag, and the effects of nonspherical gravity perturbation are neglected in both orbital and attitude motions. The Clohessy–Wiltshire (CW) equation is subsequently used to describe the relative translational dynamics [13]:
(1)x¨−2nẏ−3n2x=ux,y¨+2nẋ=uy,z¨+n2ż=uz,where ux,uy,uz denote the components of the control acceleration and are resolved into ∑t, n=μ/Rt3 denotes the orbital rate of the target with respect to ∑I, and μ denotes the gravitational constant.Eq. (1) is rewritten in the form of a state space model as follows:
(2)Ẋ=AX+BU,where X∈ℜ6 denotes the state vector, U∈ℜ3 denotes the control vector, and
(3)X=xyzẋẏż,A=0001000000100000013n20002n0000−2n0000−n2000,B=000000000100010001,U=uxuyuz.With respect to the implementation of the MPC, the discrete-time model of the relative translational dynamics is derived with a sampling periodTs as follows:
(4)Xk+1=AdXk+BdUk,where Xk and Uk denote the state and control vectors, respectively, at sampling instant k∈Z+. The matrices Ad∈ℜ6×6 and Bd∈ℜ6×3 are defined as follows [28]:
(5)Ad=eATs,Bd=∫0TseAτdτB.


## 2.3.2. Attitude Dynamics of the Target

The attitude dynamics of the target are modeled to describe the relationship between∑t and ∑b. The attitude of the target is parameterized by using the rotation quaternion q=q0,q¯TT, where the first component denotes the scalar part and the other is a three-dimensional vector. The rotation matrix associated with q, denoted by Rtbq, is given as follows [32]:
(6)Rtbq=q02−q¯Tq¯Ι3+2q¯q¯T−2q0q¯×,where Rtbq transforms a vector from ∑t into ∑b, Ι3 denotes a three-dimensional identity matrix, and ∗× denotes a skew-symmetric matrix that represents the cross product operator.The angular velocityw of ∑b relative to ∑t is defined as follows:
(7)w=wb−wt,where wb and wt denote the angular velocities of ∑b and ∑t, respectively, relative to ∑I.The first time derivative of Eq. (7) in the ∑I is given as follows:
(8)dwdtI=dwbdtI−dwtdtI.Based on Coriolis’s theorem, the following expression is obtained:(9)dwdtI=dwdtt+wt×w.Given thatdwb/dt∣I=dwb/dt∣b and dwt/dt∣I=dwt/dt∣t, a combination of Eqs. (8) and (9) yields the following expression:
(10)dwbdtb=dwdtt+dwtdtt+wt×w.The attitude dynamics of the target are given by the following expression:(11)Itdwbdtb+wb×Itwb=Nt,where It and Nt denote the inertial tensors of the target and external torques, respectively. It is assumed that the target is on a circular orbit and environmental perturbations are neglected, and dwt/dt∣t=0 and Nt=0 are subsequently obtained, respectively.Eq. (10) is substituted into Eq. (11) to yield the following expression:
(12)Itẇ+wt×w+wt+w×It⋅wt+w=0,where wt, resolved into ∑t, is given as wt=0,0,nT.Subsequently, the change in target’s attitude is given as follows [32]:
(13)q̇=12q0Ι3+q¯×−q¯×w.According to Eq. (6), the rotation matrix Rtbk+1 at the sampling instant k+1 is a function of rotation quaternion qk+1. According to Eqs. (12) and (13), qk+1 can be obtained through solving the numerical integration with the given wk and qk. Thereby, Rtbk+1 can be represented as a nonlinear function expression:
(14)Rtbk+1=fwk,qk.The fourth-order Runge–Kutta algorithm is utilized to solve the numerical integration in Eqs. (12) and (13). It is assumed that there is no attitude maneuvering of the failed target and that there is no external torque applied to the target from the servicer. The offline calculation can be adopted to predict target’s attitude, which will be used in MPC to predict the relative translational states and servicer’s control input.


## 2.3.3. Control Objective

The study focuses on two types of the on-orbit failed spacecraft. The first includes two large solar panels that are symmetrically mounted on the body as shown in Figure2. The second includes solar panels on only one side where deployment fails on the other side as shown in Figure 3.Figure 2
Failed spacecraft with two mounted solar panels.Figure 3
Failed spacecraft with solar panels deployed on only one side.Thexb-axis is along the center line of the solar panels, the yb-axis is along the symmetric axis of target’s body and is orthogonal to the xb-axis, and the zb-axis completes the dextral triad. In the study, the target is assumed as tumbling freely. It is also assumed that target’s DP is located along the yb-axis and the position vector of the DP is rDPb=0,−l,0T resolved in ∑b. The tracking error is defined as δrt=Rtbtrt−rDPb, where rt denotes servicer’s trajectories resolved in ∑t. Therefore, the objective of the approach procedure involves driving the servicer while tracking with the target DP, and this is expressed as follows:
(15)limt→∞δrt≤ε,where ε is a constant denoting the upper bound of tracking error.


## 2.3.4. Constraint Modeling

With respect to the problem of the ultraclose approach, collision avoidance is the basic requirement. Given the aforementioned two types of the failed spacecraft, a double-ellipsoid composite envelope strategy is introduced to model target’s keep-out zone in which the spacecraft body and the solar panels are modeled as two different ellipsoids based on their sizes. Furthermore, a sphere envelope with radiusre is simplified as servicer’s keep-out zone.Intuitively, collision avoidance is attained when the envelopes of the target and the servicer do not intersect. In order to determine collision avoidance, expanded ellipsoids considering the sizes of envelopes of the target and the servicer are introduced, as shown in Figure4. The original ellipsoids denote target’s double-ellipsoid envelope. It is assumed that centers of the original ellipsoids for the first type of the failed spacecraft are located at its mass center. With respect to the original ellipsoid enveloping solar panels of the second type of the failed spacecraft, its center is located at the center of the solar panels.Figure 4
The double-ellipsoid composite envelope strategy for (a) the first type of a failed spacecraft and (b) the second type of a failed spacecraft.
(a)
(b)Furthermore, the scheme for discriminating collision avoidance between the original ellipsoid and the sphere is shown in Figure5.Figure 5
Scheme for discriminating collision avoidance between the original ellipsoid and the sphere.Without loss of generality, the original ellipsoid is as follows:(16)x−d2a2+y2b2+z2c2=1,where a, b, and c are constants designed to envelop the target body or the solar panels based on their sizes and a>b>c is satisfied. If the center of ellipsoid is located at target’s mass center, d=0 is satisfied. The parametric equation of the original ellipsoid is given as follows [17]:
(17)x=d+asinθcosφ,y=bsinθsinφ,z=ccosθ,where θ∈0,π and φ∈0,2π.The expanded ellipsoid envelope is as follows:(18)x−d2a+are/c2+y2b+bre/c2+z2c+re2=1.It should be noted that the value of the minimum half-axis of the expanded ellipsoid is the summation of the radius of the sphere and the value of the minimum half-axis of the original ellipsoid. The other two half-axes of expanded ellipsoid are augmented based on the product of the radius and the ratio between the corresponding half-axis and minimum half-axis.The parametric equation of the expanded ellipsoid is given as follows:(19)x=d+a+arecsinθcosφ,y=b+brecsinθsinφ,z=c+recosθ.Proposition 1.
Given the original ellipsoid in Eq. (16) and the sphere with radius re, the sufficient condition for collision avoidance between the original ellipsoid and the sphere is defined as follows:
(20)x−d2a+are/c2+y2b+bre/c2+z2c+re2−1>0,where x,y,z denote the coordinates of the sphere center. Evidently, Eq. (20) indicates that the center of the sphere is outside the expanded ellipsoid shown in Eq. (18).Proof.
Given a pointp0 using spherical coordinates θ,φ describing its direction, there are two other corresponding points located at the original ellipsoidal surface and the expanded ellipsoidal surface (points p1 and p2, as shown in Figure 5), and their coordinates are given as follows:
(21)p1:x1=d+asinθcosφ,y1=bsinθsinφ,z1=ccosθ,p2:x2=d+a+arecsinθcosφ,y2=b+brecsinθsinφ,z2=c+recosθ.
Subsequently, the tangent planes of the two points (planesL1 and L2, as shown in Figure 5) are given as follows:
(22)L1:sinθcosφax+sinθsinφby+cosθcz−sinθcosφad−1=0,L2:sinθcosφax+sinθsinφby+cosθcz−sinθcosφad−1−rec=0.
According to Eq. (22), it is noted that the two tangent planes are parallel. Subsequently, the minimum distance from the point p2 to the tangent plane L1 is equal to the distance between tangent planes L2 and L1, and this is given as follows:
(23)rc=rec2/a2sin2θcos2φ+c2/b2sin2θsin2φ+cos2θ.
Given thata>b>c is defined, rc≥re is subsequently obtained and it denotes that the minimum distance between the original ellipsoid and the expanded ellipsoid exceeds or is equal to the radius of the sphere. Thus, when the center of the sphere lies outside the expanded ellipsoid, the collision avoidance of the sphere and the original ellipsoid is satisfied. The proof is completed.Remark 1.
Given the double-ellipsoid envelope modeling target’s keep-out zone, the collision avoidance condition for the target and servicer is the collision avoidance between the sphere and the two original ellipsoids. The modeled ellipsoid envelope is fixed in target’s body frame. Therefore, the mathematical description of collision avoidance between the servicer and the target can be conducted in target’s body frame. Although this is a sufficient but not a necessary condition to determine collision avoidance, it is derived by using the form of concave constraint, which is conducive to the implementation of MPC as subsequently discussed.In practice, all thrusters are designed with limited control capability. The infinite norm of vector is applied to describe the thrust saturation constraint, and this is given as follows [24]:
(24)U∞≤umax,where umax denotes the maximum thrust output.The relative velocity constraint is also considered to ensure that the servicer exhibits the capability to adjust the thrust output and avoid collision in case of emergency, and this is extremely important in actual on-orbit approaches. The velocity constraint is expressed as follows [27]:
(25)ẋ≤ẋmax,ẏ≤ẏmax,ż≤żmax,where ẋmax, ẏmax, and żmax denote the components of the maximum approach velocity in each direction.


## 3. Model Predictive Control Formulation

By using the dynamics and constraints discussed in Section2, the model predictive controller is designed to approach the tumbling failed spacecraft based on the basic principle of finite predictive control. Given the nonconvex constraint of the collision avoidance condition, the convexification method is discussed and the sequential convex programming is subsequently used to implement the MPC.


### 3.1. Prediction of the State Variables

The predicted relative translational state sequence generated by Eq. (4) with state Xk and control input Uk is expressed as follows:
(26)Xsk=AsXk+BsUsk.Subsequently,(27)Xsk=Xk+1kXk+2kXk+3k⋮Xk+Npk,Usk=UkkUk+1kUk+2k⋮Uk+Nc−1k,As=AdAd2Ad3⋮AdNp,Bs=Bd00⋯0AdBdBd0⋯0Ad2BdAdBdBd⋯0⋮⋮⋮⋱0AdNp−1BdAdNp−2BdAdNp−3Bd⋯AdNp−NcBd.where Xk+i∣k denotes the predicted state variable at k+i with the given information Xk, Uk+i∣k denotes the predicted control input at k+i, Uk∣k=Uk, Np denotes the prediction horizon, and Nc denotes the control horizon [12].As described in Section2.2, the predicated rotation matrix Rtbk+i at k+i can be solved offline with the given initial w0 and q0. Thus, while implementing the MPC, the rotation matrix Rtbk+i is considered as known.


### 3.2. Optimization Index

The objective of the MPC involves minimizing the tracking error between the predicted states and the desired trajectory. In order to optimize the thruster fuel, the control effort is included in the objective function. Thus, based on Eq. (15), the objective function is defined as follows:
(28)Γk=∑i=1NpRtbk+irk+ik−rDPbTQiRtbk+irk+ik−rDPb+∑i=0Nc−1UTk+ikRiUk+ik,where Qi∈ℜ3×3 denotes a positive-definite state-weighting matrix, Ri∈ℜ3×3 denotes a positive-definite control-weighting matrix, and rk+i∣k is available from the predicted state sequence Xk+i∣k.
(29)rk+ik=CdXk+ik,Cd=100000010000001000.In order to simplify the expression in Eq. (28), matrix Cs=INp⊗Cd is introduced, where ⊗ denotes the Kronecker product of two matrices [27], and INp denotes a Np-dimensional identity matrix. If we define rd=rDPbT,rDPbT,…,rDPbTT, R~tbk=diagRtbk+1,Rtbk+2,…,Rtbk+Np, Q=diagQ1,Q2,…,QNp, and R=diagR1,R2,…,RNc, the objective function in Eq. (28) is expressed as follows:
(30)Γk=UskTWkUsk+PkUsk+MkTQMk,where Wk=R~tbkCsBsTQR~tbkCsBs+R, Pk=2MkTQR~tbkCsBs, and Mk=R~tbkCsAsXsk−rd are satisfied.


### 3.3. Inequality Constraints

In this section, the constraints discussed in Section2 are reformulated for the implementation of the MPC. We reconsider the collision avoidance constraint shown in Eq. (20), and the collision avoidance between the servicer and the target is given as follows:
(31)gmxb,yb,zb=xb−dm2am+amre/cm2+yb2bm+bmre/cm2+zb2cm+re2−1>0,m=1,2,where xb,yb,zb denote the components of the vector r resolved in ∑b and subscript m denotes the two different ellipsoids that model the keep-out zone of the target. Function gm is defined as the collision threshold, and the collision avoidance is attained when gm>0. Based on Eq. (31), the collision avoidance constraint for the MPC is expressed as follows:
(32)gmiUsk=Rtbk+iHiCsAsXk+BsUsk−dm,0,0TT×Lm×Rtbk+iHiCsAsXk+BsUsk−dm,0,0T−1>0,i=1,2,…,Np,m=1,2,where Lm=diagam+amre/cm2,bm+bmre/cm2,cm+re2 and Hi is expressed as follows:
(33)Hi=T1,T2,…,Tj,…,TNpTj=I3ifj=i,Tj=E3ifj≠i,where I3 denotes a three-dimensional identity matrix and E3 denotes a three-dimensional zero matrix.After a few manipulations (see AppendixA), Eq. (32) is rewritten as follows:
(34)gmiUsk=UskTF~miUsk+G~miUsk+GmiTLmGmi−1>0,i=1,2,…,Np,m=1,2,where F~mi=FiTLmFi, G~mi=2GmiTLmFi, Fi=Rtbk+iHiCsBs, and Gmi=Rtbk+iHiCsAsXk−dm,0,0T are satisfied.Based on Eq. (24), the control input constraint is reexpressed as follows:
(35)I3Nc;−I3NcUs≤Umax,where I3Nc denotes a 3Nc-dimensional identity matrix and Umax∈ℜ6Nc is defined as Umax=umax,⋯,umaxT.The velocity constraint expressed in Eq. (25) is reexpressed as follows (see Appendix B):
(36)I3Np;−I3NpDsAsXk+BsUs≤Vmax,where the matrices Ds and Vmax∈ℜ6Np are defined as follows:
(37)Dd=000100000010000001,Ds=INp⊗Dd,vmax=ẋmaxẏmaxżmaxẋmaxẏmaxżmaxT,Vmax=vmaxT,⋯,vmaxTT.


### 3.4. MPC by Using Sequential Convex Programming

The design of the MPC for the approach problem is summarized as follows:Problem 1.
(38)minUskΓk=UskTWkUsk+PkUsk+MkTQMksubjecttoI3Nc;−I3NcUsk≤UmaxI3Np;−I3NpDsAsXk+BsUsk≤VmaxgmiUsk=UskTF~miUsk+G~miUsk+GmiTLmGmi−1>0,i=1,2,…,Np,m=1,2>0.Given that the collision avoidance constraint is a nonconvex constraint, it is difficult to satisfy the convergence and the optimality of the solution to Problem1. In order to solve Problem 1, it is subsequently converted into a convex optimization problem and sequential convex programming is utilized to implement MPC. Given the collision avoidance constraint in Eq. (38), the linearization technique is adopted to convert it into linear terms as follows:
(39)gmiUsk≈gmiUsnk+∇gmiUsnkTUsk−Usnk>0,i=1,2,…,Np,m=1,2,where Usnk denotes the n-iterated solution of Usk, and ∇gmiUsnk=2F~miUsnk+G~miT denotes the gradient of gmi at Usnk.Proposition 2.
Eq. (39) is utilized as the collision avoidance constraint as opposed to Eq. (34), and Problem 1 is subsequently converted to Problem 2 as follows:Problem 2.
(40)minUskΓk=UskTWkUsk+PkUsk+MkTQMksubjecttoI3Nc;−I3NcUsk≤UmaxI3Np;−I3NpGsAsXk+BsUsk≤VmaxgmiUsnk+∇gmiUsnkTUsk−Usnk>0,i=1,2,…,Np,m=1,2.
Subsequently, anyUsk that is feasible for Problem 2 is also feasible for Problem 1.Proof.
LetUsk denote an arbitrary feasible point for Problem 2, i.e.,
(41)gmiUsnk+∇gmiUsnkTUsk−Usnk>0.
RepresentgmiUsk by the second-order Taylor series expansion at Usnk:
(42)gmiUsk=gmiUsnk+∇gmiUsnkTUsk−Usnk+12Usk−UsnkT∇2gmiUsnkUsk−Usnk,where ∇2gmi denotes the Hessian matrix of gmi. It should be noted that ∇2gmi=2F~mi denotes a constant positive-definite matrix where F~mi=FiTLmFi is satisfied. We substitute Eq. (41) into Eq. (42), and this results in gmiUsk>0, and thus, any Usk that satisfies the condition in Eq. (39) also satisfies the condition in Eq. (34). The proof is completed.It should be noted that Problem2 is a quadratic programming problem [24]. In order to optimize the control input Usk, sequential quadratic programming is utilized and a trust region between Usnk and Usk is introduced to ensure the convergence.
(43)Usk−Usnk∞≤ρnD0,where D0∈ℜ3Nc denotes the radius of the trust region and ρ∈0,1 denotes a parameter that determines the rate of convergence. After a few manipulations (see Appendix C), Problem 2 that considers the trust region is subsequently converted into Problem 3.Problem 3.
(44)minUskΓk=UskTWkUsk+PkUsk+MkTQMksubjecttoI3Nc;−I3NcUsk≤UmaxI3Np;−I3NpGsBsUsk≤Vmax−I3Np;−I3NpGsAsXk−∇gmiUsnkTUsk<−∇gmiUsnkTUsnk+gmiUsnk,i=1,2,…,Np,m=1,2I3Nc;−I3NcUsk≤ρnD0+Usnk;ρnD0−Usnk.Remark 2.
Based on Proposition2, if there exists a solution feasible for Problem 3, the solution is also feasible for Problem 2. Furthermore, an optimized solution is attained based on sequential quadratic programming for Problem 3. The successive solution process based on Problem 3 requires a first solution Us1k to commence. In the study, Us1k is set to a zero vector and a maximum iteration number n is designed to end the sequential quadratic programming.Thus, the constrained optimization problem in the context of the MPC is converted into a sequential quadratic programming problem. The standard quadratic programming algorithm can be employed to solve this problem easily. With the given state vectorXk, rotation quaternion qk, and relative angular velocity wk at k, the MPC should be solved at each sampling instant k by using sequential quadratic programming to obtain a sequence of the control input where only the first control input is applied.


## 3.1. Prediction of the State Variables

The predicted relative translational state sequence generated by Eq. (4) with state Xk and control input Uk is expressed as follows:
(26)Xsk=AsXk+BsUsk.Subsequently,(27)Xsk=Xk+1kXk+2kXk+3k⋮Xk+Npk,Usk=UkkUk+1kUk+2k⋮Uk+Nc−1k,As=AdAd2Ad3⋮AdNp,Bs=Bd00⋯0AdBdBd0⋯0Ad2BdAdBdBd⋯0⋮⋮⋮⋱0AdNp−1BdAdNp−2BdAdNp−3Bd⋯AdNp−NcBd.where Xk+i∣k denotes the predicted state variable at k+i with the given information Xk, Uk+i∣k denotes the predicted control input at k+i, Uk∣k=Uk, Np denotes the prediction horizon, and Nc denotes the control horizon [12].As described in Section2.2, the predicated rotation matrix Rtbk+i at k+i can be solved offline with the given initial w0 and q0. Thus, while implementing the MPC, the rotation matrix Rtbk+i is considered as known.


## 3.2. Optimization Index

The objective of the MPC involves minimizing the tracking error between the predicted states and the desired trajectory. In order to optimize the thruster fuel, the control effort is included in the objective function. Thus, based on Eq. (15), the objective function is defined as follows:
(28)Γk=∑i=1NpRtbk+irk+ik−rDPbTQiRtbk+irk+ik−rDPb+∑i=0Nc−1UTk+ikRiUk+ik,where Qi∈ℜ3×3 denotes a positive-definite state-weighting matrix, Ri∈ℜ3×3 denotes a positive-definite control-weighting matrix, and rk+i∣k is available from the predicted state sequence Xk+i∣k.
(29)rk+ik=CdXk+ik,Cd=100000010000001000.In order to simplify the expression in Eq. (28), matrix Cs=INp⊗Cd is introduced, where ⊗ denotes the Kronecker product of two matrices [27], and INp denotes a Np-dimensional identity matrix. If we define rd=rDPbT,rDPbT,…,rDPbTT, R~tbk=diagRtbk+1,Rtbk+2,…,Rtbk+Np, Q=diagQ1,Q2,…,QNp, and R=diagR1,R2,…,RNc, the objective function in Eq. (28) is expressed as follows:
(30)Γk=UskTWkUsk+PkUsk+MkTQMk,where Wk=R~tbkCsBsTQR~tbkCsBs+R, Pk=2MkTQR~tbkCsBs, and Mk=R~tbkCsAsXsk−rd are satisfied.


## 3.3. Inequality Constraints

In this section, the constraints discussed in Section2 are reformulated for the implementation of the MPC. We reconsider the collision avoidance constraint shown in Eq. (20), and the collision avoidance between the servicer and the target is given as follows:
(31)gmxb,yb,zb=xb−dm2am+amre/cm2+yb2bm+bmre/cm2+zb2cm+re2−1>0,m=1,2,where xb,yb,zb denote the components of the vector r resolved in ∑b and subscript m denotes the two different ellipsoids that model the keep-out zone of the target. Function gm is defined as the collision threshold, and the collision avoidance is attained when gm>0. Based on Eq. (31), the collision avoidance constraint for the MPC is expressed as follows:
(32)gmiUsk=Rtbk+iHiCsAsXk+BsUsk−dm,0,0TT×Lm×Rtbk+iHiCsAsXk+BsUsk−dm,0,0T−1>0,i=1,2,…,Np,m=1,2,where Lm=diagam+amre/cm2,bm+bmre/cm2,cm+re2 and Hi is expressed as follows:
(33)Hi=T1,T2,…,Tj,…,TNpTj=I3ifj=i,Tj=E3ifj≠i,where I3 denotes a three-dimensional identity matrix and E3 denotes a three-dimensional zero matrix.After a few manipulations (see AppendixA), Eq. (32) is rewritten as follows:
(34)gmiUsk=UskTF~miUsk+G~miUsk+GmiTLmGmi−1>0,i=1,2,…,Np,m=1,2,where F~mi=FiTLmFi, G~mi=2GmiTLmFi, Fi=Rtbk+iHiCsBs, and Gmi=Rtbk+iHiCsAsXk−dm,0,0T are satisfied.Based on Eq. (24), the control input constraint is reexpressed as follows:
(35)I3Nc;−I3NcUs≤Umax,where I3Nc denotes a 3Nc-dimensional identity matrix and Umax∈ℜ6Nc is defined as Umax=umax,⋯,umaxT.The velocity constraint expressed in Eq. (25) is reexpressed as follows (see Appendix B):
(36)I3Np;−I3NpDsAsXk+BsUs≤Vmax,where the matrices Ds and Vmax∈ℜ6Np are defined as follows:
(37)Dd=000100000010000001,Ds=INp⊗Dd,vmax=ẋmaxẏmaxżmaxẋmaxẏmaxżmaxT,Vmax=vmaxT,⋯,vmaxTT.


## 3.4. MPC by Using Sequential Convex Programming

The design of the MPC for the approach problem is summarized as follows:Problem 1.
(38)minUskΓk=UskTWkUsk+PkUsk+MkTQMksubjecttoI3Nc;−I3NcUsk≤UmaxI3Np;−I3NpDsAsXk+BsUsk≤VmaxgmiUsk=UskTF~miUsk+G~miUsk+GmiTLmGmi−1>0,i=1,2,…,Np,m=1,2>0.Given that the collision avoidance constraint is a nonconvex constraint, it is difficult to satisfy the convergence and the optimality of the solution to Problem1. In order to solve Problem 1, it is subsequently converted into a convex optimization problem and sequential convex programming is utilized to implement MPC. Given the collision avoidance constraint in Eq. (38), the linearization technique is adopted to convert it into linear terms as follows:
(39)gmiUsk≈gmiUsnk+∇gmiUsnkTUsk−Usnk>0,i=1,2,…,Np,m=1,2,where Usnk denotes the n-iterated solution of Usk, and ∇gmiUsnk=2F~miUsnk+G~miT denotes the gradient of gmi at Usnk.Proposition 2.
Eq. (39) is utilized as the collision avoidance constraint as opposed to Eq. (34), and Problem 1 is subsequently converted to Problem 2 as follows:Problem 2.
(40)minUskΓk=UskTWkUsk+PkUsk+MkTQMksubjecttoI3Nc;−I3NcUsk≤UmaxI3Np;−I3NpGsAsXk+BsUsk≤VmaxgmiUsnk+∇gmiUsnkTUsk−Usnk>0,i=1,2,…,Np,m=1,2.
Subsequently, anyUsk that is feasible for Problem 2 is also feasible for Problem 1.Proof.
LetUsk denote an arbitrary feasible point for Problem 2, i.e.,
(41)gmiUsnk+∇gmiUsnkTUsk−Usnk>0.
RepresentgmiUsk by the second-order Taylor series expansion at Usnk:
(42)gmiUsk=gmiUsnk+∇gmiUsnkTUsk−Usnk+12Usk−UsnkT∇2gmiUsnkUsk−Usnk,where ∇2gmi denotes the Hessian matrix of gmi. It should be noted that ∇2gmi=2F~mi denotes a constant positive-definite matrix where F~mi=FiTLmFi is satisfied. We substitute Eq. (41) into Eq. (42), and this results in gmiUsk>0, and thus, any Usk that satisfies the condition in Eq. (39) also satisfies the condition in Eq. (34). The proof is completed.It should be noted that Problem2 is a quadratic programming problem [24]. In order to optimize the control input Usk, sequential quadratic programming is utilized and a trust region between Usnk and Usk is introduced to ensure the convergence.
(43)Usk−Usnk∞≤ρnD0,where D0∈ℜ3Nc denotes the radius of the trust region and ρ∈0,1 denotes a parameter that determines the rate of convergence. After a few manipulations (see Appendix C), Problem 2 that considers the trust region is subsequently converted into Problem 3.Problem 3.
(44)minUskΓk=UskTWkUsk+PkUsk+MkTQMksubjecttoI3Nc;−I3NcUsk≤UmaxI3Np;−I3NpGsBsUsk≤Vmax−I3Np;−I3NpGsAsXk−∇gmiUsnkTUsk<−∇gmiUsnkTUsnk+gmiUsnk,i=1,2,…,Np,m=1,2I3Nc;−I3NcUsk≤ρnD0+Usnk;ρnD0−Usnk.Remark 2.
Based on Proposition2, if there exists a solution feasible for Problem 3, the solution is also feasible for Problem 2. Furthermore, an optimized solution is attained based on sequential quadratic programming for Problem 3. The successive solution process based on Problem 3 requires a first solution Us1k to commence. In the study, Us1k is set to a zero vector and a maximum iteration number n is designed to end the sequential quadratic programming.Thus, the constrained optimization problem in the context of the MPC is converted into a sequential quadratic programming problem. The standard quadratic programming algorithm can be employed to solve this problem easily. With the given state vectorXk, rotation quaternion qk, and relative angular velocity wk at k, the MPC should be solved at each sampling instant k by using sequential quadratic programming to obtain a sequence of the control input where only the first control input is applied.


## 4. Simulation Results and Discussion

In order to test the performance of the proposed MPC method, numerical simulations are described in this section. It is assumed that the target is on GEO, and the simulation parameters relevant to the target are listed in Table1.Table 1
Parameters of the failed spacecraft.
Parameter
Value
Orbit parameter
h=35786km, e=0, Ω=0∘, i=0∘, w=0∘, f=0∘
Inertial moment
It=diag500,600,700kg⋅m2
Docking port
rDPb=0,−3,0Tm
Initial rotation quaternion
q0=1,0,0,0TConsidering the orbit parameters of the target described in Table1, h is the orbit height, e is the eccentricity, Ω is the right ascension of ascending node, w is the argument of perigee, and f is the true anomaly. Given the two failed spacecraft, both their body’s keep-out zones are modeled as an ellipsoid envelope with a1=2m, b1=1.8m, c1=1.5m, and d1=0m. The parameters of the ellipsoid envelope for the first failed spacecraft’s solar panels are set to a2=6m, b2=1.2m, c1=1.2m, and d2=0m. The parameters of the ellipsoid envelope for the second failed spacecraft’s solar panels are set to a2=2m, b2=1.2m, c1=1.2m, and d2=4m. The radius of the sphere envelope for the servicer spacecraft is re=0.5m.In our simulations, the sampling periodTs is set to Ts=0.1s. The weighting matrices are selected as Qi=I3 and Ri=I3. The maximum control input umax is set to umax=0.1m/s2, and the maximum relative velocity vmax is set to vmax=0.5,0.5,0.5,0.5,0.5,0.5Tm/s. The value of the prediction horizon Np is set to Np=20, and the value of control horizon Nc is set to Nc=10. The value of maximum iteration number n is set to n=5, and the trust region related parameters are set to ρ=0.9 and D0=0.04,0.04,⋯,0.04,⋯,0.04T.


### 4.1. Approaching the First Type of a Failed Spacecraft

First, the effectiveness of the proposed MPC is analyzed through the approach simulation. The initial relative translational statesX0 are set to 25,−25.3,0Tm and 0,−0.0027,0Tm/s. The initial angular velocities of the failed spacecraft are set to wbb=2.2,2.5,2.3Tdeg/s. The simulation time is set to 100 s. The simulation results are shown in Figures 6 and 7.Figure 6
Temporal response of the relative position, relative velocity, control input, and collision threshold. (a) Temporal history along thex direction. (b) Temporal history along the y direction. (c) Temporal history along the z direction. (d) Temporal history of relative velocity. (e) Temporal history of the control input. (f) Temporal history of the collision threshold.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Three-dimensional plot of the approaching trajectories resolved in the target body-fixed frame. (a) Servicer’s center relative to target’s expanded ellipsoids. (b) Servicer’s sphere relative to target’s original ellipsoids.
(a)
(b)As shown in Figures6(a)–6(c), it is evident that the desired trajectory to track the DP is time variant in target’s LVLH frame, the approach process is implemented, and the trajectory tracking error in each direction is at 10−4 m magnitude. As depicted in Figure 6(d), the relative velocity in each direction is lower than 0.5 m/s, and the relative velocity constraint is satisfied. As shown in Figure 6(e), the control input in each direction is lower than 0.1 m/s2 and the control input constraint is satisfied. As depicted in Figure 6(f), while approaching the target, the collision threshold always exceeds zero, thereby indicating that the center of the servicer lies outside the expanded two ellipsoids and that the collision constraint between the target and servicer is attained. Furthermore, the three-dimensional trajectories of servicer’s center relative to the expanded two ellipsoids of the target and servicer’s keep-out zone relative to target’s keep-out zone are shown in Figure 7. Evidently, the collision constraint is satisfied.Subsequently, the performances of MPC on target with different initial angular velocities are analyzed. In addition to the initial angular velocitieswbb=2.2,2.5,2.3Tdeg/s (case A) as discussed in the previous simulation, initial angular velocities wbb=1.2,1.0,0.8Tdeg/s (case B) and wbb=0,0,0Tdeg/s (case C) are analyzed for comparison purposes and all other simulation parameters are set as the same. Furthermore, the convergence time T and total control input ΔV are calculated at different initial angular velocities. It should be noted that the convergence time T is recorded when the distance between the servicer and the DP is less than 0.1 m. The analysis results are shown in Figure 8.Figure 8
Three-dimensional plot of the transfer trajectories at target’s different initial angular velocities resolved in the target body-fixed frame.The servicer can arrive at the DP while satisfying above control input, relative velocity, and collision avoidance constraints in the aforementioned three cases. As shown in Figure8, the angular velocity of the target significantly influences servicer’s transfer trajectory. When the target exhibits low angular velocity, the total control input of the approaching is lower. Furthermore, the influence of the control horizon and predictive horizon on the MPC is also analyzed. With respect to each control horizon Nc=5,10,20, the convergence time and total control input are calculated at different control horizons Np=20,30,40,50,60. The initial angular velocities of the target are also set to wbb=2.2,2.5,2.3Tdeg/s. The simulation results are shown in Figure 9.Figure 9
Influences of the predictive horizonNp and control horizon Nc on (a) convergence time and (b) total control input.
(a)
(b)As shown in Figure9, with increases in the predictive horizon, the convergence time for approaching increases and the total control input decreases. It is noted that the objective function in Eq. (28) is consist of the tracking error and the control effort. The increase of predictive horizon Np denotes more predicted states, leading to less state error and less control input to modify the state error. Accordingly, the required total control input is less and the convergence time will be longer. With increases in the control horizon Nc, the convergence time for approaching decreases, while the total control input and the accuracy of tracking the DP increase. This is because more predicted control inputs have been used to modify the state error in every sampling instant, leading to a shorter convergence time and higher tracking accuracy. Furthermore, the increases in the predictive horizon and control horizon, denoting more predicted states and control inputs, lead to a higher computation load when solving the underlying optimization problem in MPC. Therefore, while implementing the on-orbit missions, proper MPC parameters should be selected based on mission requirements.


### 4.2. Approaching the Second Type of a Failed Spacecraft

With respect to the simulation that involves approaching the second type of the failed spacecraft, the parameters of the MPC are set as identical to those in the previous simulation. The parameters of the second failed spacecraft’s keep-out zone are discussed at the beginning of Section4. Furthermore, in order to confirm that the proposed method is generally applicable, the designed MPC dealing with the servicer at different initial positions, such as r0=25,‐25.3,0Tm (Sat 1), r0=25,25.3,0Tm (Sat 2), r0=−25,25.3,0Tm (Sat 3), and r0=−25,−25.3,0Tm (Sat 4), is analyzed. The initial angular velocities of the target are also set to wbb=2.2,2.5,2.3Tdeg/s. Simulation results on the transfer trajectory are shown in Figure 10.Figure 10
Three-dimensional plot of the transfer trajectories with the servicer at different initial positions.As shown in Figure10, the proposed MPC exhibits the ability to drive the servicers at different positions to the DP of the target. The convergence times for each servicer are 57.9 s, 59 s, 69.2 s, and 69.2 s. The total control inputs for each servicer are 3.69 m/s, 4.23 m/s, 5.46 m/s, and 4.83 m/s. The temporal response of the collision threshold of each servicer is shown in Figure 11.Figure 11
Temporal response of the collision threshold of each servicer.As shown in Figure11, the collision avoidance constraint is evidently satisfied during the approaching process. As shown in Figure 12, required constraints, such as control input and relative velocity, are also satisfied. It should be noted that when compared with linearizing the collision avoidance constraint at the DP in which only half of the plane is the collision avoidance zone, the proposed method for handling collision avoidance constraint exhibits a better performance in which the space outside the expanded ellipsoids is the collision avoidance zone.Figure 12
Temporal response of the control input (left) and relative velocities (right) for (a) Sat 1, (b) Sat 2, (c) Sat 3, and (d) Sat 4.
(a)
(b)
(c)
(d)In previous simulation analysis, the controller uncertainties are neglected and the measurement of relative states and the control input of the system are assumed to be completely accurate. In order to evaluate the performance of the proposed MPC for real approaching missions, navigation and control noise are considered. Considering unknown but bounded navigation noise, the measurement-relative translational statesX^k can be represented as
(45)X^k=Xk+ek,X^k+1=Xk+1+ek+1,where ek=epk;evk and epk and evk are the measurement errors of position and velocity, respectively. Considering that the control output noise was proportional to the control output [33], the control input U^k applied to the servicer can be represented as
(46)U^k=I3+ΛkUk,where Λk=diagΛxk,Λyk,Λzk is a three-dimensional diagonal matrix.Combining Eqs. (4), (45), and (46), the prediction of relative states should be represented as
(47)X^k+1=AdX^k+BdI+ΛkUk+ek+1−Adek.However, navigation and control noise cannot be measured. Then, Eq. (4) with measurement X^k and control input Uk was used to predict relative states. The parameters of the MPC are set as identical to those in the previous simulation, and the initial position of the servicer is r0=25,‐25.3,0Tm (Sat 1). The initial angular velocities of the target are also set to wbb=2.2,2.5,2.3Tdeg/s.A Monte-Carlo analysis is performed to illustrate the influence of uncertainties on the MPC performance. It is assumed that elements inep, ev, and Λ follow Gaussian distribution, and their accuracies (3σ) are defined as Ep, Ev, and Eu, respectively. When developing Monte-Carlo analysis, 1000 simulations in each of two conditions such as Ep,Ev,Eu=0.005m,10−4m/s,0.01 and Ep,Ev,Eu=0.02m,10−3m/s,0.02 are conducted. Reconsidering the tracking error defined in Eq. (15), the average tracking error δr~ is considered to evaluate the influence of uncertainties, which is defined as
(48)δr~=∫tstfδrtdttf−ts,where tf=100s is the end time of simulation and ts=80s is designed to start recording the tracking error. The Monte-Carlo simulation results regarding the average tracking error δr~ are shown in Figure 13.Figure 13
Performance comparison in terms of controller uncertainties.When developing Monte-Carlo analysis, the required constraints such as control input, relative velocity, and collision avoidance in every simulation are satisfied. When there are no navigation and control noise, the average tracking error is approximately 0.5 mm. As shown in Figure13, when the uncurtains are Ep,Ev,Eu=0.005m,10−4m/s,0.01 and Ep,Ev,Eu=0.02m,10−3m/s,0.02, the average tracking errors are approximately 1.5 mm and 7 mm, respectively. With the increase of uncurtains, the average tracking error also becomes larger. The effectiveness of the MPC in handling navigation and control noise is also proved. Based on the above analyses, it is concluded that the proposed MPC strategy can drive the servicer to the DP of a tumbling target while satisfying various constraints and providing capability in handling navigation and control noises.


## 4.1. Approaching the First Type of a Failed Spacecraft

First, the effectiveness of the proposed MPC is analyzed through the approach simulation. The initial relative translational statesX0 are set to 25,−25.3,0Tm and 0,−0.0027,0Tm/s. The initial angular velocities of the failed spacecraft are set to wbb=2.2,2.5,2.3Tdeg/s. The simulation time is set to 100 s. The simulation results are shown in Figures 6 and 7.Figure 6
Temporal response of the relative position, relative velocity, control input, and collision threshold. (a) Temporal history along thex direction. (b) Temporal history along the y direction. (c) Temporal history along the z direction. (d) Temporal history of relative velocity. (e) Temporal history of the control input. (f) Temporal history of the collision threshold.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Three-dimensional plot of the approaching trajectories resolved in the target body-fixed frame. (a) Servicer’s center relative to target’s expanded ellipsoids. (b) Servicer’s sphere relative to target’s original ellipsoids.
(a)
(b)As shown in Figures6(a)–6(c), it is evident that the desired trajectory to track the DP is time variant in target’s LVLH frame, the approach process is implemented, and the trajectory tracking error in each direction is at 10−4 m magnitude. As depicted in Figure 6(d), the relative velocity in each direction is lower than 0.5 m/s, and the relative velocity constraint is satisfied. As shown in Figure 6(e), the control input in each direction is lower than 0.1 m/s2 and the control input constraint is satisfied. As depicted in Figure 6(f), while approaching the target, the collision threshold always exceeds zero, thereby indicating that the center of the servicer lies outside the expanded two ellipsoids and that the collision constraint between the target and servicer is attained. Furthermore, the three-dimensional trajectories of servicer’s center relative to the expanded two ellipsoids of the target and servicer’s keep-out zone relative to target’s keep-out zone are shown in Figure 7. Evidently, the collision constraint is satisfied.Subsequently, the performances of MPC on target with different initial angular velocities are analyzed. In addition to the initial angular velocitieswbb=2.2,2.5,2.3Tdeg/s (case A) as discussed in the previous simulation, initial angular velocities wbb=1.2,1.0,0.8Tdeg/s (case B) and wbb=0,0,0Tdeg/s (case C) are analyzed for comparison purposes and all other simulation parameters are set as the same. Furthermore, the convergence time T and total control input ΔV are calculated at different initial angular velocities. It should be noted that the convergence time T is recorded when the distance between the servicer and the DP is less than 0.1 m. The analysis results are shown in Figure 8.Figure 8
Three-dimensional plot of the transfer trajectories at target’s different initial angular velocities resolved in the target body-fixed frame.The servicer can arrive at the DP while satisfying above control input, relative velocity, and collision avoidance constraints in the aforementioned three cases. As shown in Figure8, the angular velocity of the target significantly influences servicer’s transfer trajectory. When the target exhibits low angular velocity, the total control input of the approaching is lower. Furthermore, the influence of the control horizon and predictive horizon on the MPC is also analyzed. With respect to each control horizon Nc=5,10,20, the convergence time and total control input are calculated at different control horizons Np=20,30,40,50,60. The initial angular velocities of the target are also set to wbb=2.2,2.5,2.3Tdeg/s. The simulation results are shown in Figure 9.Figure 9
Influences of the predictive horizonNp and control horizon Nc on (a) convergence time and (b) total control input.
(a)
(b)As shown in Figure9, with increases in the predictive horizon, the convergence time for approaching increases and the total control input decreases. It is noted that the objective function in Eq. (28) is consist of the tracking error and the control effort. The increase of predictive horizon Np denotes more predicted states, leading to less state error and less control input to modify the state error. Accordingly, the required total control input is less and the convergence time will be longer. With increases in the control horizon Nc, the convergence time for approaching decreases, while the total control input and the accuracy of tracking the DP increase. This is because more predicted control inputs have been used to modify the state error in every sampling instant, leading to a shorter convergence time and higher tracking accuracy. Furthermore, the increases in the predictive horizon and control horizon, denoting more predicted states and control inputs, lead to a higher computation load when solving the underlying optimization problem in MPC. Therefore, while implementing the on-orbit missions, proper MPC parameters should be selected based on mission requirements.


## 4.2. Approaching the Second Type of a Failed Spacecraft

With respect to the simulation that involves approaching the second type of the failed spacecraft, the parameters of the MPC are set as identical to those in the previous simulation. The parameters of the second failed spacecraft’s keep-out zone are discussed at the beginning of Section4. Furthermore, in order to confirm that the proposed method is generally applicable, the designed MPC dealing with the servicer at different initial positions, such as r0=25,‐25.3,0Tm (Sat 1), r0=25,25.3,0Tm (Sat 2), r0=−25,25.3,0Tm (Sat 3), and r0=−25,−25.3,0Tm (Sat 4), is analyzed. The initial angular velocities of the target are also set to wbb=2.2,2.5,2.3Tdeg/s. Simulation results on the transfer trajectory are shown in Figure 10.Figure 10
Three-dimensional plot of the transfer trajectories with the servicer at different initial positions.As shown in Figure10, the proposed MPC exhibits the ability to drive the servicers at different positions to the DP of the target. The convergence times for each servicer are 57.9 s, 59 s, 69.2 s, and 69.2 s. The total control inputs for each servicer are 3.69 m/s, 4.23 m/s, 5.46 m/s, and 4.83 m/s. The temporal response of the collision threshold of each servicer is shown in Figure 11.Figure 11
Temporal response of the collision threshold of each servicer.As shown in Figure11, the collision avoidance constraint is evidently satisfied during the approaching process. As shown in Figure 12, required constraints, such as control input and relative velocity, are also satisfied. It should be noted that when compared with linearizing the collision avoidance constraint at the DP in which only half of the plane is the collision avoidance zone, the proposed method for handling collision avoidance constraint exhibits a better performance in which the space outside the expanded ellipsoids is the collision avoidance zone.Figure 12
Temporal response of the control input (left) and relative velocities (right) for (a) Sat 1, (b) Sat 2, (c) Sat 3, and (d) Sat 4.
(a)
(b)
(c)
(d)In previous simulation analysis, the controller uncertainties are neglected and the measurement of relative states and the control input of the system are assumed to be completely accurate. In order to evaluate the performance of the proposed MPC for real approaching missions, navigation and control noise are considered. Considering unknown but bounded navigation noise, the measurement-relative translational statesX^k can be represented as
(45)X^k=Xk+ek,X^k+1=Xk+1+ek+1,where ek=epk;evk and epk and evk are the measurement errors of position and velocity, respectively. Considering that the control output noise was proportional to the control output [33], the control input U^k applied to the servicer can be represented as
(46)U^k=I3+ΛkUk,where Λk=diagΛxk,Λyk,Λzk is a three-dimensional diagonal matrix.Combining Eqs. (4), (45), and (46), the prediction of relative states should be represented as
(47)X^k+1=AdX^k+BdI+ΛkUk+ek+1−Adek.However, navigation and control noise cannot be measured. Then, Eq. (4) with measurement X^k and control input Uk was used to predict relative states. The parameters of the MPC are set as identical to those in the previous simulation, and the initial position of the servicer is r0=25,‐25.3,0Tm (Sat 1). The initial angular velocities of the target are also set to wbb=2.2,2.5,2.3Tdeg/s.A Monte-Carlo analysis is performed to illustrate the influence of uncertainties on the MPC performance. It is assumed that elements inep, ev, and Λ follow Gaussian distribution, and their accuracies (3σ) are defined as Ep, Ev, and Eu, respectively. When developing Monte-Carlo analysis, 1000 simulations in each of two conditions such as Ep,Ev,Eu=0.005m,10−4m/s,0.01 and Ep,Ev,Eu=0.02m,10−3m/s,0.02 are conducted. Reconsidering the tracking error defined in Eq. (15), the average tracking error δr~ is considered to evaluate the influence of uncertainties, which is defined as
(48)δr~=∫tstfδrtdttf−ts,where tf=100s is the end time of simulation and ts=80s is designed to start recording the tracking error. The Monte-Carlo simulation results regarding the average tracking error δr~ are shown in Figure 13.Figure 13
Performance comparison in terms of controller uncertainties.When developing Monte-Carlo analysis, the required constraints such as control input, relative velocity, and collision avoidance in every simulation are satisfied. When there are no navigation and control noise, the average tracking error is approximately 0.5 mm. As shown in Figure13, when the uncurtains are Ep,Ev,Eu=0.005m,10−4m/s,0.01 and Ep,Ev,Eu=0.02m,10−3m/s,0.02, the average tracking errors are approximately 1.5 mm and 7 mm, respectively. With the increase of uncurtains, the average tracking error also becomes larger. The effectiveness of the MPC in handling navigation and control noise is also proved. Based on the above analyses, it is concluded that the proposed MPC strategy can drive the servicer to the DP of a tumbling target while satisfying various constraints and providing capability in handling navigation and control noises.


## 5. Conclusion

A model predictive control (MPC) method for a servicer spacecraft autonomously approaching a tumbling failed spacecraft is presented in the study. The objective of the proposed MPC involves driving the servicer tracking the time-variant motion of the docking port of a target by considering collision avoidance, control input saturation, and velocity constraints. The relative translation is predicted by the CW equation, and target’s attitude is predicted by a derived attitude dynamics. Sufficient conditions for collision avoidance are derived by using the form of concave constraint. The underlying optimization program for the implementation of MPC is converted into a convex optimization problem, and this is solved by using sequential convex programming. The process of the approach is simulated to evaluate the performance of the MPC strategy, and a few main contributions are obtained as follows:(1)
Given the complex configurations of two types of the failed spacecraft, a double-ellipsoid envelope is designed to model the keep-out zone. Sufficient conditions with simple mathematical expressions for collision avoidance are derived. The proposed modeling method for collision avoidance can be extended to other collision avoidance missions
(2)
The convexification method for the collision avoidance constraints and the sequential convex programming for solving the underlying optimization problem in MPC are effective. The proposed MPC exhibits the capacity to handle the problem of tracking the time-varying docking point while handling various constraints



---

*Source: 7428535-2018-12-06.xml*