2. DYNAMICS OF SPACECRAFT RELATIVE MOTION

This section illustrates the spacecraft relative motion dynamics. The Local Vertical Local Horizontal (LVLH) coordinate (Cao et al., Reference Cao, Chen and Misra2014; Ou and Zhang, Reference Ou and Zhang2017a) is utilised to describe the spacecraft relative motion. The LVLH coordinate is centred at the centroid of the chief spacecraft, the x-axis points out radially from the centre of the Earth to the mass centre of the chief spacecraft, the y-axis is aligned in the direction of in-track motion and the z-axis completes the right-handed coordinate system. Furthermore, as the in-plane is decoupled from motion normal to the orbital plane, the six-dimensional system can be reduced into a second-order system for the motion normal to the nominal orbit plane and a fourth-order system for the in-plane motion. Then, only the coplanar problem is considered in this paper and hence the relative motion is given as:

(1)$$\left\{\matrix{\ddot{x}_t - 2n \dot{y}_t - 3 n^2 x_t = u_{x-t} \hfill \cr \ddot{y}_t + 2n \dot{x}_t = u_{y-t} \hfill}\right.$$

where r_1−t = [x _t, y _t]^T and v_1−t = [v _x−t, v _y−t]^T represent the relative positions and relative velocities in LVLH coordinates, respectively. u_1−t = [u _x−t, u _y−t]^T represents the control acceleration, μ is the gravitational constant of the Earth, a represents the semi-major axis of the orbit of the chief spacecraft and n is the angular velocity, where $n= \sqrt{\mu/a^{3}}$.

Then, the forces along both planar directions are assumed to be independent and we obtain the state dynamics matrix A and B, which is shown as:

(2)$$\dot{\bi X} = {\bi AX} + {\bi Bu}$$

where ${\bi X} = [{\bi r}_{1-t}^{\rm T} \quad {\bi v}_{1 - t}^{\rm T}]^{\rm T}$ is the state vector:

(3)$${\bi A} =\left[\matrix{0 & 0 &1 &0 \cr 0 & 0 & 0 & 1 \cr 3n^2 &0 &0 &2n \cr 0 & 0 & -2n & 0} \right] \cdot {\bi B} = \left[\matrix{0 &0 &1 &0 \cr 0 & 0 & 0 & 1} \right]^{\rm T}$$

Furthermore, as only the circular orbit of the chief spacecraft is considered in this paper, we obtain the analytical solution X (t) as:

(4)$${\bi X} (t) = {\bf \varPhi} (t_0, t) {\bi X} (t_0) + \int_{t_0}^{t} {\bf \varPhi}_v (s, t) {\bi u}(s) ds$$

where t ₀ is the initial time and the transition matrix ${\bf \varPhi} (t_{0}, t)$ is given by:

(5)$${\bf \varPhi} (t_0, t) = \left[\matrix{{\bf \varPhi}_{rr} (t_0, t) &{\bf \varPhi}_{rv} (t_0, t) \cr {\bf \varPhi}_{vr} (t_0, t) &{\bf \varPhi}_{vv} (t_0, t)} \right], {\bf \varPhi}_v (t_0, t) = \left[\matrix{{\bf \varPhi}_{rv} (t_0, t) \cr {\bf \varPhi}_{vv} (t_0, t)} \right]$$

where ${\bf \varPhi}_{rr} (t_{0}, t), {\bf \varPhi}_{rv} (t_{0}, t), {\bf \varPhi}_{vr} (t_{0}, t), {\bf \varPhi}_{vv} (t_{0}, t)$ are obtained as follows:

(6)$$\eqalign{{\bf \varPhi}_{rr} (t_0, t) & = \left[\matrix{4 -3 \cos (n (t - t_0)) & 0\cr 6 (\sin (n (t - t_0)) - n (t- t_0)) & 1} \right], \cr {\bf \varPhi}_{rv} (t_0, t) &= \displaystyle{1 \over n} \left[\matrix{\sin (n(t - t_0))& 2 (1- \cos (n(t - t_0)))\cr -2(1-\cos (n(t-t_0))) & 4 \sin (n (t-t_0)) - 3n (t - t_0)}\right], \cr {\bf \varPhi}_{vr} (t_0, t)&= \left[\matrix{3n \sin (n(t-t_0)) & 0\cr 6n (\cos (n(t-t_0))-1) & 0} \right], \cr {\bf \varPhi}_{vv} (t_0, t) &= \left[\matrix{\cos (n(t -t_0)) & 2 \sin (n(t - t_0)) \cr -2 \sin (n(t - t_0)) &4 \cos (n(t -t_0))-3} \right].}$$

3. THE ECPC METHOD

Wang et al. (Reference Wang, Bai, Ran, Zhao, Zhang and Chen2017; Reference Wang, Bai, Xing, Gianmarco, Ni and Chen2018) developed an ECPC method for spacecraft close-range safe proximity in the presence of uncertainties. Similar to the concept of an isobaric curve, the ECPC represents a curve consisting of equal collision probability points around the chief spacecraft. This proposed safe close-range proximity scheme is attractive for two reasons. First, the ECPC method does not need to solve the transcendental function after the approximate treatment, so that the computational burden can be greatly decreased while maintaining the necessary accuracy. More importantly, the effectiveness of the proposed method can be theoretically verified.

Considering the navigation and the control uncertainties, the linear covariance method is utilised to propagate the uncertainties and we obtain the covariance matrix with respect to the state error, which is given as (Yang et al., Reference Yang, Luo, Zhang and Tang2016; Luo and Yang, Reference Luo and Yang2017):

(7)$$\delta {\bi X} = {\bi X} - {\bf E} ({\bi X})$$

(8)$$\eqalign{{\bi C}_{\delta X} &={\bf E} \lcub [\delta {\bi X} - {\bf E} (\delta {\bi X})] [\delta {\bi X} - {\bf E} (\delta {\bi X})]^{\rm T} \rcub \cr &={\bf \varPhi} (t_0, t) {\bi C}_{\delta X_{0}}{\bf \varPhi} (t_0, t)^{\rm T} + \sum_{i =1}^{N} {\bf \varPhi}_v (t_i, t) {\bi C}_{\delta v_{i}} {\bf \varPhi}_v (t_i, t)^{\rm T}}$$

where N is the number of control impulses, δ X is the error of the X, and E(X) and E (δ X) represent the mean values of X and δ X, respectively. ${\bi C}_{\delta x_{0}}$ and ${\bi C}_{\delta v_{i}}$ are the covariance matrices of the initial navigation and control uncertainties in LVLH coordinates, respectively and C_δX is the uncertainty covariance matrix of the state vector.

Wang et al. (Reference Wang, Bai, Ran, Zhao, Zhang and Chen2017) developed an equal collision probability cure, which is a Gaussian form and can be described as follows:

(9)$$\left\{\matrix{if \displaystyle{x_{1-t}^2 \over \sigma_{x-t}^2} + \displaystyle{y_{1-t}^2 \over \sigma_{y-t}^2} \lt D_0 \& v_{1- t}^{paral} \gt 0 \hfill \cr V_0 ({\bi r}_{1-t}) = \lambda_0 \left\{\exp \left[ -\displaystyle{1 \over 2} \left(\displaystyle{x_{1-t}^2 \over \sigma_{x-t}^2} + \displaystyle{y_{1-t}^2 \over \sigma_{y-t}^2} \right) \right] - \exp \left(-\displaystyle{1 \over 2} D_0^2\right)\right\}, \hfill \cr else \hfill \cr V_0 ({\bi r}_{1-t}) = 0. \hfill} \right.$$

where ${\bi C}_{\delta r_{1-t}} = \hbox{diag} ([\sigma_{x-t}^{2} \quad \sigma_{y-t}^{2}])$ is the uncertainty covariance matrix of the relative position. λ₀ is a positive constant that shapes the repulsive potential magnitude, D ₀ is the radius of the hazardous zone and $v_{1-t}^{paral}$ is the value of relative parallel velocity, which is given as:

(10)$$v_{1-t}^{paral} = {\bi v}_{1-t} \cdot {\bi n}_{1-t}^{paral}$$

where ${\bi n}_{1-t}^{paral}$ is the unit vector pointing from the deputy spacecraft to the chief spacecraft and ${\bi n}_{1-t}^{paral} = -{\bi r}_{1-t}/\vert {\bi r}_{1-t}\vert$. If ${\bi v}_{1 - t}^{paral} \leq 0$, the deputy spacecraft is moving away from the chief spacecraft and we do not execute the avoidance manoeuvres. Otherwise, the deputy spacecraft is moving toward the chief spacecraft and the avoidance manoeuvres are implemented.

The hazardous zone is defined as the zone where a collision is likely to happen. Furthermore, the hazardous zone is a sphere around the chief spacecraft and has the following radius:

(11)$$D_0 = d_0 (R_0 + D_s)$$

where d ₀ is a positive constant and D _s is defined as the braking distance taken to stop once the avoidance manoeuvres are implemented and is computed by:

(12)$$D_s = \displaystyle{(v_{1-t}^{paral})^2 \over 2a_{\max}}$$

where a _max is the maximum control acceleration of the deputy spacecraft. The perpendicular relative velocity of the deputy spacecraft is obtained from:

(13)$${\bi v}_{1-t}^{perpen} {\bi n}_{1-t}^{perpen} = {\bi v}_{1-t} - {\bi v}_{1-t}^{paral} {\bi n}_{1-t}^{paral}$$

where ${\bi v}_{1-t}^{perpen}$ is the value of the perpendicular relative velocity and ${\bi n}_{1-t}^{perpen}$ is the unit vector perpendicular to ${\bi n}_{1-t}^{paral}$, which is as follows:

(14)$${\bi n}_{1-t}^{perpen} =-\displaystyle{1 \over \sqrt{x_{1-t}^2 + y_{1-t}^2}} \left[\matrix{-y_{1-t} \cr x_{1-t}} \right].$$

According to Equation (9), we define the negative gradient of the repulsive potential with respect to the position terms as the avoidance manoeuvres (Ge and Cui, Reference Ge and Cui2002):

(15)$${\bi F}_{repel} = - \nabla_{{\bi r}_{1-t}} V_0 ({\bi r}_{1-t})$$

where:

(16)$$\nabla_{{\bi r}_{1-t}} V_0 ({\bi r}_{1-t}) = \displaystyle{\partial V_0 ({\bi r}_{1-t}) \over \partial r_{1-t}} \displaystyle{\partial r_{1-t} \over \partial {\bi r}_{1-t}} + \displaystyle{\partial V_0 ({\bi r}_{1-t}) \over \partial D_0} \displaystyle{\partial D_0 \over \partial {\bi r}_{1-t}}$$

With the help of Equations (9), (15) and (16), we obtain the collision avoidance manoeuvre as:

(17)$${\bi F}_{repel} = {\bi F}_{oparal} + {\bi F}_{perpen}$$

where

(18)$$M_1 = \exp \left[-\displaystyle{1 \over 2} \left(\displaystyle{x_{1-t}^{2} \over \sigma_{x-t}^2} + \displaystyle{y_{1-t}^2 \over \sigma_{y-t}^2}\right)\right], M_2 = \exp \left(-\displaystyle{1 \over 2} D_0^2\right)$$

(19)$${\bi F}_{oparal} = -\lambda_0 r_{1-t} M_1 \left(\displaystyle{1 \over \sigma_{x-t}^2} + \displaystyle{1 \over \sigma_{y-t}^2} \right) {\bi n}_{1-t}^{paral}, {\bi F}_{perpen} = \lambda_0 D_0 M_2 d_0 \displaystyle{v_{1-t}^{paral} \over a_{\max}} \displaystyle{v_{1-t}^{perpen} \over r_{1-t}}{\bi n}_{1-t}^{perpen}.$$

4. THE MECPC METHOD

As the collision probability varies with different shapes (Bai et al., Reference Bai, Ma, Chen and Tang2016), the shape has great influence on the traditional ECPC method. However, in the ECPC method, Wang et al. (Reference Wang, Bai, Ran, Zhao, Zhang and Chen2017; Reference Wang, Bai, Xing, Gianmarco, Ni and Chen2018) assumed that the geometrical shapes of the two spacecraft were spheres but most space objects are not actually spheres. Therefore, if the spacecraft shape is assumed as a sphere in some scenarios, the safety performance of the mission deteriorates, as a collision will normally lead to mission failure. To address this problem, the MECPC method is proposed to guarantee convex polygon shape spacecraft safe proximity with other spacecraft.

Assuming that the shape of the chief spacecraft is a convex polygon, such as a square, the chief spacecraft's geometrical shape can be divided into one maximum inside envelope circle part and several arbitrary parts. Moreover, each arbitrary shape part has its corresponding minimum exterior envelope circle. Once the deputy spacecraft flies around the chief spacecraft, the function modules, which are the maximum inside the envelope circle part with respect to the main circle and the minimum exterior envelope circle part with respect to the arbitrary shape parts, jointly generate avoidance manoeuvres on the deputy spacecraft with the basic ECPC method. Figure 1 shows the defined divisions of the convex polygon shape chief spacecraft. As shown in Figure 1, the vertical coordinates system is divided into parts I II, III and IV. The convex polygon shape of the chief spacecraft is also separated into five parts, and part 5 is the maximum inside envelope circle. When the deputy spacecraft flies around the chief spacecraft, the function module is selected by:

Figure 1.

The defined divisions of the convex polygon shape chief spacecraft.

(20)$$\left\{\matrix{Part\,1\,and\,Part\,5, h = 1 & x\,\gt 0\,and\,y\,\gt 0 \hfill \cr Part\,2\,and\,Part\,5, h = 2 & x \lt 0\,and\,y\,\gt 0 \hfill \cr Part\,3\,and\,Part\,5, h=3 & x\,\lt\,0\,and\,y\,\lt 0 \hfill \cr Part\,4\,and\,Part\,5, h=4 &x\,\gt 0\,and\,y\,\lt 0 \hfill \cr Part\,5, h=5 \hfill & x = 0\,or\,y = 0 \hfill} \right.$$

Then, if the expected encounter position ${\bi r}_{ep} = \left[\matrix{x_{ep} & y_{ep}} \right]^{\rm T}$, which is located in the outside envelope surface, the following two conditions are analysed:

$${\rm (i)}\,x_{ep} = 0\,or\,y_{ep} = 0$$

Figure 2 illustrates the force analysis of the deputy spacecraft in Section I. From Figure 2, part 5 and the minimum exterior envelope circle of part 1 generates forces F_50-repel and F_10-repel, respectively. Moreover, the O ₁x ₁y ₁ coordinate system is established, which has its x ₁ axis radiating from the origin in the LVLH coordinates toward the centre of the minimum exterior envelope circle located at the centre. Rotating the x ₁ axis counter-clock-wise along the Oz axis for 90°, we obtain the y ₁ axis.

Figure 2.

The force analysis of the deputy spacecraft in Section I.

As shown in Figure 2, any points r_1−t in the LVLH coordinates are transformed to a relative position in the O ₁x ₁y ₁ coordinates, which is given as:

(21)$${\bi r}_{10-t} = {\bi W}_1 \lcub {\bi r}_{1 - t} - [x_{10} y_{10}]^{\rm T} \rcub $$

where:

(22)$${\bi W}_1 = \left[\matrix{\cos (\theta) &\sin (\theta)\cr {-}\sin (\theta) & \cos(\theta)}\right], \sin (\theta) = \displaystyle{y_{10} \over \sqrt{x_{10}^{2} + y_{10}^2}}, \cos (\theta) = \displaystyle{x_{10} \over \sqrt{x_{10}^2 + y_{10}^2}}.$$

Assuming the attitude of the chief spacecraft is fixed, the relative velocity is transformed from the LVLH coordinates to the O ₁x ₁y ₁ coordinates:

(23)$${\bi v}_{10 - t} ={\bi W}_{1}{\bi v}_{t}.$$

From Equation (21), the covariance matrices ${\bi C}_{\delta r_{1-t}}$ are defined in LVLH coordinates and are transformed to the O ₁x ₁y ₁ coordinates:

(24)$${\bi C}_{\delta r _{10-t}} = {\bi W}_{1} {\bi C}_{\delta r_{1-t}} {\bi W}_{1}^{-1}.$$

The transformation from the O ₁x ₁y ₁ frame to the O ₁x _1dy _1d frame in which ${\bi C}_{\delta r_{10-t}}$ is diagonalized is given by ${\bi G}_{1-t} {\bi W}_{1}^{-1}$:

(25)$${\bi C}_{\delta r _{10-td}} ={\bi G}_{1-t}{\bi W}_{1}^{-1} {\bi C}_{\delta r _{10-t}} \lpar {\bi G}_{1-t} {\bi W}_{1}\rpar ^{-1} = \left[\matrix{\sigma_{x-td}^2 & 0\cr 0 & \sigma_{y-td}^{2}} \right]$$

where G_1−t is an orthogonal transformation.

According to Equations (21), (23) and (25), the relative position and relative velocity in the O ₁x _1dy _1d coordinates are given as:

(26)$${\bi r}_{10-td} = {\bi G}_{1-t} \lcub {\bi r}_{1-t} - [x_{10} \quad y_{10}]^{\rm T} \rcub , {\bi r}_{10-td} = {\bi G}_{1-t} {\bi v}_{1-t}$$

where r_10−td = [x _10−td, y _10−td]^T and v_10−td = [v _x10−td, v _y10−td]^T represent the relative position and the relative velocity in the O ₁x _1dy _1d coordinates, respectively.

From Equations (9)–(26), when the deputy spacecraft flies in Section I, the collision avoidance manoeuvres generated by Part 1 in the LVLH coordinates are given as:

(27)$${\bi F}_{10-repel} = {\bi F}_{10-paral} + {\bi F}_{10-perpen}$$

where D ₀₁₀ is radius of the hazardous zone around Part 1 in the O ₁x _1dy _1d coordinates. r _10−td is the value of relative position in the O ₁x _1dy _1d coordinates, v _10−td^paral and v _10−td^perpen are the value of relative parallel velocity and relative perpendicular velocity in the O ₁x _1dy _1d coordinates and n_10−td^paral and n_10−td^perpen represent the unit vector opposite and perpendicular to r_1−td, respectively.

(28)$$M_{110d} = \exp \left[ - \displaystyle{1 \over 2} \left(\displaystyle{x_{10-td}^2 \over \sigma_{x-td}^2} + \displaystyle{y_{10-td}^2 \over \sigma_{y-td}^2} \right) \right], M_{210d} = \exp \left( - \displaystyle{1 \over 2} D_{010d}^{2} \right)$$

(29)$${\bi F}_{10-paral} =-\lambda_0 r_{1-d} M_{110d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right) {\bi n}_{10-td}^{paral}$$

(30)$${\bi F}_{10-perpen} = \lambda_0 D_{010} M_{210}d_0 \displaystyle{v_{10-td}^{paral} \over a_{\max}} \displaystyle{v_{10-d}^{perpen} \over r_{10-td}} {\bi n}_{10-td}^{perpen}.$$

With the aid of Equations (9), (10), (14) and Equations (21)–(22), we obtain:

(31)$${\bi r}_{10} = a_{10-1} {\bi n}_{1-t}^{paral} + a_{10-2} {\bi n}_{1-t}^{perpen}$$

where:

(32)$$a_{10-1} = -\displaystyle{x_{10}x_{1-t} + y_{10} y_{1-t} \over \sqrt{x_{1-t}}^2 + y_{1-t}^2}, a_{10-2} = \displaystyle{x_{10} y_{1-t} y_{10} x_{1-t} \over \sqrt{x_{1-t}^2 + y_{1-t}^2}}$$

According to Equations (9)–(32), we obtain:

(33)$${\bi n}_{10-td}^{paral} = \displaystyle{r_{1-t} + a_{10-1} \over r_{10-td}} {\bi n}_{1-td}^{paral} + \displaystyle{a_{10-2} \over r_{10-td}} {\bi n}_{1-td}^{perpen}.$$

Based on Equation (33), we have

(34)$${\bi n}_{10-td}^{perpen} = \displaystyle{a_{10-2} \over +} r_{10-td} {\bi n}_{1-td}^{paral} + \displaystyle{r_{1-t} + a_{10-1} \over r_{10-td}} {\bi n}_{1-td}^{perpen}.$$

From Equations (33) and (34), we obtain:

(35)$${\bi F}_{10-paral} = -\left[\matrix{\lambda_0 r_{10-td} M_{110d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right) \displaystyle{r_{1-t} + a_{10-1} \over r_{10-td}} \cr +\lambda_0 D_{010d} M_{210d} d_0 \displaystyle{v_{10-td}^{paral} \over a_{\max}} \displaystyle{v_{10-td}^{perpen} \over r_{10-td}} \displaystyle{a_{10-2} \over r_{10-td}}}\right] {\bi n}_{1-td}^{paral}$$

(36)$${\bi F}_{10-perpen} = \left[\matrix{\lambda_0 r_{10-td} M_{110d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right) {a_{10-2} \over r_{10-td}} \cr + \lambda_0 D_{010d} M_{210d} d_0 \displaystyle{v_{10-td}^{paral} \over a_{\max}} \displaystyle{v_{10-td}^{perpen} \over r_{10-td}} \displaystyle{r_{1-t}+ a_{10-1} \over r_{10-td}}}\right] {\bi n}_{1-td}^{perpen}$$

According to the above deduction, the collision avoidance manoeuvres generated by Part 2, Part 3 and Part 4 with respect to Section II, Section III and Section IV are given as:

(37)$${\bi F}_{i0-repel} = {\bi F}_{i0-paral} + {\bi F}_{i0-perpen}, (i=2,3,4)$$

where F _i0−repel (i = 2, 3, 4) represents the collision avoidance manoeuvres generated by Part 2, Part 3 and Part 4. r_i0 = [x _i0  y _i0]^T, (i = 2, 3, 4) represents the different relative position of the minimum exterior circle's centre with respect to the corresponding part. r_i0−td = [x _i0−td, y _i0−td]^T and v_i0−td = [v _xi0−td, v _yi0−td]^T(i = 2, 3, 4) represent the relative position and the relative velocity in the corresponding coordinates, respectively. D _0i0, (i=2, 3, 4) represents the radius of the corresponding hazardous zones, r_i0−td, (i = 2, 3, 4) are values of relative positions in the corresponding coordinates and v _i0−td^paral, (i = 2,3,4) and v _i0−td^perpen, (i = 2, 3, 4) are the values of relative parallel velocity and relative perpendicular velocity in the corresponding coordinates.

(38)$${\bi r}_{i0-td} = {\bi G}_{1-t} \lcub {\bi r}_{1-t} - [x_{i0} \quad y_{i0}]^{\rm T}\rcub , {\bi v}_{i0-td} = {\bi G}_{1-t} {\bi v}_{1-t}, (i = 2,3,4)$$

(39)$$M_{1i0d} = \exp \left[- \displaystyle{1 \over 2} \left(\displaystyle{x_{i0 - td}^2 \over \sigma_{x-td}^2} + \displaystyle{y_{i0-td}^2 \over \sigma_{y-td}^2} \right)\right], M_{2i0d} = \exp \left( - \displaystyle{1 \over 2} D_{0i0d}^2 \right), (i = 2,3,4)$$

(40)$$a_{i0-1} = -\displaystyle{x_{i0}x_{1-t} + y_{i0} y_{1-t} \over \sqrt{x_{1-t}^2 + y_{1-t}^2}}, a_{i0-2} = -\displaystyle{x_{i0} y_{1-t} + y_{i0} x_{1-t} \over \sqrt{x_{1-t}^2 + y_{1-t}^2}}, (i = 2,3,4)$$

(41)$${\bi F}_{10-paral} = -\left[\matrix{\lambda_0 r_{i0-td} M_{1i0d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right) \displaystyle{r_{1-t} + a_{i0-1} \over r_{i0-td}}\cr + \lambda_0 D_{0i0d} M_{2i0d} d_0 \displaystyle{v_{i0-td}^{paral} \over a_{\max}} \displaystyle{v_{i0-td}^{perpen} \over r_{i0-td}} \displaystyle{a_{i0-2} \over r_{i0-td}}}\right] {\bi n}_{1-td}^{paral}, (i = 2,3,4)$$

(42)$${\bi F}_{10-perpen} =- \left[\matrix{\lambda_0 r_{i0-td} M_{1i0d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right) \displaystyle{a_{i0-2} \over r_{i0-td}} \cr +\,\lambda_0 D_{0i0d} M_{2i0d} d_0 \displaystyle{v_{i0-td}^{paral} \over a_{\max}} \displaystyle{v_{i0-td}^{perpen} \over r_{i0-td}} \displaystyle{r_{1-t} + a_{i0-1} \over r_{i0-td}}} \right] {\bi n}_{1-td}^{perpen}, (i = 2,3,4)$$

From Equations (9)–(19), the collision avoidance manoeuvres generated by part 5 are given as follows:

(43)$${\bi F}_{50{\rm-}repel} = {\bi F}_{50{\rm-}paral} + {\bi F}_{50{\rm-}perpen}$$

where F _50-repel represents the collision avoidance manoeuvres generated by part 5. r_50-td = [x _50-td, y _50-td]^T and v_50-td = [v _x50-td, v _y50-td]^T represent the relative position and the relative velocity in the O ₁x _1dy _1d coordinates, respectively. D ₀₅₀ is the radius of the corresponding hazardous zone, r _50-td is the value of the relative positions and v _50-td^paral and v _50-td^perpen are the values of the relative parallel velocity and relative perpendicular velocity in the O ₁x _1dy _1d coordinates.

(44)$${\bi r}_{50{\rm-}td} = {\bi G}_{1-t} {\bi r}_{1-t}, {\bi v}_{50{\rm-}td} = {\bi G}_{1-t} {\bi v}_{1-t}$$

(45)$$M_{150d} = \exp \left[- \displaystyle{1 \over 2} \left(\displaystyle{x_{50{\rm-}td}^2 \over \sigma_{x-td}^2} + \displaystyle{y_{50{\rm-}td}^2 \over \sigma_{y-td}^2}\right)\right], M_{250d} = \exp \left( - \displaystyle{1 \over 2} D_{050d}^2\right)$$

(46)$${\bi F}_{50{\rm-}paral} = -\lambda_0 r_{50{\rm-}td} M_{150d} \left(\displaystyle{1 \over \sigma_{x-td}^2} \displaystyle{1 \over \sigma_{y-td}^2}\right) {\bi n}_{1-td}^{paral}$$

(47)$${\bi F}_{50{\rm-}perpen} = -\lambda_0 D_{050{\rm-}td} M_{250d}d_0 \displaystyle{v_{50{\rm-}td}^{paral} \over a_{\max}} \displaystyle{v_{50{\rm-}td}^{perpen} \over r_{50{\rm-}td}} {\bi n}_{1-td}^{perpen}$$

Thus, the total collision avoidance manoeuvres generated are given as:

(48)$$\eqalign{{\bi F}_{Total-repel} &= \left\{\matrix{{\bi F}_{50{\rm-}repel} + {\bi F}_{h0-repel}, \hfill &h = 1,2,3,4 \hfill \cr {\bi F}_{50{\rm-}repel}, \hfill & h=5 \hfill}\right. \cr &\quad (\hbox{ii}) x_{ep} \neq 0 and y_{ep} \neq 0}$$

Figure 3 illustrates the force analysis of the deputy spacecraft in Section I. As we can see from Figure 3, part 5 generates force F_50-repel. By the encounter point r_ep, part 1 is divided into two parts, and each part has its own minimum exterior envelope circle. Subsequently, the O _1-i x _1-i y _1-i (i=1, 2) coordinates are established, which have their x _1-i axis pointing out from the origin in the LVLH coordinates toward the centre of the minimum exterior envelope circle located at the centre. The y _1-i axis of the O _1-ix _1-i y _1-i (i=1, 2) frame has the same direction as the x _1-i axis and is rotated 90° counter-clock-wise along the Oz axis.

Figure 3.

The force analysis of the deputy spacecraft in Section I and the condition h ₁=1.

In this section, the following condition is divided into the following four parts:

(49)$$\left\{\matrix{h_1 =1 & x_{ep} \gt 0\,and\,y_{ep} \gt 0\cr h_1 =2 & x_{ep} \lt 0\,and\,y_{ep} \gt 0 \cr h_1 =3 & x_{ep} \lt 0\,and\,y_{ep} \lt 0 \cr h_1 =4 & x_{ep} \gt 0\,and\,y_{ep} \lt 0}.\right.$$

If the condition h ₁ =1 is satisfied, Figure 3 illustrates the force analysis of the deputy spacecraft in Section I. From Figure 3, the minimum exterior envelope circle of part 1-1 and part 1-2 generate the repulsive force F_11-repel and F_12-repel, respectively. Furthermore, both the O _1-1x _1-1y _1-1 coordinates and the O _1-2x _1-2y _1-2 coordinates are established as the O ₁x ₁y ₁ coordinates depicted in Figure 2.

From Equations (9)–(36), the collision avoidance manoeuvres generated with respect to the corresponding part are given as:

(50)$${\bi F}_{ij-repel} = {\bi F}_{ij-paral} + {\bi F}_{ij-perpen} (i = 1, 2, 3, 4, j =1,2)$$

where F_ij−repel (i = 1, 2, 3, 4, j = 1, 2) represent the collision avoidance manoeuvres generated by the corresponding part. r_ij = [x _ij  y _ij]^T, (i = 1, 2, 3, 4, j = 1, 2) represent the different relative position of the minimum exterior circle's centre with respect to the corresponding part. r_ij−td = [x _ij−td,   y _ij−td]^T, v_ij−td = [v _xij−td,   v _yij−td]^T(i = 1, 2, 3, 4, j = 1, 2) represent the relative position and the relative velocity in the corresponding coordinates, respectively. D _0ij, (i = 1, 2, 3, 4, j = 1, 2) are the radii of the corresponding hazardous regions, r _ij−td, (i = 1, 2, 3, 4, j = 1, 2) are values of relative positions in the corresponding coordinates and v _ij-td^paral, (i = 1, 2, 3, 4, j = 1,2) and v _ij-td^perpen, (i = 1, 2, 3, 4, j = 1,2) are the values of relative parallel velocity and relative perpendicular velocity in the corresponding coordinates.

(51)$${\bi r}_{ij-td} = {\bi G}_{1-t} \left\{{\bi r}_{1-t} - \left[\matrix{x_{ij} & y_{ij}} \right]^{\rm T}\right\}, {\bi v}_{ij-td} = {\bi G}_{1-t} {\bi v}_{1-t}, (i = 1,2,3,4, j = 1,2)$$

(52)$$M_{1ijd} = \exp \left[-\displaystyle{1 \over 2} \left(\displaystyle{x_{ij-td}^2 \over \sigma_{x-td}^2} + \displaystyle{y_{ij-td}^2 \over \sigma_{y-td}^2}\right)\right], M_{2ijd} = \exp \left(-\displaystyle{1 \over 2} D_{0ijd}^2 \right) (i = 1,2,3,4, j = 1,2)$$

(53)$$a_{ij-1} = -\displaystyle{x_{ij}x_{1-t} + y_{ij} y_{1-t} \over \sqrt{x_{1-t}^2 + y_{1-t}^2}}, a_{ij-2} = -\displaystyle{x_{ij}y_{1-t} + y_{ij} x_{1-t} \over \sqrt{x_{1-t}^2 + y_{1-t}^2}}, (i = 1,2,3,4, j = 1,2)$$

(54)$${\bi F}_{ij-paral} =- \left[\matrix{\lambda_0 r_{ij-td} M_{1ijd} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right) \displaystyle{r_{1-t} + a_{ij-1} \over r_{ij-td}}\cr +\lambda_0 D_{0ijd} M_{2ijd} d_0 \displaystyle{v_{ij-td}^{paral} \over a_{\max}} \displaystyle{v_{ij-td}^{perpen} \over r_{ij-td}} \displaystyle{a_{ij-2} \over r_{ij-td}}}\right] {\bi n}_{1-td}^{parpal}, (i = 1,2,3,4, j = 1,2)$$

(55)$${\bi F}_{ij-perpen} = - \left[\matrix{-\lambda_0 r_{ij-td} M_{1ijd} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2}\right) \displaystyle{a_{ij-2} \over r_{ij-td}} \cr + \lambda_0 D_{0ijd} M_{2ijd} d_0 \displaystyle{v_{ij-td}^{paral} \over a_{\max}} \displaystyle{v_{ij-td}^{perpen} \over r_{ij-td}} \displaystyle{r_{1-t} +a_{ij-1} \over r_{ij-td}}} \right] {\bi n}_{1-td}^{perpen}, (i = 1,2,3,4, j = 1,2).$$

Thus, the total collision avoidance manoeuvres are obtained as:

(56)$${\bi F}_{Total-repel} = \left\{\matrix{{\bi F}_{50{\rm-}repel} + {\bi F}_{h1-repel} + {\bi F}_{h2-repel}, \hfill &h = h_1 \hfill \cr {\bi F}_{50{\rm-}repel} + {\bi F}_{h0-repel} \hfill & h=1,2,3,4 \& h \neq h_1. \hfill \cr {\bi F}_{50{\rm-}repel} \hfill & h=5 \hfill}\right.$$

5. THE EFFECTIVE COLLISION AVOIDANCE CONTROL LAW AND THE DEFINITION OF THE COMPOSITE CONTROL LAW

If $x_{1-t}^{2}/ \sigma_{x-t}^{2} + y_{1-t}^{2}/ \sigma_{1-t}^{2} \lt D_{0}$, the relative parallel velocity $v_{1-t}^{paral} \gt 0$ at any time t and $\forall \vert {\bi r}_{1-t} \vert \neq\, 0$ are satisfied and the repulsive force is activated. In this section, we analyse the effectiveness of the avoidance scheme and design the composite law.

5.1. Analysis of the effective collision avoidance control law

According to Equations (9)–(56), the following conditions are satisfied:

1. 1) M _1ijd > M _2ijd > 0 (i = 1, 2, 3, 4, j = 0, 1, 2), M _150d > M _250d > 0.

2. 2) λ ₀ > 0, D _0ijd > R ₀ > 0 (i = 1, 2, 3, 4, j = 0, 1, 2), D _050d > R ₀ > 0.

3. 3) d ₀ > 0, a _max > 0, r_ij-td > 0 (i = 1, 2, 3, 4, j = 0, 1, 2), r _50-td > 0.

4. 4) v _ij-td^paral>0, v _ij-td^paral>0 (i = 1,2,3,4, j = 0,1,2), v _50-td^paral>0, v _50-td^paral>0.

As the parallel repulsive force is much larger than the perpendicular repulsive force (Ge and Cui, Reference Ge and Cui2002), we get

(57)$$\lambda_0 r_{ij-td} M_{1ij-td} \left(\displaystyle{1 \over \sigma_{x-t}^2} + \displaystyle{1 \over \sigma_{y-t}^2}\right)\gg \lambda_0 D_{0ijd} M_{2ijd} d_{0} \displaystyle{v_{ij-d}^{paral} \over a_{\max}} \displaystyle{v_{ij-d}^{perpen} \over r_{ij-d}} \gt 0 (i=1,2,3,4, j=0,1,2)$$

(58)$$\lambda_0 r_{50{\rm-}td} M_{150{\rm-}td} \left(\displaystyle{1 \over \sigma_{x-t}^2} + \displaystyle{1 \over \sigma_{y-t}^2}\right)\gg \lambda_0 D_{050d} M_{250{\rm-}td} d_{0} \displaystyle{v_{50{\rm-}d}^{paral} \over a_{\max}} \displaystyle{v_{50{\rm-}d}^{perpen} \over r_{50{\rm-}td}} \gt 0$$

Subsequently, from Equations (48) and (56), we obtain the total collision avoidance manoeuvres as follows:

(59)$${\bi F}_{Total-repel} = {\bi F}_{Total-paral} + {\bi F}_{Total-perpen}$$

where:

(60)$${\bi F}_{Total-paral} = \left\{\matrix{x_{ep} = 0 \,or\, y_{ep} = 0 \left\{\matrix{F_{50{\rm-}paral} + F_{h0-paral}, \hfill & h=1,2,3,4 \hfill \cr F_{50{\rm-}paral},\hfill & h =5 \hfill}\right. \hfill \cr x_{ep} \neq 0 \, and \, y_{ep} \neq 0 \left\{\matrix{F_{50{\rm-}paral} + F_{h1-paral} + F_{h2-paral}, \hfill & h=h_1\hfill \cr F_{50{\rm-}paral} + F_{h0-paral}, \hfill & h =1,2,3,4 \& h \neq h_1\hfill \cr F_{50{\rm-}paral}, \hfill & h=5 \hfill}\hfill \right.}\right.$$

(61)$${\bi F}_{Total-perpen} = \left\{\matrix{x_{ep} = 0 \, or \, y_{ep} = 0 \left\{\matrix{F_{50{\rm-}perpen} + F_{h0-perpen}, \hfill & h = 1,2,3,4 \hfill \cr F_{50{\rm-}perpen}, \hfill & h =5\hfill}\right. \hfill \cr x_{ep} \, \neq \, 0 \, and \, y_{ep} \, \neq \, 0 \left\{\matrix{F_{50{\rm-} perpen} + F_{h1- perpen} + F_{h2- perpen}, \hfill & h = h_1 \hfill \cr F_{50{\rm-} perpen} + F_{h0- perpen}, \hfill &h = 1,2,3,4 \, \& \, h \, \neq \, h_1 \hfill \cr F_{50{\rm-} perpen}, \hfill &h = 5 \hfill}\right.}\right.$$

Therefore, based on Equations (57) and (58), the parallel force F_{Total−paral} is along the negative direction of ${\bi n}_{1-t}^{paral}$ and prevents the deputy spacecraft from moving closer to the chief spacecraft. The perpendicular component F_{Total−operpen} of the force F_{Total−repel} steers the deputy spacecraft around the chief spacecraft.

As the parallel repulsive force is much larger than the perpendicular repulsive force (Ge and Cui, Reference Ge and Cui2002), therefore, the main analysis of the repulsive force is with respect to the parallel repulsive force. Moreover, as the force F_Total-repel contains the force F_50-paral in every condition, F_50-paral can be taken as an example to show the analysis progress in proving the effectiveness of the MECPC method. Subsequently, the parallel repulsive force generated by part 5 is differentiated with respect to the relative distance as:

(62)$$\displaystyle{\partial {\bi F}_{50{\rm-}paral} \over \partial r_{50{\rm-}td}} =-\lambda_0 M_{150d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2} \right)^2 r_{50{\rm-}td}^2 + \lambda_0 M_{150d} \left( \displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2}\right).$$

Similar to Equation (62), the auxiliary function I with respect to h (r _50-td) is defined as:

(63)$$h(r_{50{\rm-}td}) =-\lambda_0 M_{150d} \left( \displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2}\right)^2 r_{50{\rm-}td}^2 + \lambda_0 M_{150d} \left( \displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2}\right).$$

With the help of Equation (63), the defined auxiliary function is a second-order equation with respect to r _50-td and its quadratic coefficient is negative. Therefore, the maximum value of Equation (63) can be obtained as:

(64)$$H_{\max} (r_{50{\rm-}td}) = \lambda_0 M_{150d} \left(\displaystyle{1 \over \sigma_{x-td}^2} + \displaystyle{1 \over \sigma_{y-td}^2}\right) \gt 0.$$

Equation (63) has two zero points and its symmetry axis is zero. Moreover, we have r _50-td > 0. Then, x ₂>0 is assumed to be one zero point with respect to h(r _50-td), which is given as:

(65)$$X_2 = \displaystyle{\sigma_{x-td}\sigma_{y-td} \over \sqrt{\sigma_{x-td}^2 + \sigma_{y-td}^2}}.$$

As a result, if r _50-td > x ₂, h (r _50-td) < 0, if 0 < r _50-td < x ₂,   h (r _50-td) > 0.

From Equation (9), x ₃ is assumed to be the maximum value of r _50-td, which is defined as:

(66)$$r_{50{\rm-}td} \lt x_3 = \max \left(\sqrt{\sigma_{y-t}^2 D_{050} + x_{1-t}^2 \left(1 -\displaystyle{\sigma_{y-t}^2 \over \sigma_{x-t}^2}\right)}, \sqrt{\sigma_{y-t}^2 D_{050} + x_{1-t}^2 \left(1 -\displaystyle{\sigma_{y-t}^2 \over \sigma_{x-t}^2}\right)}\right).$$

According to Equations (62)–(66), if x ₂ < r _50-td < x ₃, ∂F_50-paral/∂r _50-td < 0. Moreover, F_50-paral is a decreasing function; if 0 < r _50-td < x ₂, ∂F_50-paral/∂r _50-td < 0. F_50-paral is an increasing function. Then, the critical point ζ₅₀ is the minimum relative distance between the deputy spacecraft and the centre of part 5 of the chief spacecraft and thus we obtain the critical point. When the deputy spacecraft reaches the critical point ζ₅₀, its parallel relative velocity decreases to zero.

According to the above analysis, the critical points ζ _ij (i = 1, 2, 3, 4, j = 0, 1, 2) and ζ₅₀ with respect to those parts can be obtained. Figure 4 shows the interval of speed range in the collision avoidance manoeuvre in two conditions. Using the available parameters’ values {λ ₀, d ₀}, the deputy spacecraft will not collide with the chief spacecraft when the following condition is satisfied:

Figure 4.

The interval of speed range in collision avoidance manoeuvre in two conditions.

(67)$$\left\{\matrix{\zeta_{50} \gt R_{50}\, \& \, \zeta_{i0} \gt R_{i0} (i=1,2,3,4); If\, x_{ep} = 0 \, or \, y_{ep} = 0\hfill \cr \zeta_{50} \gt R_{50}\, \& \, \zeta_{ij} \gt R_{ij} (i=1,2,3,4, j=1,2); If\, x_{ep} \neq \, 0 \, or \, y_{ep} \neq \, 0\hfill} \right.$$

where R _ij (i = 1, 2, 3, 4, j = 0, 1, 2) are the radii of the small components’ exterior envelope circle. R ₅₀ is the radius of part 5’s exterior envelope circle.

5.2. The composite control law

Considering a system as Equation (2), the LQR controller is utilised to find an optimal control law. Furthermore, the cost function (Lin, Reference Lin2007; Xing et al., Reference Xing, Yu, Wang, Zheng and Chen2016) is given as:

(68)$$J_0 = \displaystyle{1 \over 2} \int_{t_0}^{\infty} [{\bi X}^{\rm T} {\bi QX} + {\bi u}^{\rm T} {\bi Ru}]dt$$

where Q ≥ 0 and 1 ≥ R ≥ 0 are semi-positive definite symmetric matrices. According to the principle of minimum values as per Equation (68), we obtain the optimal control law as follows:

(69)$${\bi u}^{\ast} = {\bi K}_1 {\bi X}$$

where K₁ = −R⁻¹B^TS₁ is the optimal feedback gain matrix and S₁ is a matrix which satisfies the Riccati equation:

(70)$${\bi S}_1 {\bi A} + {\bi A}^{\rm T} {\bi S}_1 - {\bi S}_1 {\bi BR}^{-1} {\bi B}^{\rm T} {\bi S}_1 + {\bi Q} = 0.$$

According to Equations (2) and (68), we can obtain u* (t) by solving the Riccati algebraic equation. Then, the ILQR controller is proposed to track the reference trajectory and introduced below.

Based on Equation (59), we obtain:

(71)$${\bi F}_{Total-repel} = b_{11} {\bi n}_{1-td}^{paral} + b_{12} {\bi n}_{1-td}^{perpen}$$

where b ₁₁, b ₁₂ are the auxiliary parameters.

With the help of Equations (9)–(14), Equation (71) is computed by:

(72)$${\bi F}_{Total-repel} = {\bi K}_2 {\bi X}$$

where:

(73)$${\bi K}_2 = \displaystyle{1 \over r_{1-t}} \left[\matrix{-b_{11} & b_{12} & 0 & 0 \cr -b_{12} & b_{11} & 0 & 0} \right].$$

From Equation (73), we get:

(74)$${\bi K}_2^{\rm T} {\bi K}_2 \leq \left[\matrix{{\bi K}_3 & 0_{2\times 2}\cr 0_{2\times 2} &0_{2\times 2}} \right]$$

where:

(75)$${\bi K}_3 = \displaystyle{b_{11}^2 + b_{12}^2 \over r_{1-t}^2} {\bi I}_{2\times 2}$$

Furthermore, the upper matrix K₄ on the ${\bi K}_{2}^{\rm T} (t) {\bi K}_{2} (t)$ is defined as:

(76)$${\bi K}_4 = \left[\matrix{K_{upper} {\bi I}_{2\times 2} &0_{2\times 2}\cr 0_{2\times 2} &0_{2\times 2}} \right]$$

where:

(77)$$\displaystyle{b_{11}^2 + b_{12}^2 \over r^2} \leq K_{upper}$$

Considering a system as Equation (2), the ILQR controller aims to find an optimal control law which minimises the following cost function (Xing et al., Reference Xing, Yu, Wang, Zheng and Chen2016):

(78)$$J = \displaystyle{1 \over 2} \int_{t_0}^{\infty} \left[{\bi X}^{\rm T} \left(\displaystyle{{\bi K}_4 \over m^2} +{\bi Q} \right) {\bi X} + {\bi u}^{\rm T} {\bi Ru} \right]dt$$

where m is the mass of the deputy spacecraft. With the help of the principle of minimum values as per Equation (78), we get the optimal control as:

(79)$${\bi u}_1^{\ast} = {\bi K}_5 {\bi X}$$

where K₅ = −R⁻¹B^TS² is the optimal feedback gain matrix and S₂ is a matrix which satisfies the Riccati equation:

(80)$${\bi S}_2 {\bi A} + {\bi A}^{\rm T} {\bi S}_2 - {\bi S}_2 {\bi BR}^{-1} {\bi B}^{\rm T} {\bi S}_2 + \displaystyle{{\bi K}_4 \over m^2} + {\bi Q} = 0$$

Based on Equations (2) and (78), the design of ${\bi u}_{1}^{\ast}$ can be obtained by solving the Riccati algebraic equation. Then, the auxiliary function II (Lin, Reference Lin2007) is given as:

(81)$${\bi V}_1 ({\bi X}) = \min_{u^{\ast}} \int_0^{\infty} \left({\bi X}^{\rm T} \left(\displaystyle{{\bi K}_4 \over m^2} + {\bi Q} \right) {\bi X} + {\bi u}^{\rm T} {\bi Ru} \right) dt.$$

As the V₁ (X) is the minimum cost of the optimal control of the nominal system from some initial state X, V₁ (X) is defined as a Lyapunov function of Equation (2). By definition, V₁ (X) must satisfy the Hamilton-Jacobi-Bellman equation (Lin, Reference Lin2007):

(82)$$\min_{u^{\ast}} \left\{{\bi X}^{\rm T} \left(\displaystyle{{\bi K}_4 \over m^2} + {\bi Q} \right) {\bi X} + {\bi u}^{\rm T} {\bi Ru} + {\bi V}_{1-X}^{\rm T} ({\bi AX} + {\bi Bu})\right\} = 0$$

where V_1−X = (∂V/∂X). Since ${\bi u}_{1}^{\ast} = {\bi K}_{5} {\bi X}$ is the optimal control, it must satisfy: (1) the above minimum zero; and (2) the derivative of ${\bi X}^{\rm T} ({\bi K}_{4}/m^{2} +{\bi Q}) {\bi X} + {\bi u}^{\rm T} {\bi Ru} + {\bi V}_{1-X}^{\rm T} ({\bi AX} + {\bi Bu})$ with respect to u is zero.

(83)$${\bi X}^{\rm T} \left(\displaystyle{{\bi K}_4 \over m^2} + {\bi Q} \right) {\bi X} + {\bi X}^{\rm T} {\bi K}_5^{\rm T} {\bi RK}_5 {\bi X} + {\bi V}_{1-X}^{\rm T} ({\bi AX} + {\bi BK}_5 {\bi X}) = 0$$

(84)$$2{\bi X}^{\rm T} {\bi K}_5^{\rm T} + {\bi V}_{1-X}^{\rm T} {\bi B} = 0.$$

So, the composite controller is given as:

(85)$${\bi u}_{total} = {\bi u}_1^{\ast} + \displaystyle{{\bi F}_{Total-repel} \over m}$$

Based on the Lyapunov-based method, the stability of the overall close-system with the composite controller is verified (Lin, Reference Lin2007).

5.3. The flow of the MECPC method algorithm

With the aid of Equations (2)–(85), we obtain Algorithm 1 to describe the MECPC method. The control error of the deputy spacecraft is assumed to be an uncorrelated zero-mean Gaussian distribution (Psiaki, Reference Psiaki2011; Peynot et al., Reference Peynot, Lui, McAllister, Fitch and Sukkarieh2014). Furthermore, the covariance matrix of the control uncertainty is given as follows (Luo et al., Reference Luo, Liang, Wang and Tang2011):

(86)$$C_{\delta \Delta v_t} = \left\{\matrix{0, \hfill & if \, \vert {\bi u} \vert = 0 \hfill \cr Diag \left(\left[ (0{\cdot}02 + 0{\cdot}01 \vert u_x \vert)^2, (0{\cdot}02 + 0{\cdot}01 \vert u_y \vert)^2\right] \right), \hfill & if\, \vert {\bi u}\vert \neq \, 0}.\right.$$

Algorithm 1.

The MECPC method.

6. SIMULATION RESULTS AND ANALYSIS

To verify the effectiveness of the proposed method, a numerical simulation was carried out. Using the relative dynamic motion Equation (2), the relative motion of the deputy spacecraft was designed and propagated (Xing et al., Reference Xing, Tang, Xi and Li2007; Ou and Zhang, Reference Ou and Zhang2017b). Table 1 gives the physical parameters of the chief spacecraft and the deputy spacecraft. Table 2 lists the initial relative state vector of the deputy spacecraft in the LVLH frame. Figure 5 shows the reference relative motion of the deputy spacecraft. There are three blue asterisks, which are asterisks 1, 2, and 3. The deputy spacecraft begins at asterisk 1 with the initial relative velocity v₀ and flies around the chief spacecraft along the grey trajectory without any transfer impulse. Then, the deputy spacecraft arrives at asterisk 3 and the transfer impulse is implemented on the deputy spacecraft to reach asterisk 2. According to the pulse rendezvous for orbital transformation, we compute the control impulse as Table 3 shows. Table 3 lists the initial impulse control. The purple curves are the reference relative trajectories and the yellow curve is the boundary of the chief spacecraft's size. Moreover, the working area is defined as the area where the deputy spacecraft will execute further operations, such as space inspection and maintenance. A circle is centred on the planned terminal position asterisk 2 and has a certain radius. The working area is the overlap part between the circle and outside the yellow curve.

Figure 5.

The reference relative motion of the deputy spacecraft.

Table 1.

The physical parameters of the deputy spacecraft and chief spacecraft.

Table 2.

The initial relative state vector of the deputy spacecraft in the LVLH frame.

Table 3.

The initial impulse control in the LVLH frame.

The final condition is t _f=3, 400 s and the radius of the working area is 10 m. Both the frequency of the delta velocity manoeuvres and the integration step are 1 s. The semi-major axis and the eccentricity of the orbit of the chief spacecraft are 6,778.1336 km and 0, respectively, λ ₀ =30 and d ₀ =3. The distribution with respect to the navigation uncertainty is assumed as an uncorrelated zero-mean Gaussian (Psiaki, Reference Psiaki2011; Peynot et al., Reference Peynot, Lui, McAllister, Fitch and Sukkarieh2014). Moreover, the covariance matrix is given as follows (Luo et al., Reference Luo, Liang, Wang and Tang2011):

(87)$$C_{\delta X_0} = diag ([16, 16, 1{\cdot}6 \times 10^{-3}, 1{\cdot}6\times 10^{-3}])$$

In this simulation, Q = diag [10⁻¹⁴, 10⁻¹⁴, 0, 0];R = diag [10⁻⁸, 10⁻⁸] and the gravitational constant of Earth is 3 · 986 × 10¹⁴m ³/².

6.1. Case with ECPC method with λ₀=30

In this section, we utilise the traditional ECPC method to avoid collision. Figures 6(a) and 6(b) show the actual trajectory of the deputy spacecraft over 100 Monte Carlo runs with a LQR controller and ILQR controller, respectively. The yellow curve is the boundary of the chief spacecraft's size and the purple curve is the reference relative trajectory of the deputy spacecraft. The red curve is the boundary of the working area. Meanwhile, the green region is the actual trajectory of the deputy spacecraft over 100 Monte Carlo runs. According to Figure 6, when only the ECPC method is executed on the deputy spacecraft, the deputy spacecraft does collide with the chief spacecraft in the presence of a convex polygon shape.

Figure 6.

The trajectory of the deputy spacecraft over 100 Monte Carlo runs.

Figure 7 illustrates the values of collision probability over 100 Monte Carlo runs. The blue part is the value of collision probability and the light green part is the ratio of safe missions. From Figure 7, a collision between the deputy spacecraft and the chief spacecraft will occur. Furthermore, the collision probabilities with respect to the LQR controller and the ILQR controller are 31% and 33%. Thus, from Figures 6 and 7, we can conclude that the ECPC method does not solve the safe proximity problem in the presence of a convex polygon shape.

Figure 7.

The values of collision probability over 100 Monte Carlo runs.

For comparison, 100 independent Monte Caro runs were executed under the same conditions, and the Root Mean Square Error (RMSE) in position was calculated as a performance metric. The RMSE is given as (Chen et al., Reference Chen, Yin, Zhou, Wang, Wang and Chen2018):

(88)$$\hbox{RMSE} = \sqrt{Num^{-1} \sum_{t=1}^{Num} (z_1\hbox{--}z_2)^2}$$

where Num represents the total number of iterations. z ₁ and z ₂ represent two different series of data. In this paper, the two different series are the actual trajectory and the reference trajectory in the x and y directions.

Figure 8 depicts a comparison of RMSE behaviours in the x and y directions. Based on the ECPC method with the LQR controller, the brown curve and the pink curve represent the RMSE value in the x and y directions, respectively. Combining the ECPC method and the ILQR controller, we can obtain the blue curve and the green curve. The blue curve and the green curve are the RMSE values with respect to the x direction and y directions. From Figure 10, the RMSE values with the LQR controller are larger than those with the ILQR controller. Compared with the LQR controller, the ILQR controller obtains higher control accuracy in the influence of uncertainties. In other words, the robustness of the system is strengthened.

Figure 8.

Comparison with respect to the RMSE in two directions.

Figure 9 illustrates the associated velocity increment consumptions over 100 Monte Carlo runs. The green curve and the blue curve are the velocity increment consumptions of the deputy spacecraft with the ILQR controller and the LQR controller, respectively. The green curve fluctuates between 90 m/s and 95 m/s while the blue curve is around 89 m/s. From Figure 9, the deputy spacecraft with a ILQR controller requires more fuel. The reason is that K₄/m ² + Q in the ILQR controller is larger than Q in the LQR controller while the R in the two methods have equal values.

Figure 9.

The velocity increment consumptions over 100 Monte Carlo runs.

Figure 10 depicts the collision probability and the successful rate of arrival at the working area with different values of λ₀ over 100 Monte Carlo runs. The blue curve is the collision probability and the green curve represents the success rate. From Figure 10, if we adopt the ECPC method with different values of λ₀, collisions do also take place and the safe proximity problem in the presence of the convex polygon shape is not solved.

Figure 10.

The collision probability and the successful rate with different values of λ₀ over 100 Monte Carlo runs.

6.2. Case with MECPC method with λ₀=30

To solve the safe proximity problem in the presence of a convex polygon shape, the MECPC method is proposed to generate the avoidance manoeuvres. Figure 11 shows the actual trajectory of the deputy spacecraft over 100 Monte Carlo runs. The yellow curve and the purple curve represent the boundary of the chief spacecraft's size and the reference relative trajectory of the deputy spacecraft. The red curve is the boundary of the working area and the green zone is the actual trajectory of the deputy spacecraft over 100 Monte Carlo runs. From Figure 11, as the avoidance manoeuvres generated by the MECPC method are implemented, the deputy spacecraft successfully avoids collision with the chief spacecraft.

Figure 11.

The actual trajectory of deputy spacecraft over 100 Monte Carlo runs.

Figure 12 shows the comparison of RMSE behaviours in the x and y directions. According to the MECPC method with the LQR controller, we get the brown curve and the pink curve. The brown curve and the pink curve represent the RMSE values with respect to the x and y directions. The blue curve and the green curve are the RMSE values in the x and y directions, respectively. From Figure 12, the RMSE values with the ILQR controller are lower than the values with the LQR controller. Therefore, the ILQR controller obtains higher control accuracy.

Figure 12.

Comparison with respect to RMSE in two directions.

Figure 13 illustrates the associated velocity increment consumptions over 100 Monte Carlo runs. The blue curve and the green curve represent the velocity increment consumptions of the deputy spacecraft with the LQR controller and the ILQR controller, respectively. The green curve is about 94 m/s and the blue curve fluctuates around 90 m/s. Comparing the two curves in Figure 13, the deputy spacecraft with the ILQR controller requires more fuel.

Figure 13.

The velocity increment consumptions over 100 Monte Carlo runs.