Next Article in Journal
Model and Optimization of the Tether for a Segmented Space Elevator
Next Article in Special Issue
Multi-Layer Fault-Tolerant Robust Filter for Integrated Navigation in Launch Inertial Coordinate System
Previous Article in Journal
Real-Time Fuel Optimization and Guidance for Spacecraft Rendezvous and Docking
Previous Article in Special Issue
Maneuvering Spacecraft Orbit Determination Using Polynomial Representation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Nonlinear Covariance Analysis-Based Robust Rendezvous Trajectory Design by Improved Differential Evolution Method

1
The 54th Research Institute of CETC, Shijiazhuang 050081, China
2
Key Laboratory of Hebei Province on Unmanned System Intelligent Telemetry & Telecontrol Information Technology, Shijiazhuang 050081, China
3
School of Automation, Central South University, Changsha 410083, China
4
College of Intelligence Science and Technology, National University of Defense Technology, Changsha 410073, China
*
Author to whom correspondence should be addressed.
Submission received: 8 March 2022 / Revised: 11 May 2022 / Accepted: 18 May 2022 / Published: 21 May 2022
(This article belongs to the Special Issue Recent Advances in Spacecraft Dynamics and Control)

Abstract

:
This paper presents a robust trajectory design method for approaching and rendezvous with a space target considering multi-source uncertainties. A nonlinear covariance analysis method based on the state transition tensor is presented to formulate the propagation of uncertainties including environment parameter uncertainty, actuator error, sensor noise, navigation error and initial state dispersion of the closed-loop GN&C system. Then, the robust trajectory design problem is defined based on the quantified effect of the uncertainties, and an improved self-adaptive differential evolution algorithm is presented to solve the robust trajectory design problem with uncertainties. Finally, four groups of numerical simulations are carried out to show that the designed robust trajectories can satisfy the final state dispersion constraint under multi-source uncertainties.

1. Introduction

Rendezvous and proximity operation (RPO) with a space target is an essential precondition for the subsequent capture and on-orbit service [1]. The performance of RPO (rendezvous time, energy consumption, etc.) relies heavily on the RPO trajectory, and as a result, trajectory design has been an appealing area for decades. Trajectory design can provide impulse maneuvers needed by chaser spacecraft to approach target spacecraft and achieve goals such as flying by, capture or docking with the target. Lu et al. presented a new methodology for a spacecraft to autonomously rendezvous with a target spacecraft in an arbitrary orbit [2]. The N-impulse rendezvous problem was formalized by Zhou et al. and an improved particle swarm optimization method was proposed to solve the problem [3]. The minimum-fuel trajectory design problem was constructed by Zhao et al. and solved by the indirect optimization method [4]. A novel hierarchical optimization algorithm was proposed in Ref. [5] to speed up the optimization process of spacecraft trajectory design. An artificial emotion memory optimization method has been proposed in Ref. [6] for the trajectory design problem to avoid the local minimum solution. The trajectory design problem was addressed for Hayabusa2 rendezvousing with Ryugu [7] by minimizing the total fuel consumption. Considering the collision problem, Liu et al. proposed a trajectory design method based on the convex optimization algorithm [8]. Taking energy as the performance index, Sandberg et al. presented an autonomous trajectory design algorithm based on Pontryagin’s method [9]. Yang and Li proposed a fuel-optimal trajectory design method based on a novel irregular gravitational Lambert solver [10]. It is noteworthy that most of the RPO trajectory design methods do not take the uncertainties into consideration, but the designed trajectory will be applied in the scenario and spacecraft with uncertainties. The performance of the trajectory may be significantly degraded when considering the multi-source uncertainties.
There are multi-source uncertainties including errors in navigation, guidance and control during the RPO mission, and the uncertainties affect the actual trajectories significantly. Extra energies are needed when tracking the designed trajectories based on the ideal orbit dynamic models. Though the trajectory can be accurately tracked by the advanced tracking control methods [11,12,13,14,15], the whole performance of the GN&C system is not optimal. To solve this issue, researchers have carried out many works on the propagation of uncertainties. A comprehensive overview of the linear and nonlinear uncertainty propagation methods and related applications was provided in Ref. [16], including the Monte Carlo simulation, local linearization, polynomial chaos expansion, state transition tensors (STT), differential algebra, etc. Jones et al. employed the polynomial chaos expansion to give a solution to an uncertain differential equation considering orbit state uncertainty propagation [17]. Valli et al. used differential algebra for the nonlinear propagation of uncertainty, which can realize high-precision estimation with low computational complexity [18]. Uncertain maneuvers of the spacecraft were estimated by Zhai et al. in Ref. [19]. Yang et al. presented a modified uncertainty propagation method on the basis of the STT method [20]. Within the uncertainty propagation methods, the STT-based method is free with random samplings of the multi-source uncertainties, and thus provides an efficient scheme for uncertainty propagation in complex nonlinear dynamic systems.
Based on the uncertainty propagation methods, the uncertainty effect has been considered in the trajectory optimization. Li et al. investigated a multi-objective rendezvous optimization problem under uncertainty and defined a novel objective function containing the final state error [21]. Luo et al. proposed a robust rendezvous trajectory design method based on a novel non-dominated sorting genetic algorithm [22]. Louembet et al. developed a robust guidance algorithm considering the navigation, guidance and control uncertainties [23]. For the entry trajectory design problem, a novel robust trajectory design method was proposed in Ref. [24] by the quantification of both the dynamics of parametric uncertainty and initial state uncertainty. The Mars entry trajectory of Tianwen-1 was designed based on the estimation of uncertain aerodynamic states [25]. The uncertainties of the RPO process have been analyzed and quantized in the open-loop system, and the effect of the uncertainties was evaluated. However, the uncertainties of the closed-loop rendezvous system are rarely analyzed.
Motivated by the foregoing discussions, this work proposes a novel robust rendezvous trajectory design method for the closed-loop RPO system based on the STT method and improved differential evolution algorithm. The contributions of this work are summarized as follows:
  • Uncertainties, including environment parameter uncertainty, actuator error, sensor noise, navigation error and initial state dispersion, of the RPO GN&C closed-loop system are modeled and the propagation, update and correction equations of the uncertainties are formalized based on the STT method.
  • The effect of the uncertainties is rapidly analyzed and predicted quantitatively, and a robust trajectory design problem is constructed considering the uncertainties’ effect.
  • An improved differential evolution algorithm is proposed to solve the robust trajectory design problem for the space rendezvous mission.
The remainder of this paper is organized as follows. Section 2 provides the rendezvous problem models including the dynamics model, environment model, actuator model, sensor model and navigation algorithm model. In Section 3, the effect of the multi-source uncertainties is quantified based on the STT method and the robust performance of the closed-loop GN&C system is evaluated. In Section 4, the robust trajectory design problem is formulated and finally solved by an improved differential evolution algorithm. The simulation results of the robust trajectory design method compared to the traditional optimal trajectory design method and nonlinear Monte Carlo method are presented in Section 5. Conclusions are drawn in Section 6.

2. Rendezvous Dynamic Modeling under Multi-Source Uncertainties

In this section, the rendezvous dynamic models between the service spacecraft (denoted as server for brevity) and the space target are constructed. The constructed models include the orbit dynamic model, environment model, sensor model, actuator model and the navigation algorithm model.

2.1. Orbit Dynamic Model

In the reference inertial frame, the position and velocity vectors of the server are denoted as r c i R 3 and v c i R 3 , and the position and velocity vectors of the considered space target are denoted as r t i R 3 and v t i R 3 , wherein the superscript i indicates the inertial frame, and the subscripts c and t indicate the parameters are associated with the server and the space target, respectively. Then, the orbit dynamic models of the server and target are given as
r ˙ c i = v c i v ˙ c i = μ r c i 3 r c i + a c r c i , v c i , p e n v + u i p Δ v
r ˙ t i = v t i v ˙ t i = μ r t i 3 r t i + a t r t i , v t i , p e n v
where μ is the gravitational coefficient, u denotes the thrust acceleration, and a R 3 is the acceleration due to Earth’s oblateness. p e n v and p Δ v are the uncertainty parameters relevant to the environment model and actuator model, respectively, and will be defined later.

2.2. Environment Model

In the rendezvous process, the uncertainty of the Earth’s oblateness is considered in the environment model and defined as p e n v = η J 2 T , b J 2 T T R 6 , where η J 2 R 3 and b J 2 R 3 are the scale-factor uncertainty and bias uncertainty, respectively. Then, the acceleration a caused by the Earth’s oblateness is modeled as
a = I 3 × 3 + diag η J 2 a ^ + b J 2
where a ^ R 3 is the theoretical value of the acceleration caused by the Earth’s oblateness.

2.3. Actuator Model

In the rendezvous process, all maneuvers are assumed to be impulsive, and the actuator uncertainty parameters are defined as p Δ v = η Δ v T , ϵ Δ v T , b Δ v T T R 9 , where η Δ v R 3 , ϵ Δ v R 3 and b Δ v R 3 are the scale-factor uncertainty, misalignment uncertainty and bias uncertainty, respectively. For the j-th impulsive maneuver Δ v j , the actuator model is constructed as
Δ v j = δ T ϵ Δ v I 3 × 3 + diag η Δ v Δ v ^ j + b Δ v + w Δ v
where Δ v ^ j R 3 is the calculated control instruction given by the guidance algorithm; w Δ v R 3 is a zero mean Gaussian white noise with a constant variance and satisfies
E w Δ v ( t j 1 ) w Δ v T ( t j 2 ) = Q w j δ t j 1 t j 2
where t j 1 and t j 2 are the maneuver time of j 1 -th and j 2 -th maneuvers; Q w j δ t j 1 t j 2 is the standard Kronecker delta function.

2.4. Sensor Model

During the rendezvous process, the state vectors x c = r c i T , v c i T T and x t = r t i T , v t i T T cannot be accurately obtained. Sensors such as the ranging radar and GPS sensors are requisite and the uncertainties of the sensors should be considered and modeled. The sensor uncertainty parameters are defined as p s e n = η s e n T , ϵ s e n T , b s e n T T R 18 , where η s e n R 6 , ϵ s e n R 6 and b s e n R 6 are the scale-factor uncertainty, misalignment uncertainty and bias uncertainty, respectively. For a state vector x ( t ) at time t k , the measurement model of the sensors is constructed as
x ˜ k = δ T ϵ s e n I 6 × 6 + diag η s e n x ( t k ) + b s e n + w s e n
where x ˜ k R 6 is the discrete measurement value of x ( t ) at time t k , and w s e n R 6 is the zero-mean measurement noise, which satisfies
E w s e n ( t k ) w s e n T ( t k ) = Q s e n δ t k t k

2.5. True Dynamic Model

With the constructed models in Section 2.1, Section 2.2, Section 2.3 and Section 2.4, the uncertainty parameters p e n v , p Δ v and p s e n are integrated as p = p e n v T , p Δ v T , p s e n T T = p 1 , p 2 , , p n p T R n p . Without loss of generality, uncertainty parameter p l ( l = 1 , , n p ) is modeled as the first-order Markov process given as
p ˙ l = p l τ l + ω p l , l = 1 , , n p
where τ l is the known time constant and ω p l satisfies
E [ ω p l ( t ) ω p l ( t ) ] = σ p l 2 δ ( t t )
To simplify the presentation, the general expression for the true dynamic model can be reorganized as
x ˙ = f x , u , t
where x : = x c T , x t T , p T T R n p + 12 is the state vector of the true dynamic model.

2.6. Navigation Algorithm Model

After acquiring the measurement values of the system states by the sensors, a navigation algorithm is always needed. Without loss of generality, the most frequently used Kalman algorithm is employed to obtain the navigation states. The navigation state x ^ R n ^ ( n ^ = n p + 12 ) of the true state x is defined as
x ^ = x ^ c T , x ^ t T , p ^ T T
which consists of six chaser states, six target states, and n p parameter states p ^ and given as
p ^ = p ^ e n v T , p ^ l i d a r T , p ^ Δ v T T
The propagation equations of the navigation state are presented as
r ^ ˙ c i = v ^ c i v ^ ˙ c i = μ r ^ c i 3 r ^ c i + a c r ^ c i , v ^ c i , p ^ e n v + u i p ^ Δ v
r ^ ˙ t i = v ^ t i v ^ ˙ t i = μ r ^ t i 3 r ^ t i + a t r ^ t i , v ^ t i , p ^ e n v
p ^ ˙ l = p ^ l τ l , l = 1 , 2 , 3 , , n p
To simplify the presentation, propagation equations of the navigation state x ^ in Equations (13)–(15) are integrated as
x ^ ˙ = f ^ x ^ , u , t
Based on Equation (16), the propagation equation of the navigation state covariance can be given as
P ^ ˙ = F ^ x ^ P ^ + P ^ T F ^ x ^ T + S ^ ω
where P ^ is the navigation state of the covariance matrix, F ^ x ^ is the state-transition matrix of system (16), and S ^ ω is the uncertainty parameter vector.
After acquiring new sensor states, the navigation state and covariance matrix are updated based on
x ^ k + = x ^ k + K ^ k [ z ˜ k z ˜ ^ k ]
P ^ k + = ( I n ^ × n ^ K ^ k H ^ k ) P ^ k ( I n ^ × n ^ K ^ k H ^ k ) T + K ^ k R ^ ν K ^ k T
wherein the symbols + and − represent the navigation state and covariance matrix after and before this update, z ˜ ^ k is the measurement value of the sensors, H ^ k is the measurement sensitivity matrix, and K ^ k is the Kalman filter gain and is presented as
K ^ k = P ^ k H ^ k T ( H ^ k P ^ k H ^ k T + R ^ ν ) 1
where R ^ v is the measurement noise covariance.
When the j-th maneuver occurs, the navigation state and covariance matrix will be corrected. The navigation state x ^ j is corrected according to
x ^ j + c = x ^ j c + d ^ ( x ^ j c , Δ v ^ j )
where the symbols + c and c represent the state after and before the correction, and d ^ ( x ^ j c , Δ v ^ j ) is given as
d ^ ( x ^ j c , Δ v ^ j ) = δ T ^ ϵ Δ v I 3 × 3 + diag η ^ Δ v Δ v ^ j + b ^ Δ v
The covariance matrix is corrected according to
P ^ j + c = [ I + D ^ x ^ ] P ^ j c [ I + D ^ x ^ ] T + S ^ ω j
where D ^ x ^ is the Jacobian matrix of d ^ ( x ^ j c , Δ v ^ j ) with respect to x ^ .

3. Uncertainty Effect Prediction and Robustness Evaluation of GN&C System

On the basis of the rendezvous orbit dynamic models formalized in Section 2, the prediction of the uncertainties’ effect is quantified based on the STT method. Then, the robust performance of the GN&C system against the uncertainties is evaluated quantitatively.

3.1. Uncertainty Effect Prediction Based on STT Method

The models for propagating the true state x as well as the navigation state x ^ can be reorganized based on Equations (10) and (16) as
x ( t ) = Φ ( t ; x j + c , t j )
x ^ ( t ) = Φ ^ ( t ; x ^ j + c , t j )
where Φ and Φ ^ are functions satisfying the conditions
d Φ d t = f ( t , Φ ( t ; x j + c , t j ) )
d Φ ^ d t = f ( t , Φ ^ ( t ; x ^ j + c , t j ) )
In order to analyze the true state x and navigation state x ^ simultaneously, an augmented state vector X is defined as X = x T , x ^ T T , and then the effect of the uncertainties is defined as
δ X = δ x δ x ^ = x x ¯ x ^ x ^ ¯
where x ¯ = E x , x ^ ¯ = E x ^ . Apparently, δ X contains the dispersions of the true state and navigation state, which are caused by the uncertainties. Therefore, the property of δ X can be used to analyze the effect of uncertainties.
This work employs the STT method [20] to obtain an analytical and nonlinear approximation to the solution of δ X . The states’ deviation along the reference states x ¯ and x ^ ¯ are presented based on Einstein’s summation notation as
δ x j ( t ) = p = 1 m 1 p ! Φ j , k 1 k p ( t ) δ x k 1 ( t j + c ) δ x k p ( t j + c )
δ x ^ j ( t ) = p = 1 m 1 p ! Φ ^ j , k 1 k p ( t ) δ x ^ k 1 ( t j + c ) δ x ^ k p ( t j + c )
where the subscript k is the component of each tensor, m is the order of the expansion, and the Φ is the STT of the order p, which can be calculated by
Φ ˙ j , k 1 k p = h A j , k 1 , , A j , k 1 k p ; Φ j , k 1 , , Φ j , k 1 k p
where
A j , k 1 k p = p f j x k 1 x k p
and − indicates that A i , k 1 k p is evaluated over the reference trajectory.
Then, the propagation, update and correction equations of δ X can be written in compact form as
δ X j ( t ) = p = 1 m 1 p ! Φ j , k 1 k p a ( t ) δ X k 1 ( t j + c ) δ X k p ( t j + c )
δ X k + = A k δ X k + B k υ k
δ X j + c = D j δ X j c + N j w j
where
Φ j , k 1 k p a ( t ) = Φ j , k 1 k p ( t ) 0 6 × 6 0 6 × 6 Φ ^ j , k 1 k p ( t )
A k = I 6 × 6 0 6 × 6 K k H k , I 6 × 6 K k H k , B k = 0 6 × 6 K k
D j = I 6 × 6 , BM 0 6 × 6 , I 6 × 6 + BM , N j = B 0 6 × 3
with B = 0 3 × 3 , I 3 × 3 T .
Define the covariance C A of δ X as C A = E [ δ X δ X T ] . As E δ X T = 0 ; then, the propagation, update and correction equations of C A are obtained as
C A i j ( t ) = p = 1 m q = 1 m 1 p ! q ! Φ i , k 1 k p a ( t ) Φ j , s 1 s q a ( t ) · E [ δ X k 1 ( t j + c ) δ X k p ( t j + c ) δ X s 1 ( t j + c ) δ X s q ( t j + c ) ] m i ( t ) m j ( t )
C A ( t k + ) = A k C A ( t k ) A k T + B k R υ B k T
C A ( t j + c ) = D j C A ( t j c ) D j T + N j Q w j N j T
where the mean m i ( t ) can be given as
m i ( t ) = p = 1 m 1 p ! Φ i , k 1 k p a ( t ) E [ δ X k 1 ( t j + c ) δ X k p ( t j + c ) ]

3.2. Robustness Evaluation

The robustness of the GN&C system against multi-source uncertainties can be evaluated by analyzing the covariance matrix of the final true state dispersion. The covariance matrix is obtained by the augmented state covariance matrix as
D t r u e = E δ x t δ x T t = I 6 × 6 0 6 × 6 C A I 6 × 6 0 6 × 6 T
The covariance of the final true state dispersion is caused by the multi-source uncertainties such as the actuator errors and navigation errors. The covariance matrix of the navigation state dispersion δ x ^ ( t ) can also be obtained as
D n a v = E [ δ x ^ ( t ) δ x ^ T ( t ) ] = 0 6 × 6 I 6 × 6 C A 0 6 × 6 I 6 × 6 T
As the desired control values are calculated by the navigation states, which contain dispersions, extra control is needed to correct the errors and this can be evaluated by the control covariance matrix D c o m j as
D c o m j = E δ u j δ u j T = G ^ j D n a v G ^ j T

4. Robust Rendezvous Trajectory Design

In order to take the uncertainties’ effect into consideration in the RPO problem, a robust rendezvous trajectory design method is presented in this section. The objective of the robust rendezvous trajectory design is to propose the optimal maneuver plan that not only considers the total maneuver consumption (quantified by the sum of all the expected maneuvers Δ v ( t j ) 2 ), but also the trajectory robustness against the multi-source uncertainties (quantified by the 3- σ Δ v dispersion and final 3- σ position dispersion).

4.1. Problem Formulation

The robust trajectory design problem is stated as follows:
min j = 1 n ( Δ v ( t j ) 2 + 3 σ Δ v ( t j ) )
s . t . x ˙ = f x , Δ v , t
x t 0 = x 0
x t f = x d
σ r t f σ r r e q
where
σ Δ v ( t j ) = t r a c e D c o m j
σ r t f = t r a c e M r D t r u e t f M r
and D t r u e and D c o m j are given in Equations (43) and (45), σ v ( t j ) is the dispersion of the j-th maneuver, σ r t f is the position dispersion at the final time, M r is the mapping matrices from the true state to the position states, and σ r r e q is the requirement of final state dispersion.
Remark 1.
In the traditional RPO studies, the rendezvous trajectory is usually optimized without considering the effect of uncertainties. In this work, the effect of multi-source uncertainties is considered by σ Δ v ( t j ) and σ Δ v ( t j ) , which are both calculated based on the STT method. If the trajectory optimization problem in Equations (46)–(50) is solved, then the designed trajectory not only has good performance in fuel consumption, but also has strong robustness against the multi-source uncertainties.

4.2. Improved Self-Adaptive Differential Evolution Algorithm

To solve the stochastic optimization problem, an improved self-adaptive differential evolution algorithm for robust trajectory design is presented.
In the differential evolution algorithm, the evolution scaling factor F has an important impact on the optimization searching process. F is always set as a random value between 0 and 1 in the traditional differential evolution algorithm, which leads to the same search range and low efficiency. The robust trajectory design is a complex optimization problem; a self-adaptive search range can improve the search efficiency. To make the algorithm more efficient, a self-adaptive evolution scaling factor is designed and given as
F i , j = F i , j ( Q ( g 1 ) b e s t ) × F i , j ( g ) × r a n d ( 0 , 1 )
where i, j and g represent the sequence numbers of variables, individuals and generations; F i , j ( Q ( g 1 ) b e s t ) is relevant to the best fitness function of the last generation. To ensure that the factors F i , j ( Q ( g 1 ) b e s t ) change from C (C is related with Q ( g 1 ) b e s t ) to 0.1 with the fitness function of the last generation decrease, F i , j ( Q ( g 1 ) b e s t ) is designed as
F i , j ( Q ( g 1 ) b e s t ) = l g ( Q ( g 1 ) b e s t + 10 0.1 ) F i , j ( g ) ( 0.1 , C )
where F i , j ( g ) is an adaptive function related to the generations. To guarantee that the factors F i , j ( g ) vary from 1 to 0.1 with the number of generations increasing, F i , j ( g ) is designed in this work as
F i , j ( g ) = ( cos ( g K × π ) + 11 9 ) × 9 20 F i , j ( g ) ( 0.1 , 1 )
The proposed self-adaptive evolution scaling factor in Equation (55) can ensure that the search range is the whole solution set at the beginning and converges to a smaller set around the best solution, which can obviously improve the search efficiency.
Moreover, the search process of the traditional differential evolution algorithm is sometimes trapped in a local optimal solution. To solve this problem, a random mutant U i , n ( g ) is introduced to enhance its ability to jump out of the local optimal solution. When the search process is trapped in a local optimal solution, the mutant will try to find an individual with better performance, and then restart a new evolution process around this mutant. The random mutant is given as  
U i , n ( g ) = B i L + r a n d j ( 0 , 1 ) × ( B i U B i L ) , g = 1 , , K
where B i L and B i U are lower and upper bounds of the i-th variables.
The proposed improved self-adaptive differential evolution algorithm for robust trajectory design is summarized in Algorithm 1.    
Algorithm 1: Improved self-adaptive differential evolution algorithm for robust trajectory design
Aerospace 09 00277 i001
Remark 2.
It can be obtained from Step 4 in Algorithm 1 that the self-adaptive differential evolution algorithm is not actually applied when the final 3-σ position dispersion constraint in Equation (50) is satisfied under the uncertainties. This means that even applying Δ v ( t j ) of the deterministic trajectory calculated without uncertainties, the final states of the system with uncertainties satisfy the final 3-σ position dispersion constraint. There is no need to start the optimization algorithm as the current maneuver plan is sufficient to meet the needs of the RPO mission.

5. Simulation Results

Numerical simulations are carried out to verify the effectiveness of the novel robust rendezvous trajectory design method based on the STT method and improved differential evolution method.

5.1. Simulation Conditions

The proposed robust rendezvous trajectory design method is verified by a concrete orbit transfer problem. The measured initial relative states are given as x 0 = 0 m , 10 , 000 m , 200 m , 0 m / s , 0 m / s , 0 m / s T , and the expected final relative states are set as x f = 0 m , 200 m , 0 m , 0 m / s , 0 m / s , 0 m / s T . Different desired transfer times are set as t f = 2000 s, 6000 s, 12,000 s and free-time to produce different simulation examples. The deterministic optimal solutions to this scenario without uncertainties for t f = 2000 s, 6000 s and 12,000 s, respectively, are optimized and listed in Table 1.
Different uncertainties are considered in this section. In order to clearly present the uncertainties’ effect on different levels, three typical GN&C systems are defined in Table 2 and denoted as a high-cost system, nominal-cost system and low-cost system, respectively. Here, “high/nominal/low cost” means that the system is equipped with expensive/nominal/cheap sensors and actuators, so the uncertainty levels are different. Different levels of uncertainties (3- σ ) are also given in Table 2, including the initial state dispersions, measurement errors, actuator execution errors and final dispersion requirements. The initial position dispersion ranges from 10 m to 1 km, and the initial velocity dispersion ranges from 0.01 m/s to 1 m/s. The measurement errors range from 0.1 m to 10 m. The actuator execution errors range from 1 mm/s to 100 mm/s. The final position dispersion requirement ranges from 3 m to 30 m.

5.2. Robust Trajectory Design Results

In the proposed self-adaptive evolution algorithm, the maximum number of iterations is set as 100 × number of variables, population size is 1000, and the stopping criterion is that the change in the fitness function is smaller than 10 6 .
Figure 1 shows the optimal fitness function with respect to the number of maneuvers when t f = 2000 s. Though only two maneuvers are needed in the optimal deterministic trajectory, it can be seen from the figure that three maneuvers are needed in the robust optimal trajectory for the high-cost system and four maneuvers are needed in the robust optimal trajectories for the nominal-cost and low-cost systems.
Table 3 presents the simulation results when t f = 2000 s. For the total Δ v representing the fuel consumption, it is obtained that more Δ v is required to satisfy the constraint of final position dispersion when there are stronger uncertainties. This is because extra actuator execution errors are injected into the system due to extra maneuvers. For the optimal maneuver times, it is seen that the high-cost system requires only one mid-course correction at 576 s, the nominal-cost system requires two mid-course corrections—one early at 423 s and one late at 1504 s—and the low-cost system requires two very late mid-course corrections at 1406 s and 1880 s. The late maneuvers for the low-cost system are not surprising since injecting large maneuver execution errors earlier in the trajectory creates a larger final position dispersion.
Figure 2 and Figure 3 show the in-plane and out-of-plane optimal trajectories (notation cvx indicates the optimal trajectory solved by convex optimization without uncertainties) and the associated maneuver locations (denoted as M1, M2, etc.) when t f = 2000 s. It can be seen from the two figures that the designed robust optimal trajectories are basically the same as the deterministic optimal trajectory, which means that the main objective of the robust trajectory optimization algorithm is to correct the trajectory against the uncertainties rather than replan a new trajectory.
Figure 4 presents the Monte Carlo analysis for the low-cost system when t f = 2000 s. A total of 1000 independent Monte Carlo simulations were carried out with the stochastic uncertainties and are shown in the figure with blue and solid lines. The uncertainty effect on the 3- σ position dispersion is also analyzed by the proposed STT-based method and the result is shown in the figure with red and dashed line. Meanwhile, the yellow and dotted line in the figure denotes the mission constraint of the final position dispersion. It is verified that the proposed STT-based uncertainty effect evaluation method can provide the accurate upper bound of the Monte Carlo results. Moreover, all Monte Carlo trajectories satisfy the final position dispersion constraint, which indicates that the proposed robust trajectory design method is effective when considering strong uncertainties.
Figure 5 shows the optimal fitness function with respect to the number of maneuvers when t f = 6000 s. Four maneuvers are needed in the robust optimal trajectories for both the high- and nominal-cost systems, while five maneuvers are needed in the robust optimal trajectories for the low-cost system.
Table 4 presents the simulation results when t f = 6000 s. For the optimal maneuver times, it is seen that the high-cost system requires two mid-course corrections at 1127 s and 4735 s, the nominal cost system requires two mid-course corrections—one at 3868 s and one late at 5162 s—and the low-cost system requires three mid-course corrections at 1377 s, 4741 s and 5829 s. It is noteworthy that the designed maneuvers are completely the same in the deterministic trajectory and the high-cost system trajectory. This is because even applying Δ v ( t j ) of the deterministic trajectory into the high-cost system, the final 3- σ position dispersion constraint is still satisfied and, as a result, the self-adaptive differential evolution algorithm is not run to avoid unnecessary calculations. A detailed explanation can be found in Algorithm 1 and Remark 2.
Figure 6 and Figure 7 show the in-plane and out-of-plane optimal trajectories when t f = 6000 s, as well as the maneuver locations. It can be seen from the figures that the designed robust optimal trajectory of the high-cost system is basically the same as the deterministic optimal trajectory due to the same maneuvers, while the designed robust optimal trajectories of the nominal and low-cost systems are different from the deterministic optimal trajectory. A conclusion is drawn that the proposed method is capable of replanning a new trajectory when the uncertainty effect is great to the rendezvous system.
Figure 8 presents the Monte Carlo analysis for the high-cost system when t f = 6000 s. It is verified that the proposed STT-based uncertainty effect evaluation method can provide the accurate upper bound of the Monte Carlo results. Moreover, all Monte Carlo trajectories satisfy the final position dispersion constraint.
Figure 9 shows the optimal fitness function with respect to the number of maneuvers when t f = 12,000 s. Four maneuvers are needed in the robust optimal trajectories for both the high- and nominal-cost systems, while five maneuvers are needed in the robust optimal trajectories for the low-cost system.
Table 5 presents the simulation results when t f = 12,000 s. For the optimal maneuver times, it is interesting to note that the time of the first maneuver is not 0 in all stochastic cases. This is because there are initial position and velocity dispersions at t = 0 s. For a long transfer time, rather than using a poor measured state with large dispersion to compute the maneuver, it would be better to wait until the navigation algorithm has run for some time and the measured states are more accurate.
Figure 10 and Figure 11 show the in-plane and out-of-plane optimal trajectories when t f = 12,000 s, as well as the optimal maneuver locations. In this case, the optimal maneuvers change the desired optimal trajectory.
Figure 12 presents the Monte Carlo analysis for the high-cost system when t f = 12,000 s. It is verified that the proposed STT-based uncertainty effect evaluation method can provide the accurate upper bound of the Monte Carlo results. Moreover, all Monte Carlo trajectories satisfy the final position dispersion constraint.
Table 6 presents the simulation results when the final time is free. Four maneuvers are needed in the robust optimal trajectories for both the high- and nominal-cost systems, while five maneuvers are needed in the robust optimal trajectories for the low-cost system. As expected, the optimal fitness functions in all systems are smaller than the fixed final time case t f = 12,000 s. Theoretically, when the system is free with uncertainties, Δ v should approach zero when the final time is free, and the final time should approach infinity. When uncertainties are considered, however, an excessively long final time will lead to a large final state dispersion. It can be observed that the optimized final time is shorter when the uncertainties are stronger, which provides clear evidence of the uncertainty effect on the rendezvous system.
Figure 13 and Figure 14 show the in-plane and out-of-plane trajectories when t f = free, as well as the optimal maneuver locations.

6. Conclusions

A novel robust trajectory design method has been proposed in this work for the rendezvous and proximity operations mission in Earth orbit. Optimal maneuver plans have been provided that not only satisfy the final state dispersion constraints but minimize the effect of multi-source uncertainties. The main contributions and advantages are twofold: on the one hand, multi-source uncertainties are modeled of the closed-loop rendezvous system and the effect of the uncertainties is quantized and evaluated based on the state transition tensor method. On the other hand, the effect of the uncertainties is considered in the robust rendezvous trajectory design problem, which is effectively solved by an improved differential evolution algorithm. Numerical simulations clearly show that the multi-source uncertainties, including environment parameter uncertainty, actuator error, sensor noise, navigation error and initial state dispersion, may have significant effects on the designed trajectory. By employing the proposed method, the effect of uncertainties is effectively considered and the optimized robust rendezvous trajectory satisfies the mission requirements and has strong robustness against uncertainties. In future work, complex constraints including a no-fly zone constraint, collision avoidance constraint and docking axis constraint will be considered in the proposed robust trajectory design framework.

Author Contributions

Conceptualization, K.J. and Z.Y.; methodology, K.J.; software, K.J.; validation, K.J., Z.Y. and Y.L.; resources, Y.L. and X.Z.; writing—original draft preparation, K.J.; writing—review and editing, Z.Y.; visualization, Y.Z.; supervision, Y.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (62103446, 11902346).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Li, W.J.; Cheng, D.Y.; Liu, X.G.; Wang, Y.B.; Shi, W.H.; Tang, Z.X.; Gao, F.; Zeng, F.M.; Chai, H.Y.; Luo, W.B.; et al. On-orbit service (OOS) of spacecraft: A review of engineering developments. Prog. Aerosp. Sci. 2019, 108, 32–120. [Google Scholar] [CrossRef]
  2. Lu, P.; Liu, X. Autonomous trajectory planning for rendezvous and proximity operations by conic optimization. J. Guid. Control. Dyn. 2013, 36, 375–389. [Google Scholar] [CrossRef]
  3. Zhou, H.; Wang, X.; Cui, N. Fuel-optimal multi-impulse orbit transfer using a hybrid optimization method. IEEE Trans. Intell. Transp. Syst. 2019, 21, 1359–1368. [Google Scholar] [CrossRef]
  4. Zhao, S.; Zhang, J.; Xiang, K.; Qi, R. Target sequence optimization for multiple debris rendezvous using low thrust based on characteristics of SSO. Astrodynamics 2017, 1, 85–99. [Google Scholar] [CrossRef]
  5. He, H.; Shi, P.; Zhao, Y. Hierarchical optimization algorithm and applications of spacecraft trajectory optimization. Aerospace 2022, 9, 81. [Google Scholar] [CrossRef]
  6. Shengnan, F.; Liang, W.; Qunli, X. Trajectory optimization of a reentry vehicle based on artificial emotion memory optimization. J. Syst. Eng. Electron. 2021, 32, 668–680. [Google Scholar] [CrossRef]
  7. Tsuda, Y.; Takeuchi, H.; Ogawa, N.; Ono, G.; Kikuchi, S.; Oki, Y.; Ishiguro, M.; Kuroda, D.; Urakawa, S.; Okumura, S.I. Rendezvous to asteroid with highly uncertain ephemeris: Hayabusa2’s Ryugu-approach operation result. Astrodynamics 2020, 4, 137–147. [Google Scholar] [CrossRef]
  8. Liu, X.; Yang, H.; Li, S. Collision-free trajectory design for long-distance hopping transfer on asteroid surface using convex optimization. IEEE Trans. Aerosp. Electron. Syst. 2021, 57, 3071–3083. [Google Scholar] [CrossRef]
  9. Sandberg, A.; Sands, T. Autonomous trajectory generation algorithms for spacecraft slew maneuvers. Aerospace 2022, 9, 135. [Google Scholar] [CrossRef]
  10. Yang, H.; Li, S. Fuel-optimal asteroid descent trajectory planning using a lambert solution-based costate initialization. IEEE Trans. Aerosp. Electron. Syst. 2020, 56, 4338–4352. [Google Scholar] [CrossRef]
  11. Hu, Q.; Chen, W.; Guo, L. Fixed-time maneuver control of spacecraft autonomous rendezvous with a free-tumbling target. IEEE Trans. Aerosp. Electron. Syst. 2018, 55, 562–577. [Google Scholar] [CrossRef]
  12. Mancini, M.; Bloise, N.; Capello, E.; Punta, E. Sliding mode control techniques and artificial potential field for dynamic collision avoidance in rendezvous maneuvers. IEEE Control. Syst. Lett. 2019, 4, 313–318. [Google Scholar] [CrossRef]
  13. Wei, C.; Wu, X.; Xiao, B.; Wu, J.; Zhang, C. Adaptive leader-following performance guaranteed formation control for multiple spacecraft with collision avoidance and connectivity assurance. Aerosp. Sci. Technol. 2022, 120, 107266. [Google Scholar] [CrossRef]
  14. Shao, X.; Hu, Q.; Shi, Y. Adaptive pose control for spacecraft proximity operations with prescribed performance under spatial motion constraints. IEEE Trans. Control. Syst. Technol. 2020, 29, 1405–1419. [Google Scholar] [CrossRef]
  15. Wei, C.; Gui, M.; Zhang, C.; Liao, Y.; Dai, M.Z.; Luo, B. Adaptive appointed-time consensus control of networked Euler-Lagrange systems with connectivity preservation. IEEE Trans. Cybern. 2021, in press. [Google Scholar] [CrossRef]
  16. Luo, Y.z.; Yang, Z. A review of uncertainty propagation in orbital mechanics. Prog. Aerosp. Sci. 2017, 89, 23–39. [Google Scholar] [CrossRef]
  17. Jones, B.A.; Parrish, N.; Doostan, A. Postmaneuver collision probability estimation using sparse polynomial chaos expansions. J. Guid. Control. Dyn. 2015, 38, 1425–1437. [Google Scholar] [CrossRef] [Green Version]
  18. Valli, M.; Armellin, R.; Di Lizia, P.; Lavagna, M. Nonlinear mapping of uncertainties in celestial mechanics. J. Guid. Control. Dyn. 2012, 36, 48–63. [Google Scholar] [CrossRef]
  19. Zhai, G.; Wang, Y.; Zhao, Q. Tracking the maneuvering spacecraft propelled by swing propulsion of constant magnitude. J. Syst. Eng. Electron. 2020, 31, 370–382. [Google Scholar] [CrossRef]
  20. Yang, Z.; Luo, Y.Z.; Zhang, J. Nonlinear semi-analytical uncertainty propagation of trajectory under impulsive maneuvers. Astrodynamics 2019, 3, 61–77. [Google Scholar] [CrossRef]
  21. Li, H.Y.; Luo, Y.Z.; Tang, G.J. Optimal multi-objective linearized impulsive rendezvous under uncertainty. Acta Astronaut. 2010, 66, 439–445. [Google Scholar] [CrossRef]
  22. Luo, Y.; Yang, Z.; Li, H. Robust optimization of nonlinear impulsive rendezvous with uncertainty. Sci. China Phys. Mech. Astron. 2014, 57, 731–740. [Google Scholar] [CrossRef]
  23. Louembet, C.; Arzelier, D.; Deaconu, G. Robust rendezvous planning under maneuver execution errors. J. Guid. Control. Dyn. 2014, 38, 76–93. [Google Scholar] [CrossRef] [Green Version]
  24. Jiang, X.; Li, S. Mars entry trajectory planning using robust optimization and uncertainty quantification. Acta Astronaut. 2019, 161, 249–261. [Google Scholar] [CrossRef]
  25. Wei, H.; Rao, W.; Chen, G.; Wang, G.; Zou, X.; Li, Q.; Hu, Y. Tianwen-1 Mars entry vehicle trajectory and atmosphere reconstruction preliminary analysis. Astrodynamics 2022, 6, 81–91. [Google Scholar] [CrossRef]
Figure 1. Optimal fitness function with respect to the number of maneuvers when t f = 2000 s.
Figure 1. Optimal fitness function with respect to the number of maneuvers when t f = 2000 s.
Aerospace 09 00277 g001
Figure 2. Optimal trajectories (in-plane) when t f = 2000 s.
Figure 2. Optimal trajectories (in-plane) when t f = 2000 s.
Aerospace 09 00277 g002
Figure 3. Optimal trajectories (out-of-plane) when t f = 2000 s.
Figure 3. Optimal trajectories (out-of-plane) when t f = 2000 s.
Aerospace 09 00277 g003
Figure 4. Monte Carlo analysis for low-cost system when t f = 2000 s.
Figure 4. Monte Carlo analysis for low-cost system when t f = 2000 s.
Aerospace 09 00277 g004
Figure 5. Optimal fitness function with respect to the number of maneuvers when t f = 6000 s.
Figure 5. Optimal fitness function with respect to the number of maneuvers when t f = 6000 s.
Aerospace 09 00277 g005
Figure 6. Optimal trajectories (in-plane) when t f = 6000 s.
Figure 6. Optimal trajectories (in-plane) when t f = 6000 s.
Aerospace 09 00277 g006
Figure 7. Optimal trajectories (out-of-plane) when t f = 6000 s.
Figure 7. Optimal trajectories (out-of-plane) when t f = 6000 s.
Aerospace 09 00277 g007
Figure 8. Monte Carlo analysis for high-cost system when t f = 6000 s.
Figure 8. Monte Carlo analysis for high-cost system when t f = 6000 s.
Aerospace 09 00277 g008
Figure 9. Optimal fitness function with respect to the number of maneuvers when t f = 12,000 s.
Figure 9. Optimal fitness function with respect to the number of maneuvers when t f = 12,000 s.
Aerospace 09 00277 g009
Figure 10. Optimal trajectories(in-plane) when t f = 12,000 s.
Figure 10. Optimal trajectories(in-plane) when t f = 12,000 s.
Aerospace 09 00277 g010
Figure 11. Optimal trajectories (out-of-plane) when t f = 12,000 s.
Figure 11. Optimal trajectories (out-of-plane) when t f = 12,000 s.
Aerospace 09 00277 g011
Figure 12. Monte Carlo analysis for high-cost system when t f = 12,000 s.
Figure 12. Monte Carlo analysis for high-cost system when t f = 12,000 s.
Aerospace 09 00277 g012
Figure 13. Optimal trajectories (in-plane) when t f = free.
Figure 13. Optimal trajectories (in-plane) when t f = free.
Aerospace 09 00277 g013
Figure 14. Optimal trajectories (out-of-plane) when t f = free.
Figure 14. Optimal trajectories (out-of-plane) when t f = free.
Aerospace 09 00277 g014
Table 1. Deterministic maneuver data.
Table 1. Deterministic maneuver data.
t f (s)Total Δ v (m/s)Number of ManeuversManeuver Times (s) Δ v (m/s)Maneuver Positions [ x y z ] (km)
20009.4802 0 2000 4.737 4.743 [0 −10 0.2]
[0 0 0]
60001.2674 0 1127 4735 6000 0.414 0.221 0.160 0.472 [0 −10 0.2]
[−0.545 −9.950 0.012]
[−0.744 −0.410 0.075]
[0 0 0]
12,0000.6654 0 1359 10,183 12,000 0.133 0.202 0.071 0.259 [0 −10 0.2]
[−0.227 −9.918 −0.020]
[−0.639 −0.815 0.073]
[0 0 0]
Table 2. Different levels of uncertainties (3- σ ) for high-, nominal- and low-cost systems.
Table 2. Different levels of uncertainties (3- σ ) for high-, nominal- and low-cost systems.
UncertaintyValue (3- σ )Unit
High CostNominal CostLow Cost
Initial State Dispersion
Position101001000m
Velocity0.010.11m/s
Measurement Error
Noise0.1110m
Misalignment0.1110mrad
Scale factor1005001000ppm
Bias0.010.030.05m
Actuator Execution Error
Noise0.0010.010.1m/s
Misalignment0.115mrad
Scale factor1005001000ppm
Bias0.115mm/s
Final Dispersion Requirement
Position31530m
Table 3. Optimal trajectory results when t f = 2000 s.
Table 3. Optimal trajectory results when t f = 2000 s.
Case 1, t f = 2000 s
SystemDeterministicHigh CostNominalLow Cost
Fitness (m/s)NA9.537 10.009 14.181
Δ v (m/s)9.4809.481 9.488 10.053
3- σ Δ v (m/s)NA0.057 0.520 4.128
3- σ r (m)NA2.999 14.968 20.246
Opt. no. of maneuvers2344
Maneuver times (s) 0 2000 0 576 2000 0 423 1504 2000 80 1406 1888 2000
Table 4. Optimal trajectory results when t f = 6000 s.
Table 4. Optimal trajectory results when t f = 6000 s.
Case 2, t f = 6000 s
DeterministicHigh CostNominalLow Cost
Fitness (m/s)NA1.3321.8205.892
Δ v (m/s)1.2671.2671.4161.354
3- σ Δ v (m/s)NA0.0650.4044.538
3- σ r (m)NA2.90214.99929.939
Opt. no. of maneuvers4445
Maneuver times (s) 0 1127 4735 6000 0 1127 4735 6000 180 3868 5162 6000 62 1377 4741 5829 6000
Table 5. Optimal trajectory results when t f = 12,000 s.
Table 5. Optimal trajectory results when t f = 12,000 s.
Case 3, t f = 12,000 s
DeterministicHigh CostNominalLow Cost
Fitness (m/s)NA0.7961.2585.437
Δ v (m/s)0.6650.7270.7850.918
3- σ Δ v (m/s)NA0.0690.4734.519
3- σ r (m)NA2.99414.98629.959
Opt. no. of maneuvers4445
Maneuver times (s) 0 1359 10,183 12,000 858 7343 10,572 12,000 171 7021 11,162 12,000 161 4222 10,527 11,828 12,000
Table 6. Stochastic optimal performance results, t f = free.
Table 6. Stochastic optimal performance results, t f = free.
Case 4, t f = free
SystemHigh CostNominalLow Cost
Final time (s)17,04316,96410,086
Fitness (m/s)0.6511.0735.424
Δ v (m/s)0.6020.6091.018
3- σ Δ v (m/s)0.0490.4644.407
3- σ r (m)2.99414.96729.987
Opt. no. of maneuvers445
Maneuver times (s) 590 11,009 15,615 17,043 179 12,629 16,132 16,964 110 4214 8752 9914 10,086
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jin, K.; Yin, Z.; Lei, Y.; Zhang, Y.; Zhang, X. Nonlinear Covariance Analysis-Based Robust Rendezvous Trajectory Design by Improved Differential Evolution Method. Aerospace 2022, 9, 277. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace9050277

AMA Style

Jin K, Yin Z, Lei Y, Zhang Y, Zhang X. Nonlinear Covariance Analysis-Based Robust Rendezvous Trajectory Design by Improved Differential Evolution Method. Aerospace. 2022; 9(5):277. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace9050277

Chicago/Turabian Style

Jin, Kai, Zeyang Yin, Yaolin Lei, Yuanlong Zhang, and Xiaolong Zhang. 2022. "Nonlinear Covariance Analysis-Based Robust Rendezvous Trajectory Design by Improved Differential Evolution Method" Aerospace 9, no. 5: 277. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace9050277

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop