# Time and Covariance Threshold Triggered Optimal Uncooperative Rendezvous Using Angles-Only Navigation

Authors: Yue You; Hua Wang; Christophe Paccolat; Volker Gass; Jean-Philippe Thiran; Jiu Ren Li
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

A time and covariance threshold triggered optimal maneuver planning method is proposed for orbital rendezvous using angles-only navigation (AON). In the context of Yamanaka-Ankersen orbital relative motion equations, the square root unscented Kalman filter (SRUKF) AON algorithm is developed to compute the relative state estimations from a low-volume/mass, power saving, and low-cost optical/infrared camera’s observations. Multi-impulsive Hill guidance law is employed in closed-loop linear covariance analysis model, based on which the quantitative relative position robustness and relative velocity robustness index are defined. By balancing fuel consumption, relative position robustness, and relative velocity robustness, we developed a time and covariance threshold triggered two-level optimal maneuver planning method, showing how these results correlate to past methods and missions and how they could potentially influence future ones. Numerical simulation proved that it is feasible to control the spacecraft with a two-line element- (TLE-) level uncertain, 34.6% of range, initial relative state to a 100 m v-bar relative station keeping point, at where the trajectory dispersion reduces to 3.5% of range, under a 30% data gap per revolution on account of the eclipse. Comparing with the traditional time triggered maneuver planning method, the final relative position accuracy is improved by one order and the relative trajectory robustness and collision probability are obviously improved and reduced, respectively.


---

## Body

## 1. Introduction

There are hundreds of thousands of debris in Earth’s orbit; while the biggest are tracked through ground radar, the large majority of pieces remain invisible. The 2009 collision between Iridium 33 and Kosmos-2251 destroyed an active satellite worth tens of millions of dollars [1]. There is international recognition that 5–10 pieces of the currently existing large debris should be removed each year by 2020. ESA’s CleanSpace program [2] is targeted at developing these capabilities, with the specific target of EnviSat, an 8-ton satellite that poses a high risk of a catastrophic collision. CleanSpace One (CSO) is EPFL’s hugely ambitious response to the problem of orbital debris [3]. The most capable nanosatellite ever built, CSO will autonomously intercept, capture, and deorbit a targeted object, EPFL’s own “Swiss Cube,” a 1 kg satellite launched in 2009. The disruptive technologies developed and demonstrated through this project will form the baseline for the upcoming active debris removal (ADR) efforts in Switzerland and abroad. The technologies will also serve for other uncooperative rendezvous mission, such as on-orbit servicing and inspection mission, potentially leading to the creation of autonomous orbital inspection and repair drones.The practice of single-optical/infrared camera based on AON provides a low-volume/mass, power saving, and low-cost solution for the nanosatellite in ADR. For these reasons, autonomous vehicles using AON are currently an active area of research [4–14]. The growing interest in vision-based autonomous rendezvous and docking has produced a series of experimental spacecraft in an attempt to develop rendezvous and proximity operations technology that would be more appropriate for small satellite. PRISMA [15] is the only in-orbit test bed of angles-only navigation, OHB Sweden [16], CNES [17], and DLR [18] imitated uncooperative rendezvous using AON separately, but the uncertainty is accuracy known with GPS and RF. Besides, all the rendezvous processes are ground-in-loop control; real AON based autonomous uncooperative rendezvous needs to be further studied.A fundamental challenge faced by AON is the inherent difficulty in accurately determining the range to target [7]. In rigorous terms, the navigation problem becomes observable only in the presence of orbit control maneuvers, which change the relative orbit geometry of the observation [11]. Woffinden and Geller [6, 7] and Grzymisch and Fichter [14] separately derived an observability criterion for AON and optimized a single-pulse maneuver to minimize estimation error in relative range direction. Spurmann [19] proposed the spiraling maneuver approach of AON to cover the transition gap from absolute to relative navigation. Jaehwan and Hyochoong [20] implemented the observability constraint at each rendezvous maneuver optimization node in which the cost function was aimed at minimizing fuel consumption; their method minimized fuel consumption while ensuring sufficient relative observability during the whole rendezvous maneuver. All the previous researches made observability as an optimization cost function or constraint, being offline and open loop [21]. But for uncooperative rendezvous especially long-time autonomous mission, the modeling error and control error exists, the adjacent maneuvers interval constraints, and data gap induced by solar eclipse must be taken into concern as well [22]. Besides observability, safety which is related to trajectory dispersion and fuel efficiency since satellites usually have a tight mass budget should be tradeoff at the meantime. Tang et al. [23] and Luo et al. [24–26] took the navigation and control error into consideration and defined the open-loop optimal robust rendezvous planning method, making tradeoff between fuel consumption, rendezvous time, and trajectory robustness. Li et al. [27] furtherly constructed closed-loop multiple objective optimization problem (MOOP) considering the position robustness, velocity robustness, and fuel. Both the previous maneuver planning methods employed the heavy computational burden physical planning and nondominated sorting genetic algorithm (NSGA-II), which is not applicable for onboard considering limit computational resource.In this study, the relative rendezvous trajectory is optimized for AON between chaser (active satellite or space robot) and target (uncooperative satellite or space debris). The main purpose of this paper is to design an online rendezvous maneuver planning method which provides better observability and robustness for the whole mission period. The second purpose is to improve fuel consumption efficiency. There is a tradeoff between the relative position robustness, relative velocity robustness, and fuel consumption.The remainder of the paper is organized as follows. Problem Formulation defines the problem by addressing the basics of relative dynamic and optical navigation. Closed-Loop Linear Covariance Analysis introduces the analytical linear covariance method to quantify the robustness measure for closed-loop relative trajectory dispersion. In Optimal Maneuver Planning, a rendezvous maneuver design strategy is proposed to minimize fuel consumption while seeking high position and velocity robustness. Using the proposed strategy, a numerical optimization technique is implemented to design the rendezvous trajectory. The designed trajectory is verified by performing Monte Carlo simulation. Conclusion presents concluding remarks for this study.


## 2. Problem Formulation

### 2.1. Dynamic Modeling

The origin of a rotating LVLH (Local Vertical Local Horizontal) reference frame (Ot-xyz), which is used for relative motion description, is collocated with the debris c.m. (c.m. is short for center of mass). The relative position and velocity of the chaser c.m. with respect to the target c.m. in the LVLH coordinates are denoted by r and v, respectively. The relative motion equations for general elliptical orbits are the well-known Tschauner-Hempel (TH) equations [28] whose homogeneous solution is known as the Yamanaka-Ankersen state transition matrix [29](1)
Φ
t
,
t
0
=
Φ
r
r
t
,
t
0
Φ
r
v
t
,
t
0
Φ
v
r
t
,
t
0
Φ
v
v
t
,
t
0
=
Φ
f
Φ
-
1
f
0
,where the respective expressions of Φ-1f0 and Φf are as follows:(2)
Φ
-
1
f
0
=
1
λ
2
-
3
s
k
+
e
2
k
2
c
-
2
e
0
-
s
k
+
1
k
0
0
-
3
e
+
c
k
-
s
0
-
c
k
+
1
k
+
e
0
0
3
k
-
λ
2
e
s
0
k
2
0
0
-
3
e
s
k
+
1
k
2
-
2
+
e
c
λ
2
-
e
s
k
+
1
k
0
0
0
0
0
0
λ
2
c
-
λ
2
s
0
0
0
0
λ
2
s
λ
2
c
,
Φ
f
=
s
c
2
-
3
e
s
I
0
0
0
s
′
c
′
-
3
e
s
′
I
+
s
k
2
0
0
0
c
1
+
1
k
-
s
1
+
1
k
-
3
k
2
I
1
0
0
-
2
s
e
-
2
c
-
3
1
-
2
e
s
I
0
0
0
0
0
0
0
c
s
0
0
0
0
-
s
cwith k=1+ecos⁡f, c=kcos⁡f, s=ksin⁡f, I=∫f0f1/k2df=μ2/h3t-t0, and λ=1-e2 for shorthand notation. In addition, e is the eccentricity of the target orbit, f is the true anomaly, μ is the gravitational constant of Earth, and h is the norm of target’s angular momentum. Here, c′ and s′ indicate the first derivatives with respect to f, respectively.These solutions can be written in discrete form for impulsive inputΔvk as(3)xt=Φt,t0x0+δtkΔvk,where x0=[r0,v0]T and xt=rt,vtT are chaser’s relative states at initial moment t0 and time t, respectively. δtk is the Dirac function.


### 2.2. Camera Observation Modeling

The camera measurement frame (Om-xmymzm) is assumed to be aligned with the focal plane of the camera. Its orientation with respect to the chaser body frame is assumed to be known and constant. The pixel location of the debris c.m. is used to form an line of sight (LOS) vector from the debris c.m. to the camera, which is expressed in the camera frame and denoted by l=xmymzmT, as shown in Figure 1. Because the transformation from LVLH to the camera measurement frame can be calculated using knowledge of inertial attitude, position, and velocity, the LOS measurement expressed in the LVLH frame can be used to formulate the angle/bearing-only measurement equation providing the measurement angles to the debris, which is explicitly written as follows:(4)Z=εθ=hx=arctan⁡xx2+y2+zρt-zxρt-z-zx2+y2arctan⁡-yx2+z2,where ε and θ are the respective elevation and azimuth measurement angles. ρt is the average radius of the debris circular motion.Figure 1
Relative observation geometry.The measurement sensitive matrixH can be further achieved by differentiating (4) with respect to state x; that is,(5)h11=xρt-zx2+y2x2+y2+ρt-z2-zx2+z2,h12=yρt-zx2+y2x2+y2+ρt-z2,h13=x2+y2x2+y2x2+y2+ρt-z2-xx2+z2,h21=xyx2+z2x2+y2+z2,h22=-x2+z2x2+z2x2+y2+z2,h23=yzx2+z2x2+y2+z2.


### 2.3. Angles-Only Navigation

SRUKF [30–35] is employed as measurements filtering scheme; it utilizes a deterministic “sampling” approach to calculate mean and covariance terms of the relative state. 2L+1 (L is the state dimension) sigma points are propagated through the true nonlinearity, without approximation, and then a weighted mean and covariance are taken. This approach results in approximations that are accurate to the third order (Taylor series expansion) for Gaussian inputs for all nonlinearities. For non-Gaussian inputs, approximations are accurate to at least the second order [31]. In contrast, the linearization approach of the EKF results only in first-order accuracy.The implementation for AON is deduced as follows.Initialization (6) x ^ 0 = E x 0 , S 0 = c h o l E x 0 - x ^ 0 x 0 - x ^ 0 T ,where chol· stands for the Cholesky decomposition.Fork∈1,…,∞, calculate the sigma points:(7)χk-1=x^k-1x^k-1+ηSkx^k-1-ηSk.Time Update (8) χ ^ k = f χ ^ k - 1 , u k - 1 + w k - 1 , x ^ k - = ∑ i = 0 2 L W i m χ ^ i , k , S ^ k - = q r W 1 c χ ^ 1 : 2 L , k - x ^ k - Q k - 1 , S ^ k - = c h o l u p d a t e S ^ k - , χ ^ 0 , k - x ^ k - , W 0 c , Z ^ k = h χ ^ k + v k - 1 , z ^ k = ∑ i = 0 2 L W i m Z ^ i , kwe use the shorthand notation qr· to denote a QR decomposition of a matrix and cholupdate· to denote Cholesky factor updating.Measurement Update (9) S z ^ k = q r W 1 c Z ^ 1 : 2 L , k - z ^ k R k - 1 , S z ^ k = c h o l u p d a t e S z ^ k , Z ^ 0 , k - z ^ k , W 0 c , P x ^ k z ^ k = ∑ i = 0 2 L W i c χ ^ i , k - x ^ k - Z ^ i , k - z ^ k T , K k = P x ^ k z ^ k / S z ^ k T S z ^ k , x ^ k = x ^ k - + K k z k - z ^ k , U = K k S z ^ k , S z ^ k = c h o l u p d a t e S z ^ k , U , - 1 ,where Wi is a set of scalar weights, Wim is the mean weight, and Wic is the covariance weight, which are calculated as follows:(10)W0m=λλ+L,W0c=λλ+L+1-α2+β,Wim=Wic=12L+λ,i=1,…,2L,where λ=Lα2-1 and η=L+λ. The constant α determines the spread of the sigma points around x^ and is usually set to 10-4≤α≤1. β is used to incorporate prior knowledge of the distribution of x (for Gaussian distributions, β=2 is optimal). Also note that we define the linear algebra operation of adding a column vector to a matrix, that is, A±u as the addition of the vector to each column of the matrix.Furthermore, unlike the EKF, no explicit derivatives (i.e., Jacobians or Hessians) need to be calculated. SRUKF employs three powerful linear algebra techniques—QR decomposition,Cholesky factor updating, andefficient least squares. It directly propagates the Cholesky factor, S (the lower triangular matrix given by SST=P), thereby avoiding the need to refactor at each time step. The computational complexity of UKF, SRUKF, and EKF is of the same order (i.e., OL3) for state estimation; however, SRUKF is approximately 20% faster than UKF and approximately 10% faster than EKF [33].


## 2.1. Dynamic Modeling

The origin of a rotating LVLH (Local Vertical Local Horizontal) reference frame (Ot-xyz), which is used for relative motion description, is collocated with the debris c.m. (c.m. is short for center of mass). The relative position and velocity of the chaser c.m. with respect to the target c.m. in the LVLH coordinates are denoted by r and v, respectively. The relative motion equations for general elliptical orbits are the well-known Tschauner-Hempel (TH) equations [28] whose homogeneous solution is known as the Yamanaka-Ankersen state transition matrix [29](1)
Φ
t
,
t
0
=
Φ
r
r
t
,
t
0
Φ
r
v
t
,
t
0
Φ
v
r
t
,
t
0
Φ
v
v
t
,
t
0
=
Φ
f
Φ
-
1
f
0
,where the respective expressions of Φ-1f0 and Φf are as follows:(2)
Φ
-
1
f
0
=
1
λ
2
-
3
s
k
+
e
2
k
2
c
-
2
e
0
-
s
k
+
1
k
0
0
-
3
e
+
c
k
-
s
0
-
c
k
+
1
k
+
e
0
0
3
k
-
λ
2
e
s
0
k
2
0
0
-
3
e
s
k
+
1
k
2
-
2
+
e
c
λ
2
-
e
s
k
+
1
k
0
0
0
0
0
0
λ
2
c
-
λ
2
s
0
0
0
0
λ
2
s
λ
2
c
,
Φ
f
=
s
c
2
-
3
e
s
I
0
0
0
s
′
c
′
-
3
e
s
′
I
+
s
k
2
0
0
0
c
1
+
1
k
-
s
1
+
1
k
-
3
k
2
I
1
0
0
-
2
s
e
-
2
c
-
3
1
-
2
e
s
I
0
0
0
0
0
0
0
c
s
0
0
0
0
-
s
cwith k=1+ecos⁡f, c=kcos⁡f, s=ksin⁡f, I=∫f0f1/k2df=μ2/h3t-t0, and λ=1-e2 for shorthand notation. In addition, e is the eccentricity of the target orbit, f is the true anomaly, μ is the gravitational constant of Earth, and h is the norm of target’s angular momentum. Here, c′ and s′ indicate the first derivatives with respect to f, respectively.These solutions can be written in discrete form for impulsive inputΔvk as(3)xt=Φt,t0x0+δtkΔvk,where x0=[r0,v0]T and xt=rt,vtT are chaser’s relative states at initial moment t0 and time t, respectively. δtk is the Dirac function.


## 2.2. Camera Observation Modeling

The camera measurement frame (Om-xmymzm) is assumed to be aligned with the focal plane of the camera. Its orientation with respect to the chaser body frame is assumed to be known and constant. The pixel location of the debris c.m. is used to form an line of sight (LOS) vector from the debris c.m. to the camera, which is expressed in the camera frame and denoted by l=xmymzmT, as shown in Figure 1. Because the transformation from LVLH to the camera measurement frame can be calculated using knowledge of inertial attitude, position, and velocity, the LOS measurement expressed in the LVLH frame can be used to formulate the angle/bearing-only measurement equation providing the measurement angles to the debris, which is explicitly written as follows:(4)Z=εθ=hx=arctan⁡xx2+y2+zρt-zxρt-z-zx2+y2arctan⁡-yx2+z2,where ε and θ are the respective elevation and azimuth measurement angles. ρt is the average radius of the debris circular motion.Figure 1
Relative observation geometry.The measurement sensitive matrixH can be further achieved by differentiating (4) with respect to state x; that is,(5)h11=xρt-zx2+y2x2+y2+ρt-z2-zx2+z2,h12=yρt-zx2+y2x2+y2+ρt-z2,h13=x2+y2x2+y2x2+y2+ρt-z2-xx2+z2,h21=xyx2+z2x2+y2+z2,h22=-x2+z2x2+z2x2+y2+z2,h23=yzx2+z2x2+y2+z2.


## 2.3. Angles-Only Navigation

SRUKF [30–35] is employed as measurements filtering scheme; it utilizes a deterministic “sampling” approach to calculate mean and covariance terms of the relative state. 2L+1 (L is the state dimension) sigma points are propagated through the true nonlinearity, without approximation, and then a weighted mean and covariance are taken. This approach results in approximations that are accurate to the third order (Taylor series expansion) for Gaussian inputs for all nonlinearities. For non-Gaussian inputs, approximations are accurate to at least the second order [31]. In contrast, the linearization approach of the EKF results only in first-order accuracy.The implementation for AON is deduced as follows.Initialization (6) x ^ 0 = E x 0 , S 0 = c h o l E x 0 - x ^ 0 x 0 - x ^ 0 T ,where chol· stands for the Cholesky decomposition.Fork∈1,…,∞, calculate the sigma points:(7)χk-1=x^k-1x^k-1+ηSkx^k-1-ηSk.Time Update (8) χ ^ k = f χ ^ k - 1 , u k - 1 + w k - 1 , x ^ k - = ∑ i = 0 2 L W i m χ ^ i , k , S ^ k - = q r W 1 c χ ^ 1 : 2 L , k - x ^ k - Q k - 1 , S ^ k - = c h o l u p d a t e S ^ k - , χ ^ 0 , k - x ^ k - , W 0 c , Z ^ k = h χ ^ k + v k - 1 , z ^ k = ∑ i = 0 2 L W i m Z ^ i , kwe use the shorthand notation qr· to denote a QR decomposition of a matrix and cholupdate· to denote Cholesky factor updating.Measurement Update (9) S z ^ k = q r W 1 c Z ^ 1 : 2 L , k - z ^ k R k - 1 , S z ^ k = c h o l u p d a t e S z ^ k , Z ^ 0 , k - z ^ k , W 0 c , P x ^ k z ^ k = ∑ i = 0 2 L W i c χ ^ i , k - x ^ k - Z ^ i , k - z ^ k T , K k = P x ^ k z ^ k / S z ^ k T S z ^ k , x ^ k = x ^ k - + K k z k - z ^ k , U = K k S z ^ k , S z ^ k = c h o l u p d a t e S z ^ k , U , - 1 ,where Wi is a set of scalar weights, Wim is the mean weight, and Wic is the covariance weight, which are calculated as follows:(10)W0m=λλ+L,W0c=λλ+L+1-α2+β,Wim=Wic=12L+λ,i=1,…,2L,where λ=Lα2-1 and η=L+λ. The constant α determines the spread of the sigma points around x^ and is usually set to 10-4≤α≤1. β is used to incorporate prior knowledge of the distribution of x (for Gaussian distributions, β=2 is optimal). Also note that we define the linear algebra operation of adding a column vector to a matrix, that is, A±u as the addition of the vector to each column of the matrix.Furthermore, unlike the EKF, no explicit derivatives (i.e., Jacobians or Hessians) need to be calculated. SRUKF employs three powerful linear algebra techniques—QR decomposition,Cholesky factor updating, andefficient least squares. It directly propagates the Cholesky factor, S (the lower triangular matrix given by SST=P), thereby avoiding the need to refactor at each time step. The computational complexity of UKF, SRUKF, and EKF is of the same order (i.e., OL3) for state estimation; however, SRUKF is approximately 20% faster than UKF and approximately 10% faster than EKF [33].


## 3. Closed-Loop Linear Covariance Analysis

Linear covariance analysis method [36, 37] is employed to deduce an analytical closed-loop control error analysis model in this paper. In order to give a clear mathematical description, thereal,estimated, andnominal relative states are noted as x, x^, and x-, respectively. We define an augmented error state vector X consisting of the true dispersions, δx=x-x-, and navigated dispersions, δx^=x^-x, as follows:(11)X=δxδx^.The covariance of X is Ω=EXXT and initialized as(12)
Ω
0
=
P
0
0
6
×
6
0
6
×
6
P
^
0
,where P0 is the a priori covariance matrix of true trajectory dispersion and P^0 is the a priori covariance matrix of navigated dispersions.The complete closed-loop control process of autonomous rendezvous includes the following three sections; they aretime propagate,measure update, andcontrol correct, which are detailedly deduced in this section.


### 3.1. Time Propagate

Fork∈1,…,∞, the augmented state vector and its covariance are propagated as follows:(13)Xk-=Ψk∣k-1Xk-1+,Ωk∣k-1=Ψk∣k-1Ωk-1+cΨk∣k-1T+Wk-1Qk-1Wk-1T,where Wk-1=I6×606×6 and Ψk∣k-1 is the augmented state transfer matrix, which is defined as(14)
Ψ
k
∣
k
-
1
=
Φ
t
k
,
t
k
-
1
0
6
×
6
0
6
×
6
Φ
t
k
,
t
k
-
1
.


### 3.2. Measure Update

The updating of augmented state vector and covariance due to observation is as follows:(15)Xk-=Xk∣k-1-+KkZk-HkXk∣k-1-,Ωk-c=CkΩk∣k-1CkT+EkRkEkT,where the Kalman gain, Kk, is achieved in (9).


### 3.3. Control Correct

Here the common used multi-impulsive Hill guidance law is employed for rendezvous maneuver calculation. Usually the total rendezvous time is given as constant; when the maneuver time array,T=t1t2⋯tfT, is settled, then the rendezvous maneuver strategy is established. At every maneuver moment, tk, the current estimated relative state x^k=r^kv^kT and desired final relative state x-f=r-fv-fT are submitted into the following function to calculate the maneuver vector:(16)Δv^k=Φrv-1tf,tkr-f-Φrrtf,tkr^k-v^k.To achieve high velocity control accuracy, the last control is applied at terminal; that is,(17)Δv^f=-v^f.Because of maneuver correction, the augmented state vector and the corresponding covariance are revised as follows:(18)Xk+=DΔv^kΔV^x-k-c+Dx-k-c+IXk-,Ωk+c=MkΩk-cMkT+NkPkδvNkT,where Pkδv is the maneuver execution error covariance matrix; the specific expressions of Mk and Nk are as follows:(19)
M
k
=
D
x
-
k
-
c
+
I
6
×
6
D
Δ
v
^
k
Δ
V
^
x
-
k
-
c
0
6
×
6
D
Δ
v
^
k
Δ
V
^
x
-
k
-
c
+
D
x
-
k
-
c
+
I
6
×
6
,
N
k
=
I
6
×
6
0
6
×
6
,where DΔv^k=∂d/∂vΔv^k, Dx-k-c=∂d/∂Xx-k-c, and ΔV^x-k-c=∂v^/∂Xx-k-c.Finally, we achieve the complete analytical closed-loop linear covariance analysis model, which is feasible for onboard use to analyze the trajectory dispersion and evaluate the navigation performance as well as online maneuver planning.


## 3.1. Time Propagate

Fork∈1,…,∞, the augmented state vector and its covariance are propagated as follows:(13)Xk-=Ψk∣k-1Xk-1+,Ωk∣k-1=Ψk∣k-1Ωk-1+cΨk∣k-1T+Wk-1Qk-1Wk-1T,where Wk-1=I6×606×6 and Ψk∣k-1 is the augmented state transfer matrix, which is defined as(14)
Ψ
k
∣
k
-
1
=
Φ
t
k
,
t
k
-
1
0
6
×
6
0
6
×
6
Φ
t
k
,
t
k
-
1
.


## 3.2. Measure Update

The updating of augmented state vector and covariance due to observation is as follows:(15)Xk-=Xk∣k-1-+KkZk-HkXk∣k-1-,Ωk-c=CkΩk∣k-1CkT+EkRkEkT,where the Kalman gain, Kk, is achieved in (9).


## 3.3. Control Correct

Here the common used multi-impulsive Hill guidance law is employed for rendezvous maneuver calculation. Usually the total rendezvous time is given as constant; when the maneuver time array,T=t1t2⋯tfT, is settled, then the rendezvous maneuver strategy is established. At every maneuver moment, tk, the current estimated relative state x^k=r^kv^kT and desired final relative state x-f=r-fv-fT are submitted into the following function to calculate the maneuver vector:(16)Δv^k=Φrv-1tf,tkr-f-Φrrtf,tkr^k-v^k.To achieve high velocity control accuracy, the last control is applied at terminal; that is,(17)Δv^f=-v^f.Because of maneuver correction, the augmented state vector and the corresponding covariance are revised as follows:(18)Xk+=DΔv^kΔV^x-k-c+Dx-k-c+IXk-,Ωk+c=MkΩk-cMkT+NkPkδvNkT,where Pkδv is the maneuver execution error covariance matrix; the specific expressions of Mk and Nk are as follows:(19)
M
k
=
D
x
-
k
-
c
+
I
6
×
6
D
Δ
v
^
k
Δ
V
^
x
-
k
-
c
0
6
×
6
D
Δ
v
^
k
Δ
V
^
x
-
k
-
c
+
D
x
-
k
-
c
+
I
6
×
6
,
N
k
=
I
6
×
6
0
6
×
6
,where DΔv^k=∂d/∂vΔv^k, Dx-k-c=∂d/∂Xx-k-c, and ΔV^x-k-c=∂v^/∂Xx-k-c.Finally, we achieve the complete analytical closed-loop linear covariance analysis model, which is feasible for onboard use to analyze the trajectory dispersion and evaluate the navigation performance as well as online maneuver planning.


## 4. Optimal Maneuver Planning

### 4.1. Optimization Fitness Function

The AON rendezvous maneuver planning is a typical multiobjective optimization problem; to make a tradeoff among fuel consumption, position control error, and velocity control error, the multiobjective optimization fitness functions should be established properly. Using the closed-loop linear covariance built in Section3, we achieve the extended error covariance at time tk:(20)Ωk=EδXkδXkTEδXkδX^kTEδX^kδXkTEδX^kδX^kT=PkδXδXPkδXδX^PkδX^δXPkδX^δX^,where PkδXδX is the trajectory dispersion covariance at time tk, which can be expressed as(21)PkδXδX=PrrPrvPvrPvv,where Prr is the relative position dispersion and represents the ellipsoid of trajectory dispersion in geometry. Its eigenvalues σx, σy, and σz represent the three semimajor axes of the ellipsoid, respectively. We take the geometric sum of eigenvalues as the robustness metric of relative position dispersion(22)σr=σx2+σy2+σz2.Similarly, we defined the robustness metric of relative velocity dispersion as(23)σv=σx˙2+σy˙2+σz˙2.Finally, the multiobjective optimization (MOO) fitness functions are achieved:(24)minf1Δv=wΔv∑i=1nΔv^iminf2σr=wσr∑k=0Ntk-t0tf-t0ησrkminf3σv=wσv∑k=0Ntk-t0tf-t0ησvk.The weight factors, wΔv, wσr, and wσv, are used to maintain three objectives at same quantity order, which improves the optimization effectiveness. As the reasonable quantities of ∑i=1nΔv^i, ∑k=0Ntk-t0/tf-t0ησrk and ∑k=0Ntk-t0/tf-t0ησvk are O10, O1000, and O1, respectively, ∑i=1nΔv^i is recommended to be selected as the reference value; as a result, wΔv=1, wσr=0.01, and wσv=10. η is an integer, which tunes the weight of errors at different time; larger η means bigger weight is addressed on the final part of time; it improves the dispersion of trajectories, that is, improving the trajectory robustness and safety. N is a positive integer, which tunes the overall trajectory dispersion; larger N means bigger weight is addressed on the trajectory dispersion of the whole mission period.The key performance required for onboard optimal maneuver planning is proper computational complexity; the feasible solution is to transfer the MOOP to a single-objective optimization problem (SOOP) by making sum of the three fitness functions:(25)f=min∑i=13fiif∑Δv∑Δv^≥ξ,f=wf.It is important to minimize the fuel consumption while ensuring sufficient robustness since a small-size satellite has a tight mass budget. We use a penalty factor, w, to address attention to the fuel consumption in fitness function, when the true fuel consumption departs further than the threshold, 1<ξ<2; from the nominal value, the penalty factor w will be multiplied. In order to make this punishment effective, w should be a positive number at least one order bigger than 1; here w=10 is recommended. In this paper, the sequence quadratic program (SQP) algorithm is used for optimization calculation.


### 4.2. Time Triggered Maneuver Planning

Generally, the total rendezvous duration,tf, is given; when maneuver time array is determined the rendezvous scheme is established; that is, maneuvers are triggered by time. At each maneuver moment, the current estimated relative state is seen as the initial relative state, aiming at the desired relative station, multi-impulsive Hill guidance law is usually used. Regularly, the last maneuver is applied at the terminal, in order to confirm high relative velocity accuracy.Spacecraft circling the Earth on a low earth orbit (LEO) experiences a 30% data gap per revolution on account of the eclipse. For long duration (longer than one orbit period) mission, the eclipse should be taken into concerned especially for optical/infrared camera based AON scenarios. The eclipse can be located at any relative time since the rendezvous beginning time by phasing, so here we assume there arem observable interval during the whole rendezvous process; that is, 0,T-0.3T,T,2T-0.3T,…,(m-1)T,mT-0.3T; all the maneuvers are applied in observable interval, and there is ni (i∈1,2,…,m) maneuvers in the ith observable interval. Based on the previous analysis, the general time triggered maneuver profile planning model is mathematically defined as follows:(26)t=t1,1t1,2⋯t1,n1⋮⋮⋮⋮tm,1tm,2⋯tm,nm,tm,1≥m-1T+τ,tm,nm≤mT-0.3T,tm,k+1-tm,k≥τ,k∈1,2,…,nm-1,where T is the orbit period of target. n1,…,nm,m∈N, m is the revolution number since approaching started (m=1). The eclipse occurs between tm,nm and tm+1,1. To confirm enough information has been accumulated before maneuver, the first maneuver is arranged to τ minutes after the navigation started. To confirm all the maneuvers are under illuminated, the total maneuver duration should be shorter than (m-0.3)T (consecutive illuminated duration). To confirm that enough time has been left for attitude adjustment and maneuver preparation, the time interval between the adjacent maneuvers is longer than τ.Theoretically, long duration maneuver is more fuel efficient, but under big uncertainty as angles-only navigation, the trajectory disperses very fast, so real fuel consumption spreads over a large area around nominal value. More importantly, the collision probability rises sharply as the relative range decreases. So, the time triggered general handover maneuver planning method is regularly feasible for ground-in-loop mission, short duration mission, or well-equipped autonomous mission, in which cases the relative state is accurately estimated. As a result, the conflictive fuel consumption, control accuracy, and safety challenge angles-only navigation.


### 4.3. Covariance Threshold Triggered Maneuver Planning

In fact, the primary cause that time triggered maneuver planning method cannot apply to long duration autonomous rendezvous of AON is that time triggered planning method generates rendezvous scheme in offline mode; the real trajectory dispersion and navigation error are not used timely.In order to overcome the drawback of time triggered maneuver planning method mentioned previously, we developed a novel and smarter onboard autonomous maneuver planning method, as shown in Figure2. This method is divided into two steps, or two levels. Firstly, use the time triggered optimal maneuver planning method to achieve a basic optimal maneuver profile, t=t1t2⋯tN. Secondly, use the closed-loop linear covariance predication to judge whether midway amending maneuver should be inserted. If the current control error is bigger than the threshold, which can be changed according to the mission requirement of robustness, an amending maneuver will be inserted into the middle of the latest applied maneuver and the following maneuver; the trigger conditions are as follows:(27)∑i=13Ωki,i>σ-r,N≥Nmin,Δtk∣k+1≥2τ,where σ-r is the desired threshold of relative position control error relative to range; according to the one-percent rule [38] which is widely adopted in optical navigation in space relative position estimation accuracy should be maintained at less than one percent of the relative range between satellites. Obey the one-percent rule, and take the 3σ situation; it is suggested to set σ-r=3%. Nmin is the appointed already applied maneuver number; if it is too small, the maneuver correction effect is weak, and if it is too big, the amending effect is weak; according to rule of thumb, Nmin should be around half of total nominal maneuver number. τ is the interval constraint of adjacent maneuvers, left for filtering, maneuver preparation, and attitude adjustment; according to rule of thumb, τ is usually several minutes.Figure 2
Chart flow of time and covariance threshold triggered maneuver planning.Transfer time, fuel consumption, passive safety, and other constraints can be taken into consideration in the basic maneuver profile, using the time triggered optimal maneuver planning layer; however it is beyond the scope of this paper. This paper mainly gives midway amending strategy. We use the covariance threshold triggered method, that is, a positive safety method, to improve the trajectory dispersion and safety. Every time when the amending maneuver is applied, the left basic maneuver series will not be reoptimized again, although it will be better, but the onboard computer is hard to bear in current mission; when one day the onboard computer becomes powerful enough, this problem will be conquered.


### 4.4. Numerical Simulation

To validate the maneuver planning methods, the SwissCube TLE data from March 23, 2016, at UTC 13:45:16 was chosen to serve as the initial inertia state of the target, that is,a0 = 7086121.3370 m, e0 = 0.0007, i0 = 1.7176 rad, Ω0 = 5.0423 rad, ω0 = 5.9745 rad, and f0 = 0.31075 rad. Both the time and covariance triggered maneuver planning method and time triggered maneuver planning method are separately used to guide the chaser from the initial 10 km relative station keeping point to the final 100 m relative station keeping point. Ten waypoints are located on the v-bar; they are 8 km, 7.5 km, 6 km, 5 km, 4 km, 3 km, 2.6 km, 2 km, 1 km, and 500 m, respectively, ahead of the debris. The parameters of the optimization algorithm are selected as wΔv = 1, wσr = 0.01, wσv = 10, w = 10, Nmin = 2, σ-r = 3%, ξ = 1.2, and τ = 6 min.Uncooperative means no knowledge and/or data exchange between the chaser and the target. But artificial objects orbiting Earth whose diameter is bigger than 10 cm can be tracked by United States Space Surveillance Network of NORAD (North American Aerospace Defense Command), and its TLE data is public, available from the website. Its TLE data can serve as the prior information of debris’ inertial orbit parameter, and its uncertainty is at about kilometer level [39]. In the meantime, chaser’s inertial orbit parameter can be achieved through the onboard GPS or tracked by the ground tracking station; the accuracy is high. Using the uncertain debris’ inertial orbit parameter and accurate chaser’s inertial orbit parameter, the estimation of initial relative state between the chaser and debris is achieved. As analyzed in Stefan’s Master Dissertation, the 3σ relative state estimation uncertain is around kilometer level; here the initial relative state navigation covariance matrix is selected as P^0 = diag[(1800 m)2, (1200 m)2, (1200 m)2, (1.8 m/s)2, (1.2 m/s)2, (1.2 m/s)2], the TLE-accuracy-level initial relative position dispersion covariance matrix is P0 = diag[(2000 m)2, (2000 m)2, (2000 m)2, (1.0 m/s)2, (1.0 m/s)2, (1.0 m/s)2], and measure noise covariance matrix is R = diag[(0.003 rad)2, (0.003 rad)2]. For the system noise covariance matrix Q, refer to Eq. (41) in [40] with T=10 s. The measure frequency is set to be 0.1 Hz and a 30% data gap per period is considered because of solar eclipse.The observation is only carried out when the target is illuminated (assuming we are working with visible light). Figure3 shows relative position estimated error variation with 30% eclipse per orbital period, that is, every last 0.3 period is eclipse, for example, 2.7 to 3. In eclipse, the local vertical and horizontal position estimation error increase obviously. The orbital manoeuvers and the orbit calculations are however independent of whether we are in eclipse or not and can be carried out independently. But in our paper, in order to ensure high control accuracy, we assumed the orbital manoeuvers are carried out in illumination.Figure 3
Relative position estimated error variation with 30% eclipse in one-period Monte Carlo shooting.In the nominal PRISMA-like long-time v-bar hopping approach motion, four impulsive maneuvers are planned in each hopping loop. The total flight duration is 11 SwissCube orbital periods; its orbital period isP=2πa03/μ=5936.17 s, so the total flight duration is tf=11P=18.14 hours, and the nominal fuel consumption is 3.477 m/s. The Monte Carlo shooting results using different optimized maneuver profiles are shown in Table 1, and the sampled Monte Carlo shooting trajectories using different optical maneuver profiles are shown in Figures 4(a) and 4(b).Table 1
Monte Carlo shooting results under TLE level accuracy initial relative state uncertainty.
Strategy 3 σ r f (m) 3 σ v f (m/s) Maneuver Number ∑ Δ v / N (m/s) ∑ Δ v Std. (m/s) SQP TT 16.554 0.0061 44 7.291 2.879 SQP CT 3.477 0.0057 68 8.025 2.395 “TT” is the shorthand notation of time triggered optimal maneuver planning method, and “CT” is the shorthand notation of time and covariance threshold triggered optimal maneuver planning method.Figure 4
Monte Carlo shooting.
(a)
Time triggered optical maneuver profile (b)
Time and covariance threshold triggered optical maneuver profileIn the figures, the red solid line denotes the nominal relative trajectory, the gray “hair plots” are 100 runs Monte Carlo shootings, the black solid ellipses stand for3σ error ellipse of the corresponding nominal position, and blue dots are the real sampled positions at the meantime.In Table1, the time triggered optimal maneuver profile consists of 44 impulsive maneuvers—that is, four maneuvers per hopping loop. The time and covariance threshold triggered optimal maneuver profile consists of 68 impulsive maneuvers. The maneuver number of each hopping loop is [8,8,4,4,4,4,5,7,8,8,8], respectively. Nevertheless, the covariance triggered maneuvers cost 0.7340 m/s more fuel, the final relative control accuracy is improved by one order when compared with the time triggered maneuvers, and the trajectory robustness is improved obviously. More importantly, the covariance triggered trajectories all maintain a safe distance from the debris (original point); however, the debris is enveloped by the time triggered trajectories, which means that the chaser may collide with the debris, which is forbidden. Moreover, the semimajor axis of position control error ellipses lies near the sight direction, which, owing to the drawback of AON, cannot provide relative range information.In the current manuscript, the sun angle constraint with respect to the optical sensor field of view does not have to be considered. However, when the sun angle constraint is considered, it is equivalent to extend the solar eclipse duration, but how is the algorithms performance during eclipse? For low earth orbit spacecraft, the eclipse rate of one orbital period varies from 0% (when sun light is vertical to the orbital plane) to about 40% (when sun light is coplanar with the orbital plane). In one orbital period, when eclipse increases and the total maneuver duration decreases; refer to rule of thumb; the total fuel consumption will increase, which is verified through Monte Carlo shooting and shown in Figure5.Figure 5
Eclipse rate influence on mean fuel consumption.Figure6 gives the statistical mean of final relative position under various eclipse rate illumination conditions based on 100-time Monte Carlo shooting. For final relative position control accuracy, due to the maneuver correction, when eclipse rate is small (eclipse rate < 0.3), there is no obvious difference; however if the eclipse is big, the accuracy will be worse because of poor navigation accuracy. From another point of view, 30 percent eclipse in one orbital period is acceptable for visible camera based angles-only navigation.Figure 6
Eclipse rate influence on final relative position control error.


## 4.1. Optimization Fitness Function

The AON rendezvous maneuver planning is a typical multiobjective optimization problem; to make a tradeoff among fuel consumption, position control error, and velocity control error, the multiobjective optimization fitness functions should be established properly. Using the closed-loop linear covariance built in Section3, we achieve the extended error covariance at time tk:(20)Ωk=EδXkδXkTEδXkδX^kTEδX^kδXkTEδX^kδX^kT=PkδXδXPkδXδX^PkδX^δXPkδX^δX^,where PkδXδX is the trajectory dispersion covariance at time tk, which can be expressed as(21)PkδXδX=PrrPrvPvrPvv,where Prr is the relative position dispersion and represents the ellipsoid of trajectory dispersion in geometry. Its eigenvalues σx, σy, and σz represent the three semimajor axes of the ellipsoid, respectively. We take the geometric sum of eigenvalues as the robustness metric of relative position dispersion(22)σr=σx2+σy2+σz2.Similarly, we defined the robustness metric of relative velocity dispersion as(23)σv=σx˙2+σy˙2+σz˙2.Finally, the multiobjective optimization (MOO) fitness functions are achieved:(24)minf1Δv=wΔv∑i=1nΔv^iminf2σr=wσr∑k=0Ntk-t0tf-t0ησrkminf3σv=wσv∑k=0Ntk-t0tf-t0ησvk.The weight factors, wΔv, wσr, and wσv, are used to maintain three objectives at same quantity order, which improves the optimization effectiveness. As the reasonable quantities of ∑i=1nΔv^i, ∑k=0Ntk-t0/tf-t0ησrk and ∑k=0Ntk-t0/tf-t0ησvk are O10, O1000, and O1, respectively, ∑i=1nΔv^i is recommended to be selected as the reference value; as a result, wΔv=1, wσr=0.01, and wσv=10. η is an integer, which tunes the weight of errors at different time; larger η means bigger weight is addressed on the final part of time; it improves the dispersion of trajectories, that is, improving the trajectory robustness and safety. N is a positive integer, which tunes the overall trajectory dispersion; larger N means bigger weight is addressed on the trajectory dispersion of the whole mission period.The key performance required for onboard optimal maneuver planning is proper computational complexity; the feasible solution is to transfer the MOOP to a single-objective optimization problem (SOOP) by making sum of the three fitness functions:(25)f=min∑i=13fiif∑Δv∑Δv^≥ξ,f=wf.It is important to minimize the fuel consumption while ensuring sufficient robustness since a small-size satellite has a tight mass budget. We use a penalty factor, w, to address attention to the fuel consumption in fitness function, when the true fuel consumption departs further than the threshold, 1<ξ<2; from the nominal value, the penalty factor w will be multiplied. In order to make this punishment effective, w should be a positive number at least one order bigger than 1; here w=10 is recommended. In this paper, the sequence quadratic program (SQP) algorithm is used for optimization calculation.


## 4.2. Time Triggered Maneuver Planning

Generally, the total rendezvous duration,tf, is given; when maneuver time array is determined the rendezvous scheme is established; that is, maneuvers are triggered by time. At each maneuver moment, the current estimated relative state is seen as the initial relative state, aiming at the desired relative station, multi-impulsive Hill guidance law is usually used. Regularly, the last maneuver is applied at the terminal, in order to confirm high relative velocity accuracy.Spacecraft circling the Earth on a low earth orbit (LEO) experiences a 30% data gap per revolution on account of the eclipse. For long duration (longer than one orbit period) mission, the eclipse should be taken into concerned especially for optical/infrared camera based AON scenarios. The eclipse can be located at any relative time since the rendezvous beginning time by phasing, so here we assume there arem observable interval during the whole rendezvous process; that is, 0,T-0.3T,T,2T-0.3T,…,(m-1)T,mT-0.3T; all the maneuvers are applied in observable interval, and there is ni (i∈1,2,…,m) maneuvers in the ith observable interval. Based on the previous analysis, the general time triggered maneuver profile planning model is mathematically defined as follows:(26)t=t1,1t1,2⋯t1,n1⋮⋮⋮⋮tm,1tm,2⋯tm,nm,tm,1≥m-1T+τ,tm,nm≤mT-0.3T,tm,k+1-tm,k≥τ,k∈1,2,…,nm-1,where T is the orbit period of target. n1,…,nm,m∈N, m is the revolution number since approaching started (m=1). The eclipse occurs between tm,nm and tm+1,1. To confirm enough information has been accumulated before maneuver, the first maneuver is arranged to τ minutes after the navigation started. To confirm all the maneuvers are under illuminated, the total maneuver duration should be shorter than (m-0.3)T (consecutive illuminated duration). To confirm that enough time has been left for attitude adjustment and maneuver preparation, the time interval between the adjacent maneuvers is longer than τ.Theoretically, long duration maneuver is more fuel efficient, but under big uncertainty as angles-only navigation, the trajectory disperses very fast, so real fuel consumption spreads over a large area around nominal value. More importantly, the collision probability rises sharply as the relative range decreases. So, the time triggered general handover maneuver planning method is regularly feasible for ground-in-loop mission, short duration mission, or well-equipped autonomous mission, in which cases the relative state is accurately estimated. As a result, the conflictive fuel consumption, control accuracy, and safety challenge angles-only navigation.


## 4.3. Covariance Threshold Triggered Maneuver Planning

In fact, the primary cause that time triggered maneuver planning method cannot apply to long duration autonomous rendezvous of AON is that time triggered planning method generates rendezvous scheme in offline mode; the real trajectory dispersion and navigation error are not used timely.In order to overcome the drawback of time triggered maneuver planning method mentioned previously, we developed a novel and smarter onboard autonomous maneuver planning method, as shown in Figure2. This method is divided into two steps, or two levels. Firstly, use the time triggered optimal maneuver planning method to achieve a basic optimal maneuver profile, t=t1t2⋯tN. Secondly, use the closed-loop linear covariance predication to judge whether midway amending maneuver should be inserted. If the current control error is bigger than the threshold, which can be changed according to the mission requirement of robustness, an amending maneuver will be inserted into the middle of the latest applied maneuver and the following maneuver; the trigger conditions are as follows:(27)∑i=13Ωki,i>σ-r,N≥Nmin,Δtk∣k+1≥2τ,where σ-r is the desired threshold of relative position control error relative to range; according to the one-percent rule [38] which is widely adopted in optical navigation in space relative position estimation accuracy should be maintained at less than one percent of the relative range between satellites. Obey the one-percent rule, and take the 3σ situation; it is suggested to set σ-r=3%. Nmin is the appointed already applied maneuver number; if it is too small, the maneuver correction effect is weak, and if it is too big, the amending effect is weak; according to rule of thumb, Nmin should be around half of total nominal maneuver number. τ is the interval constraint of adjacent maneuvers, left for filtering, maneuver preparation, and attitude adjustment; according to rule of thumb, τ is usually several minutes.Figure 2
Chart flow of time and covariance threshold triggered maneuver planning.Transfer time, fuel consumption, passive safety, and other constraints can be taken into consideration in the basic maneuver profile, using the time triggered optimal maneuver planning layer; however it is beyond the scope of this paper. This paper mainly gives midway amending strategy. We use the covariance threshold triggered method, that is, a positive safety method, to improve the trajectory dispersion and safety. Every time when the amending maneuver is applied, the left basic maneuver series will not be reoptimized again, although it will be better, but the onboard computer is hard to bear in current mission; when one day the onboard computer becomes powerful enough, this problem will be conquered.


## 4.4. Numerical Simulation

To validate the maneuver planning methods, the SwissCube TLE data from March 23, 2016, at UTC 13:45:16 was chosen to serve as the initial inertia state of the target, that is,a0 = 7086121.3370 m, e0 = 0.0007, i0 = 1.7176 rad, Ω0 = 5.0423 rad, ω0 = 5.9745 rad, and f0 = 0.31075 rad. Both the time and covariance triggered maneuver planning method and time triggered maneuver planning method are separately used to guide the chaser from the initial 10 km relative station keeping point to the final 100 m relative station keeping point. Ten waypoints are located on the v-bar; they are 8 km, 7.5 km, 6 km, 5 km, 4 km, 3 km, 2.6 km, 2 km, 1 km, and 500 m, respectively, ahead of the debris. The parameters of the optimization algorithm are selected as wΔv = 1, wσr = 0.01, wσv = 10, w = 10, Nmin = 2, σ-r = 3%, ξ = 1.2, and τ = 6 min.Uncooperative means no knowledge and/or data exchange between the chaser and the target. But artificial objects orbiting Earth whose diameter is bigger than 10 cm can be tracked by United States Space Surveillance Network of NORAD (North American Aerospace Defense Command), and its TLE data is public, available from the website. Its TLE data can serve as the prior information of debris’ inertial orbit parameter, and its uncertainty is at about kilometer level [39]. In the meantime, chaser’s inertial orbit parameter can be achieved through the onboard GPS or tracked by the ground tracking station; the accuracy is high. Using the uncertain debris’ inertial orbit parameter and accurate chaser’s inertial orbit parameter, the estimation of initial relative state between the chaser and debris is achieved. As analyzed in Stefan’s Master Dissertation, the 3σ relative state estimation uncertain is around kilometer level; here the initial relative state navigation covariance matrix is selected as P^0 = diag[(1800 m)2, (1200 m)2, (1200 m)2, (1.8 m/s)2, (1.2 m/s)2, (1.2 m/s)2], the TLE-accuracy-level initial relative position dispersion covariance matrix is P0 = diag[(2000 m)2, (2000 m)2, (2000 m)2, (1.0 m/s)2, (1.0 m/s)2, (1.0 m/s)2], and measure noise covariance matrix is R = diag[(0.003 rad)2, (0.003 rad)2]. For the system noise covariance matrix Q, refer to Eq. (41) in [40] with T=10 s. The measure frequency is set to be 0.1 Hz and a 30% data gap per period is considered because of solar eclipse.The observation is only carried out when the target is illuminated (assuming we are working with visible light). Figure3 shows relative position estimated error variation with 30% eclipse per orbital period, that is, every last 0.3 period is eclipse, for example, 2.7 to 3. In eclipse, the local vertical and horizontal position estimation error increase obviously. The orbital manoeuvers and the orbit calculations are however independent of whether we are in eclipse or not and can be carried out independently. But in our paper, in order to ensure high control accuracy, we assumed the orbital manoeuvers are carried out in illumination.Figure 3
Relative position estimated error variation with 30% eclipse in one-period Monte Carlo shooting.In the nominal PRISMA-like long-time v-bar hopping approach motion, four impulsive maneuvers are planned in each hopping loop. The total flight duration is 11 SwissCube orbital periods; its orbital period isP=2πa03/μ=5936.17 s, so the total flight duration is tf=11P=18.14 hours, and the nominal fuel consumption is 3.477 m/s. The Monte Carlo shooting results using different optimized maneuver profiles are shown in Table 1, and the sampled Monte Carlo shooting trajectories using different optical maneuver profiles are shown in Figures 4(a) and 4(b).Table 1
Monte Carlo shooting results under TLE level accuracy initial relative state uncertainty.
Strategy 3 σ r f (m) 3 σ v f (m/s) Maneuver Number ∑ Δ v / N (m/s) ∑ Δ v Std. (m/s) SQP TT 16.554 0.0061 44 7.291 2.879 SQP CT 3.477 0.0057 68 8.025 2.395 “TT” is the shorthand notation of time triggered optimal maneuver planning method, and “CT” is the shorthand notation of time and covariance threshold triggered optimal maneuver planning method.Figure 4
Monte Carlo shooting.
(a)
Time triggered optical maneuver profile (b)
Time and covariance threshold triggered optical maneuver profileIn the figures, the red solid line denotes the nominal relative trajectory, the gray “hair plots” are 100 runs Monte Carlo shootings, the black solid ellipses stand for3σ error ellipse of the corresponding nominal position, and blue dots are the real sampled positions at the meantime.In Table1, the time triggered optimal maneuver profile consists of 44 impulsive maneuvers—that is, four maneuvers per hopping loop. The time and covariance threshold triggered optimal maneuver profile consists of 68 impulsive maneuvers. The maneuver number of each hopping loop is [8,8,4,4,4,4,5,7,8,8,8], respectively. Nevertheless, the covariance triggered maneuvers cost 0.7340 m/s more fuel, the final relative control accuracy is improved by one order when compared with the time triggered maneuvers, and the trajectory robustness is improved obviously. More importantly, the covariance triggered trajectories all maintain a safe distance from the debris (original point); however, the debris is enveloped by the time triggered trajectories, which means that the chaser may collide with the debris, which is forbidden. Moreover, the semimajor axis of position control error ellipses lies near the sight direction, which, owing to the drawback of AON, cannot provide relative range information.In the current manuscript, the sun angle constraint with respect to the optical sensor field of view does not have to be considered. However, when the sun angle constraint is considered, it is equivalent to extend the solar eclipse duration, but how is the algorithms performance during eclipse? For low earth orbit spacecraft, the eclipse rate of one orbital period varies from 0% (when sun light is vertical to the orbital plane) to about 40% (when sun light is coplanar with the orbital plane). In one orbital period, when eclipse increases and the total maneuver duration decreases; refer to rule of thumb; the total fuel consumption will increase, which is verified through Monte Carlo shooting and shown in Figure5.Figure 5
Eclipse rate influence on mean fuel consumption.Figure6 gives the statistical mean of final relative position under various eclipse rate illumination conditions based on 100-time Monte Carlo shooting. For final relative position control accuracy, due to the maneuver correction, when eclipse rate is small (eclipse rate < 0.3), there is no obvious difference; however if the eclipse is big, the accuracy will be worse because of poor navigation accuracy. From another point of view, 30 percent eclipse in one orbital period is acceptable for visible camera based angles-only navigation.Figure 6
Eclipse rate influence on final relative position control error.


## 5. Conclusion

In this study, the relative rendezvous trajectory is optimized for angles-only navigation between chaser and target, especially for the upcoming small satellites based uncooperative rendezvous mission, such as active debris removal, on-orbit servicing, and inspection. The main purpose of this paper is to design an online rendezvous maneuver planning method which provides better observability and robustness for the whole mission period. The second purpose is to improve fuel consumption efficiency. Considering that there is a tradeoff between the relative position robustness, relative velocity robustness, and fuel consumption, it is important to minimize fuel consumption while ensuring sufficient observability since satellites usually have a tight mass budget.In the context of Yamanaka-Ankersen orbital relative motion equations, the square root unscented Kalman filter (SRUKF) AON algorithm is developed to compute the relative state estimations from a low-volume/mass, power saving, and low-cost optical/infrared camera’s observations. Multi-impulsive Hill guidance law is employed in closed-loop linear covariance analysis model, based on which the quantitative relative position robustness and relative velocity robustness index are defined. By balancing fuel consumption, relative position robustness, and relative velocity robustness, we developed a time and covariance threshold triggered two-level optimal maneuver planning method. Monte Carlo shooting showed that it is feasible to control the TLE level of an uncertain3σr0 = 34.642% of range initial relative state to a 100 m v-bar relative station keeping point reduced to 3σrf = 3.477% of range, under a 30% data gap per revolution on account of the eclipse.More efficient optimization algorithm and SRUKF sigma point sampling strategy will be explored to develop onboard guidance, navigation, and control algorithm for small satellites autonomous rendezvous mission.



---

*Source: 5451908-2017-01-24.xml*