Model Predictive Control to Autonomously Approach a Failed Spacecraft

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 keepout 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.


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][15][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][19][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][24][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][27][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][30][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: Section 2 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.

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 Figure 1.
As shown in Figure 1, r ∈ ℜ 3 denotes the vector directed from the center of mass (c.m.) of the target to the c.m. of the servicer.Additionally, R t and R s denote the distances of the target and servicer, respectively, relative to the center of the Earth, and R e denotes the radius of the Earth.

Inertial
Coordinate System ∑ I .The origin of the inertial coordinate system o I − x I y I z I is centered on Earth, the z I -axis is along the rotational axis, and the x I -axis points toward the vernal equinox.The y I -axis completes the right-handed orthogonal.

Target Local Vertical-Local Horizontal (LVLH) Coordinate
System ∑ t .The origin of the LVLH coordinate system o − x t y t z t lies on the c.m. of the target.The x t -axis is in the opposite direction to Earth's center, the y t -axis is along the direction of flight, and the z t -axis completes the right-handed orthogonal.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]: where u x , u y , u z denote the components of the control acceleration and are resolved into ∑ t , n = μ/R 3 t 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: where X ∈ ℜ 6 denotes the state vector, U ∈ ℜ 3 denotes the control vector, and With respect to the implementation of the MPC, the discrete-time model of the relative translational dynamics is derived with a sampling period T s as follows: where X k and U k denote the state and control vectors, respectively, at sampling instant k ∈ Z + .The matrices A d ∈ ℜ 6×6 and B d ∈ ℜ 6×3 are defined as follows [28]: 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 = q 0 , q T T , where the first component denotes the scalar part and the other is a threedimensional vector.The rotation matrix associated with q, denoted by R tb q , is given as follows [32]: where R tb q 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 velocity w of ∑ b relative to ∑ t is defined as follows: where w b and w t 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: Based on Coriolis's theorem, the following expression is obtained: Given that dw b /dt The attitude dynamics of the target are given by the following expression: where I t and N t 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 dw t /dt | t = 0 and N t = 0 are subsequently obtained, respectively.Eq. ( 10) is substituted into Eq.( 11) to yield the following expression: where w t , resolved into ∑ t , is given as w t = 0, 0, n T .Subsequently, the change in target's attitude is given as follows [32]: According to Eq. ( 6), the rotation matrix R tb k + 1 at the sampling instant k + 1 is a function of rotation quaternion q k + 1 .According to Eqs. ( 12) and ( 13), q k + 1 can be obtained through solving the numerical integration with the given w k and q k .Thereby, R tb k + 1 can be represented as a nonlinear function expression: 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.

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 Figure 2. The second includes solar panels on only one side where deployment fails on the other side as shown in Figure 3.
The x b -axis is along the center line of the solar panels, the y b -axis is along the symmetric axis of target's body and is orthogonal to the x b -axis, and the z b -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 y b -axis and the position vector of the DP is r b DP = 0, −l, 0 T resolved in ∑ b .The tracking error is defined as δr t = R tb t r t − r b DP , where r t 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: where ε is a constant denoting the upper bound of tracking error.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 Figure 4.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.
Furthermore, the scheme for discriminating collision avoidance between the original ellipsoid and the sphere is shown in Figure 5.
Without loss of generality, the original ellipsoid is as follows: 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]:

International Journal of Aerospace Engineering
The expanded ellipsoid envelope is as follows: It should be noted that the value of the minimum halfaxis 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: Proposition 1.Given the original ellipsoid in Eq. ( 16) and the sphere with radius r e , the sufficient condition for collision avoidance between the original ellipsoid and the sphere is defined as follows: 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 point p 0 using spherical coordinates θ, φ describing its direction, there are two other corresponding points located at the original ellipsoidal surface and the expanded ellipsoidal surface (points p 1 and p 2 , as shown in Figure 5), and their coordinates are given as follows: Subsequently, the tangent planes of the two points (planes L 1 and L 2 , as shown in Figure 5) are given as follows: According to Eq. ( 22), it is noted that the two tangent planes are parallel.Subsequently, the minimum distance from the point p 2 to the tangent plane L 1 is equal to the distance between tangent planes L 2 and L 1 , and this is given as follows: Given that a > b > c is defined, r c ≥ r e 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,  International Journal of Aerospace Engineering 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]: where u max 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 onorbit approaches.The velocity constraint is expressed as follows [27]:

25
where x max , y max , and z max denote the components of the maximum approach velocity in each direction.

Model Predictive Control Formulation
By using the dynamics and constraints discussed in Section 2, 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.

Prediction of the State
Variables.The predicted relative translational state sequence generated by Eq. ( 4) with state X k and control input U k is expressed as follows:

27
where X k + i | k denotes the predicted state variable at k + i with the given information X k , U k + i | k denotes the predicted control input at k + i, U k | k = U k , N p denotes the prediction horizon, and N c denotes the control horizon [12].As described in Section 2.2, the predicated rotation matrix R tb k + i at k + i can be solved offline with the given initial w 0 and q 0 .Thus, while implementing the MPC, the rotation matrix R tb k + 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: International Journal of Aerospace Engineering where Q i ∈ ℜ 3×3 denotes a positive-definite state-weighting matrix, R i ∈ ℜ 3×3 denotes a positive-definite controlweighting matrix, and In order to simplify the expression in Eq. ( 28), matrix C s = I Np ⊗ C d is introduced, where ⊗ denotes the Kronecker product of two matrices [27], and I Np denotes a N p -dimensional identity matrix.If we define , the objective function in Eq. ( 28) is expressed as follows:

Inequality Constraints.
In this section, the constraints discussed in Section 2 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
where x b , y b , z b 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 g m is defined as the collision threshold, and the collision avoidance is attained when g m > 0. Based on Eq. ( 31), the collision avoidance constraint for the MPC is expressed as follows:

32
where and H i is expressed as follows:

33
where I 3 denotes a three-dimensional identity matrix and E 3 denotes a three-dimensional zero matrix.After a few manipulations (see Appendix A), Eq. ( 32) is rewritten as follows:

34
where Based on Eq. ( 24), the control input constraint is reexpressed as follows: where I 3N c denotes a 3N c -dimensional identity matrix and The velocity constraint expressed in Eq. ( 25) is reexpressed as follows (see Appendix B): where the matrices D s and V max ∈ ℜ 6N p are defined as follows: min Given that the collision avoidance constraint is a nonconvex constraint, it is difficult to satisfy the convergence and the optimality of the solution to Problem 1.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
where U n s k denotes the n-iterated solution of U s k , and Proposition 2. Eq. (39) utilized as the collision avoidance constraint as opposed to Eq. (34), and Problem 1 is subsequently converted to Problem 2 as follows: Problem 2.
Subsequently, any U s k that is feasible for Problem 2 is also feasible for Problem 1.
Proof.Let U s k denote an arbitrary feasible point for Problem 2, i.e., Represent g i m U s k by the second-order Taylor series expansion at U n s k : where ∇ 2 g i m denotes the Hessian matrix of g i m .It should be noted that ∇ 2 g i m = 2F m i denotes a constant positivedefinite matrix where F m i = F i T L m F i is satisfied.We substitute Eq. (41) into Eq.( 42), and this results in g i m U s k > 0, and thus, any U s k that satisfies the condition in Eq. (39) also satisfies the condition in Eq. ( 34).The proof is completed.
It should be noted that Problem 2 is a quadratic programming problem [24].In order to optimize the control input U s k , sequential quadratic programming is utilized and a trust region between U n s k and U s k is introduced to ensure the convergence.
where D 0 ∈ ℜ 3N c 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.
min 44 9 International Journal of Aerospace Engineering Remark 2. Based on Proposition 2, 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 U 1 s k to commence.In the study, U 1 s k 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 vector X k , rotation quaternion q k , and relative angular velocity w k 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.

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 Table 1.
Considering the orbit parameters of the target described in Table 1, h is the orbit height, e is the eccentricity, Ω is the right ascension of ascending node, w is the argument of 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 a In our simulations, the sampling period T s is set to T s = 0 1 s.The weighting matrices are selected as Q i = I 3 and R i = I 3 .The maximum control input u max is set to u max = 0 1 m/s 2 , and the maximum relative velocity v max is set to v max = 0 5, 0 5, 0 5, 0 5, 0 5, 0 5 T m/s.The value of the prediction horizon N p is set to N p = 20, and the value of control horizon N c is set to N c = 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 D 0 = 0 04, 0 04, ⋯, 0 04, ⋯, 0 04 T .

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 states X 0 are set to 25, −25 3, 0 T m and 0, −0 0027, 0 T m/s.The initial angular velocities of the failed spacecraft are set to w b b = 2 2, 2 5, 2 3 T deg/s.The simulation time is set to 100 s.The simulation results are shown in Figures 6 and 7.
As shown in Figures 6(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/s 2 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 velocities w b b = 2 2, 2 5, 2 3 T deg/s (case A) as discussed in the previous simulation, initial angular velocities w b b = 1 2, 1 0, 0 8 T deg/s (case B) and w b b = 0, 0, 0 T deg/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.
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 Figure 8, 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 N c = 5, 10, 20 , the convergence time and total control input are calculated at different control horizons N p = 20, 30, 40, 50, 60 .The initial angular velocities of the target are also set to w b b = 2 2, 2 5, 2 3 T deg/s.The simulation results are shown in Figure 9.
As shown in Figure 9, with increases in the predictive horizon, the convergence time for approaching increases and the total control input decreases.It is noted that the

Parameter Value
Orbit parameter h = 35786 km, e = 0,  28) is consist of the tracking error and the control effort.The increase of predictive horizon N p 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 N c , 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 onorbit 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 Section 4. 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 r 0 = 25, -25 3, 0 T m (Sat 1), r 0 = 25, 25  As shown in Figure 10, 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.
As shown in Figure 11, 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.
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 states X k can be represented as   where e k = e p k ; e v k and e p k and e v k 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 Û k applied to the servicer can be represented as where Λ k = diag Λ x k , Λ y k , Λ z k is a threedimensional diagonal matrix.
Combining Eqs. ( 4), (45), and (46), the prediction of relative states should be represented as However, navigation and control noise cannot be measured.Then, Eq. ( 4) with measurement X k and control input U k 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 r 0 = 25, -25 3, 0 T m (Sat 1).The initial angular velocities of the target are also set to w b b = 2 2, 2 5, 2 3 T deg/s.A Monte-Carlo analysis is performed to illustrate the influence of uncertainties on the MPC performance.It is assumed that elements in e p , e v , and Λ follow Gaussian distribution, and their accuracies (3σ) are defined as E p , E v , and E u , respectively.When developing Monte-Carlo analysis, 1000 simulations in each of two conditions such as E p , E v , E u = 0 005 m, 10 −4 m/s, 0 01 and E p , E v , E u = 0 02 m, 10 −3 m/s, 0 02 are conducted.Reconsidering the tracking error defined in Eq. ( 15), the average tracking   International Journal of Aerospace Engineering error δr is considered to evaluate the influence of uncertainties, which is defined as where t f = 100 s is the end time of simulation and t s = 80 s is designed to start recording the tracking error.The Monte-Carlo simulation results regarding the average tracking error δr are shown in Figure 13.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 Figure 13, when the uncurtains are E p , E v , E u = 0 005 m, 10 −4 m/s, 0 01 and E p , E v , E u = 0 02 m, 10 −3 m/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.

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 timevariant 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 Eq. ( 32) can be rewritten as follows: Then, define G m i as G m i = R tb k + i H i C s A s X k − d m , 0, 0 T and F i as F i = R tb k + i H i C s B s , Eq. (A.1) can be rewritten as follows: Then, define G m i as G m i = 2G m i T L m F i and F m i as F m i = F i T L m F i , Eq. (A.2) can be expanded as follows: It is noted that Eq. (A.3) is the same as Eq.(34).
B. The Derivation of Eq. (36) Eq. ( 25) can be rewritten as follows: where the matrices v max and D d are defined as follows: For the implementation of MPC, the velocity sequence constraints can be represented as follows: where Substituting Eq. (B.4) into Eq.(B.3), Eq. (B.3) can be rewritten as follows: It is noted that Eq. (B.5) is the same as Eq.(36).
C. The Derivation of Eq. (44) Eq. ( 40) can be rewritten as follows: Considering the inequality constraint about the trust region, Eq. (43) can be rewritten as follows: It is noted that Eq. (C.3) is the same as Eq. ( 44).

2. 3 .
Target Body-Fixed Coordinate System ∑ b .The origin of the target body-fixed coordinate system o − x b y b z b is in the c.m. of the target.Without loss of generality, it is assumed that x b , y b , and z b are aligned with its principal axes.

Figure 4 :
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.

Figure 5 :
Figure 5: Scheme for discriminating collision avoidance between the original ellipsoid and the sphere.

max T 37 3 . 4 .
MPC by Using Sequential Convex Programming.The design of the MPC for the approach problem is summarized as follows: International Journal of Aerospace Engineering Problem 1.
1 = 2 m, b 1 = 1 8 m, c 1 = 1 5 m, and d 1 = 0m.The parameters of the ellipsoid envelope for the first failed spacecraft's solar panels are set to a 2 = 6 m, b 2 = 1 2 m, c 1 = 1 2 m, and d 2 = 0m.The parameters of the ellipsoid envelope for the second failed spacecraft's solar panels are set to a 2 = 2 m, b 2 = 1 2 m, c 1 = 1 2 m, and d 2 = 4 m.The radius of the sphere envelope for the servicer spacecraft is r e = 0 5 m.

Figure 6 :
Figure 6: Temporal response of the relative position, relative velocity, control input, and collision threshold.(a) Temporal history along the x 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.

Figure 7 :
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.

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.

Figure 9 :Figure 10 :
Figure 9: Influences of the predictive horizon N p and control horizon N c on (a) convergence time and (b) total control input.

Figure 11 :Figure 12
Figure 11: Temporal response of the collision threshold of each servicer.

1 ,
v max = x max y max z max x max y max z max T B 2 | b and dw t /dt | I = dw t /dt | t , a combination of Eqs.(8) and (9) yields the following expression:

Table 1 :
Parameters of the failed spacecraft.
Then, substituting Eq. (C.2) into Eq.(C.1),Problem 2 that considers the trust region is subsequently converted into Problem 3:U s k Γ k = U s k T W k U s k + P k U s k + M k T QM k subject to I 3N c ; −I 3N c U s k ≤ U max I 3N p ; −I 3N p G s B s U s k ≤ V max − I 3N p ; −I 3N p G s A s X k − ∇g i min