Optimal planning approaches with multiple impulses for rendezvous based on hybrid genetic algorithm and control method

Abstract

In this article, we focus on safe and effective completion of a rendezvous and docking task by looking at planning approaches and control with fuel-optimal rendezvous for a target spacecraft running on a near-circular reference orbit. A variety of existent practical path constraints are considered, including the constraints of field of view, impulses, and passive safety. A rendezvous approach is calculated by using a hybrid genetic algorithm with those constraints. Furthermore, a control method of trajectory tracking is adopted to overcome the external disturbances. Based on Clohessy–Wiltshire equations, we first construct the mathematical model of optimal planning approaches of multiple impulses with path constraints. Second, we introduce the principle of hybrid genetic algorithm with both stronger global searching ability and local searching ability. We additionally explain the application of this algorithm in the problem of trajectory planning. Then, we give three-impulse simulation examples to acquire an optimal rendezvous trajectory with the path constraints presented in this article. The effectiveness and applicability of the tracking control method are verified with the optimal trajectory above as control objective through the numerical simulation.

On-orbit service, autonomous rendezvous, and docking play the part of key technology in future space activities such as orbital life extension of geosynchronous earth orbit (GEO) satellites, capture of passive satellites, and the space debris removal. ¹ These technologies require close-proximity operation including the inspection phase and the final approach phase. ² Regarding these space missions, the difficulty of research is compounded by problems related to planning approaches and control with fuel-optimal rendezvous in the close-proximity operation phases.

Facing the trajectory design of close-range rendezvous, scholars have made fruitful achievements using primer vector theory proposed by Lawden ³ with the first-order necessary conditions for propellant-optimal trajectories. Prussing and Chiu ⁴ obtained the optimal, time-fixed, multiple-impulse solutions for both coplanar and a restricted class of noncoplanar circle-to-circle rendezvous. For the fixed-time problem, Carter and colleagues^5,6 present new necessary and sufficiency conditions that do not require global calculation of the primer vector. However, path constraints, such as the field of view (FOV) angle, impulse, and trajectory safety, are not considered in those studies. The eventual choice of acquisition strategy for close-range rendezvous will depend on those path constraints. ⁷

Luo et al. ⁸ investigated multi-objective impulsive rendezvous optimization with trajectory safety and FOV constraints. Epenoy ⁹ formalized minimum-fuel rendezvous as an optimal control problem subject to a collision-avoidance constraint and dealt with it by smoothed exact penalty approach. Luo et al. ¹⁰ established an optimization model including several different kinds of impulse constraints for the minimum-time multiple-impulse rendezvous. Treating rendezvous safety as a constraint, Breger and How ¹¹ present a method of on-orbit generation focused on optimal safe trajectory. Although much work has been done on optimal impulsive rendezvous with one or two constraints, researches are seldom concerning on optimal rendezvous trajectory design with all those three path constraints. The goal of this article is to design minimum-fuel trajectory by hybrid genetic algorithm (GA) to cope with all these constraints.

Therefore, based on the materials and surveys, we consider a series of practical-project constraints, such as FOV, maximum impulses’ magnitude, minimum time between imposing impulses, and passive safety. For a practical rendezvous project, the optical sensors used for relative navigation in the close-range rendezvous pose the limited FOV constraint. The thrust magnitude and the engine working time of the thrust engine onboard the chaser spacecraft raise the impulse constraints. Safety constraint in the context of spacecraft rendezvous refers to collision avoidance in case of system failure. ¹²

In this article, we design and present an optimal safety approach using a hybrid GA to minimize the overall consumption of fuel. Furthermore, with the optimal trajectory as our control objective, we adopt the tracking control method and thereby overcome the disturbances and errors existing in the practical condition, satisfying the terminal-states requirement of the rendezvous mission. The simulation results reveal that the hybrid GA which contains the sequential quadratic programming (SQP) can effectively deal with the rendezvous constraints and can improve both local search ability and search accuracy of the traditional GA, thus enhancing its credibility.

Analysis of relative motion

When studying the rendezvous trajectory for target spacecraft running on a near-circular reference orbit, it is common to adopt Clohessy–Wiltshire (C–W) equations in the local vertical, local horizontal (LVLH) frame. This system is defined as the coordinate frame $O-\mathrm{xyz}$ with right-handed Cartesian coordinates, where origin $O$ is attached to the target spacecraft’s center of mass, the y-axis runs along the vector $r$ to Earth’s center of mass, the x-axis follows the target orbit circumference, vertical with y-axis, along the speed direction, and the z-axis completes the right-handed frame, as shown in Figure 1. ¹³

Figure 1.

Target orbit coordinate system.

In accordance with the two-body dynamics theory, one subtracts between the absolute positions of the two spacecrafts. Assuming that Earth is a uniform-quality sphere, spacecrafts are affected by only central gravity, not by the perturbation forces from environment. The distance between the target and the chaser spacecraft is small compared to the distance between the target and the mass center of the Earth for the close-range rendezvous. The linearized equations for relative motion with respect to an unperturbed elliptical reference orbit with a linear differential gravity field can be obtained ¹⁴

$$\{\begin{array}{c}\hfill \stackrel{··}{x}+2\stackrel{·}{\theta }\stackrel{·}{y}-{\stackrel{·}{\theta }}^{2}x+\stackrel{··}{\theta }y+\mu x/r_t^3=a_x\hfill \\ \multicolumn{1}{c}{\stackrel{··}{y}-2\stackrel{·}{\theta }\stackrel{·}{x}-{\stackrel{·}{\theta }}^{2}y-\stackrel{··}{\theta }x-2\mu y/r_t^3=a_y}\\ \multicolumn{1}{c}{\stackrel{··}{z}+\mu z/r_t^3=a_z}\end{array}$$

(1)

where $r_t$ means the relative position between Earth’s mass center and target spacecraft; $\stackrel{·}{\theta }$ stands for orbit angular velocity of the target spacecraft; $\stackrel{··}{\theta }$ indicates orbit angular acceleration of the target spacecraft; $\mu$ is the constant of Earth gravitation; $a_x$ , $a_y$ , and $a_z$ stand for maneuvering acceleration along the direction of three axes in LVLH coordinate system, respectively; and $x$ , $y$ , and $z$ mean the relative positions between the two spacecrafts along the three axes in LVLH coordinate system, respectively.

C–W equations are the specific forms for the target spacecraft running on near-circular orbit, that is, $e\approx 0$ .

So, $\stackrel{·}{\theta }=(\omega {(1+e\cos \theta )}^2)/{(1-e^2)}^{3/2}=\omega$ and $\stackrel{··}{\theta }=(-2{\omega }^{2}e\sin \theta {(1+e\cos \theta )}^3)/{(1-e^2)}^3=0$ ; then substituting these expressions into formula (1), we get the C–W equations ¹⁵

$$\{\begin{array}{c}\hfill \stackrel{··}{x}+2\omega \stackrel{·}{y}=a_x\hfill \\ \multicolumn{1}{c}{\stackrel{··}{y}-2\omega \stackrel{·}{x}-3{\omega }^{2}y=a_y}\\ \multicolumn{1}{c}{\stackrel{··}{z}+{\omega }^{2}z=a_z}\end{array}$$

(2)

where $e$ , $\omega$ , and $\theta$ stand for eccentricity, mean angular velocity, and true anomaly of target spacecraft orbit, respectively.

The linearized equations, as shown in equation (2), are written in six-dimensional matrix forms

$$\stackrel{·}{\mathit{X}}=\mathit{AX}+\mathit{BU}$$

(3)

where $\mathit{X}={[\begin{array}{cccccc}\hfill x\hfill & \hfill y\hfill & \hfill z\hfill & \hfill \stackrel{·}{x}\hfill & \hfill \stackrel{·}{y}\hfill & \hfill \stackrel{·}{z}\hfill \end{array}]}^{\mathrm{T}}$ , $\mathit{U}={[\begin{array}{ccc}\hfill a_x\hfill & \hfill a_y\hfill & \hfill a_z\hfill \end{array}]}^{\mathrm{T}}$ , $\mathit{A}(\theta )=\left[\begin{array}{cc}\hfill {\mathit{O}}_{3\times 3}\hfill & \hfill {\mathit{I}}_{3\times 3}\hfill \\ \multicolumn{1}{c}{{\mathit{F}}_1}& \hfill {\mathit{F}}_2\hfill \end{array}\right]$ , $\mathit{B}=\left[\begin{array}{c}\hfill {\mathit{O}}_{3\times 3}\hfill \\ \multicolumn{1}{c}{{\mathit{I}}_{3\times 3}}\end{array}\right]$ , ${\mathit{F}}_1=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill \stackrel{··}{\theta }\hfill \\ \multicolumn{1}{c}{0}& \hfill 3{\omega }^2\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill -{\omega }^2\hfill \end{array}\right]$ , ${\mathit{F}}_2=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill -2\omega \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{2\omega }& \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill -{\omega }^2\hfill \end{array}\right]$ , and assuming initial time is $t_0$ , equation (3) has the following general homogenous solution for relative motion without control for any time $t$

$$\mathit{X}(t)=\mathbf{\Phi }(\Delta t)\mathit{X}(t_0)$$

(4)

where $\mathbf{\Phi }(\Delta t)=\left[\begin{array}{cccccc}\hfill 1\hfill & \hfill 6(s-\omega \Delta t)\hfill & \hfill 0\hfill & \hfill \frac{4s-3\omega \Delta t}{\omega }\hfill & \hfill \frac{2(c-1)}{\omega }\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 4-3c\hfill & \hfill 0\hfill & \hfill \frac{2(1-c)}{\omega }\hfill & \hfill \frac{s}{\omega }\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill c\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{s}{\omega }\hfill \\ \multicolumn{1}{c}{0}& \hfill 6\omega (c-1)\hfill & \hfill 0\hfill & \hfill 4c-3\hfill & \hfill -2s\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 3\omega s\hfill & \hfill 0\hfill & \hfill 2s\hfill & \hfill c\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill -\omega s\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill c\hfill \end{array}\right]$ , $\Delta t=t-t_0$ , $s=\sin \omega \Delta t$ , and $c=\cos \omega \Delta t$ .

The state transition matrix $\mathbf{\Phi }(\Delta t)$ is divided into four parts

$$\mathbf{\Phi }(\Delta t)=\left[\begin{array}{cc}\hfill {\mathbf{\Phi }}_{11}\hfill & \hfill {\mathbf{\Phi }}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{\Phi }}_{21}}& \hfill {\mathbf{\Phi }}_{22}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill {\mathbf{\Phi }}_r\hfill & \hfill {\mathbf{\Phi }}_v\hfill \end{array}\right]$$

(5)

Using the C–W equations, the position and velocity vector of the chaser at time $t$ can be expressed

$$\{\begin{array}{c}\hfill \mathit{r}(t)={\mathbf{\Phi }}_{\mathrm{rr}}(\Delta t){\mathit{r}}_0+{\mathbf{\Phi }}_{\mathrm{rv}}(\Delta t){\mathit{v}}_0\hfill \\ \multicolumn{1}{c}{\mathit{v}(t)={\mathbf{\Phi }}_{\mathrm{vr}}(\Delta t){\mathit{r}}_0+{\mathbf{\Phi }}_{\mathrm{vv}}(\Delta t){\mathit{v}}_0}\end{array}$$

(6)

where $\mathbf{r}={[\begin{array}{ccc}\hfill x\hfill & \hfill y\hfill & \hfill z\hfill \end{array}]}^T$ and $\mathbf{v}={[\begin{array}{ccc}\hfill \stackrel{·}{x}\hfill & \hfill \stackrel{·}{y}\hfill & \hfill \stackrel{·}{z}\hfill \end{array}]}^T$ .

The problem of rendezvous trajectory planning with multiple impulses

Path constraints

FOV angle

For a practical rendezvous project, the optical sensors used for relative navigation in the close-range rendezvous pose an additional operational constraint for the approach trajectory design, that is, a limited FOV. ¹⁶ Figure 2 describes the constraint of field view on the rendezvous trajectory. Abscissa represents the direction of the x-axis in the relative coordinate system, and ordinate represents the direction of the y-axis. The field angle located on the trajectory of the chaser spacecraft is defined as $\arctan (y/x)$ . Since the in-plane and out-of-plane motions are decoupled in equation (4) and the out-of-plane motion can be canceled by a single impulse, only the in-plane motion is considered in the constraints.

Figure 2.

Field of view constraint of measurement facility on rendezvous trajectory.

If the impulse maneuver $\Delta v_i$ is applied at $t_i$ , then the position is $r_i$ and velocity is $v_i^{-}$ . A superscript “−” sign is used to denote a condition preceding an impulsive maneuver. The time between two impulses is defined as $t$ ( $t_i\le t\le t_{i+1}$ ). Based on equation (6) and these conventions, we can acquire the chaser’s position after ith burn as

$$\mathbf{r}(t)={\mathbf{\Phi }}_{11}(t-t_i){\mathbf{r}}_i+{\mathbf{\Phi }}_{12}(t-t_i){\mathbf{v}}_i^{-}+{\mathbf{\Phi }}_{12}(t-t_i)\Delta {\mathbf{v}}_i$$

(7)

Assuming that $\Delta t=t-t_i$ is applied, equations (4) and (5) are substituted in equation (7), so we obtain the following equation

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{c}\hfill x\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{z}\end{array}\right]=\left[\begin{array}{c}\hfill x_i+6(s-\omega \Delta t)y_i\hfill \\ \multicolumn{1}{c}{(4-3c)y_i}\\ \multicolumn{1}{c}{{\mathrm{cz}}_i}\end{array}\right]+\left[\begin{array}{c}\hfill \frac{4s-3\omega \Delta t}{\omega }{v}_{\mathrm{ix}}^{-}+\frac{2(c-1)}{\omega }{v}_{\mathrm{iy}}^{-}\hfill \\ \multicolumn{1}{c}{\frac{2(1-c)}{\omega }{v}_{\mathrm{ix}}^{-}+\frac{s}{\omega }{v}_{\mathrm{iy}}^{-}}\\ \multicolumn{1}{c}{\frac{s}{\omega }{v}_{\mathrm{iz}}^{-}}\end{array}\right]}\\ \multicolumn{1}{c}{+\left[\begin{array}{c}\hfill \frac{4s-3\omega \Delta t}{\omega }\Delta v_{\mathrm{ix}}+\frac{2(c-1)}{\omega }\Delta v_{\mathrm{iy}}\hfill \\ \multicolumn{1}{c}{\frac{2(1-c)}{\omega }\Delta v_{\mathrm{ix}}+\frac{s}{\omega }\Delta v_{\mathrm{iy}}}\\ \multicolumn{1}{c}{\frac{s}{\omega }\Delta v_{\mathrm{iz}}}\end{array}\right]}\end{array}$$

(8)

Assuming that the maximum angle of field view ${\alpha }_{\mathrm{Fov}}$ is applied, the position of chaser spacecraft $\mathbf{r}(t)$ needs to satisfy the constraint as follows

$$\left|\frac{y}{x}\right|\le \tan \frac{{\alpha }_{F0V}}{2}$$

(9)

We can obtain the boundary of impulse from equation (8). The trajectory of the chaser spacecraft will satisfy the constraint of field view attached to the target spacecraft if the impulse range is limited.

The nonlinear expression of the angle $\alpha =\arctan 2(y,x)$ makes it difficult to solve the maximum angle of field view ${\alpha }_{max}$ by analytical methods. Therefore, we calculate that problem with the numerical method. In this article, we apply a simple method wherein we divide the entire time equally into $m$ sections, and the angle of field view at the time $t_m$ can be expressed as $\alpha =\arctan 2(y,x)$ and ${\alpha }_{max}=max({\alpha }_{t_1},{\alpha }_{t_2},\dots ,{\alpha }_{t_m})\le (({\alpha }_{\mathrm{FOV}})/2)$ .

Impulse constraints

The impulse constraints considered here include the minimum time between imposing impulses and maximum impulses’ magnitude. The time between imposing impulses’ constraint can be modeled as

$$t_m-t_{m-1}>\tau (m=1,\dots ,N-1)$$

(10)

The maximum impulses’ magnitude constraint can be modeled as

$$max(\left|\Delta {\mathbf{v}}_0\right|,\dots ,\left|\Delta {\mathbf{v}}_{N-1}\right|)\le \Delta v_{max}$$

(11)

where $\tau$ represents the minimum time between two impulses, $\Delta v_{max}$ denotes the thruster’s maximum impulse, and $N$ is the number of impulses.

Passive safety

The crucial element of safety must be considered in designing the rendezvous approach. To ensure that in the case of guidance maneuver failure, especially the last maneuver, the chaser spacecraft disappears in a flight-forbidden zone, we adopt the passive safety mode. Thereby, the chaser spacecraft can free-fly with the speed of failure point and avoid colliding with the target spacecraft. Inherent in the design of the passive safety mode, in fact, is comprehension of the relationship between position and velocity of the chaser spacecraft at failure point. That relationship is studied in this article via the introduction of a cuboid flight-forbidden zone, ¹⁷ while accounting for the engineering application. We will assume that the lengths of the cuboid are 2a, 2b, and 2c along the coordinate axes. “Limit track” occurs when the trajectory of free-flying is tangent with the flight-forbidden zone. Based on the difference of failure point, limit track can be divided into two types, and the flight trajectory belongs to two types at the same time which is called critical trajectory, as shown in Figure 3. Type 1 is defined as the flight trajectory that is tangent with the boundary of flight-forbidden zone AB. Type 2 is defined as the trajectory that goes through the vertex of flight-forbidden zone B, located in the fourth quadrant. $P_0$ is located along the x-axis; when the chaser spacecraft maneuver failed at that point, the free-flight trajectory is tangent with the boundary of flight-forbidden-zone AB, and the tangency point is B. In that situation, the free-flight trajectory is called critical trajectory. When the failure point is located between the point $P_0$ and the boundary of flight-forbidden zone, as is shown in point $P_1$ in Figure 3, the free-flight trajectory is called type 1. When the failure point is located beyond $P_0$ , as is shown in point $P_2$ in Figure 3, the free-flight trajectory is called type 2.

Figure 3.

Cuboid flight-forbidden zone.

By substituting the state of point $P_0$ ( $y_{P_0}=0$ , ${\stackrel{·}{y}}_{P_0}=0$ ) and the state of point B ( $x=a$ , $y=-b$ , $\stackrel{·}{x}=0$ ) into the solutions of C–W equations, we arrive at the following equation

$$\{\begin{array}{c}\hfill x_p=a+2b(\sqrt{7}-3\arccos \frac{3}{4})\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_{P_0}=-2\omega b}\end{array}$$

(12)

We will assume that the position and velocity of the chaser spacecraft at its failure point are $x_P$ and ${\stackrel{·}{x}}_P$ , respectively. When $x_P<x_{P_0}$ , the critical trajectory is type 1, and the state of tangency point is $x=a$ , $y=-b\prime$ , $\stackrel{·}{x}=0$ . When we adopt the method from above, the resultant equation is

$$\{\begin{array}{c}\hfill x_P=a+2b\prime (\sqrt{7}-3\arccos \frac{3}{4})\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_{P1}=\frac{-(x_P-a)\omega }{(\sqrt{7}-3\arccos \frac{3}{4})}}\end{array}$$

(13)

When the critical trajectory is type 2, the vertex state of the flight-forbidden zone is B. This means that we can derive the following equation from C–W equation solutions

$$\begin{array}{cc}\hfill a\hfill & \hfill =x_P+4{\stackrel{·}{x}}_{P2}\sqrt{1-{(1+b\omega /2{\stackrel{·}{x}}_{P2})}^2}/\omega \hfill \\ \multicolumn{1}{c}{}& \hfill -3{\stackrel{·}{x}}_{P2}\arccos (1+b\omega /2{\stackrel{·}{x}}_{P2})/\omega \hfill \end{array}$$

(14)

Looking at equation (13), it is easy to see that the absolute value of velocity is decreased in a nonlinear way as the absolute value of position increases. We can acquire the safety rendezvous approach only if the velocity of the guidance maneuver at the failure point is less than that of the critical trajectory at the same point. That is

$${\stackrel{·}{x}}_{N-1}^{-}\le \left|{\stackrel{·}{x}}_{P2}\right|$$

(15)

Here, the superscript − indicates the state before an impulse. This means that the $x$ component of relative velocity before the last maneuver should less than $\left|{\stackrel{·}{x}}_{P2}\right|$ .

Mathematical model of multiple-impulse optimal rendezvous

The bulk of the problem concerning optimal rendezvous with multiple impulses lies in calculating the impulses’ vector and the time of every burn in a way that satisfies all path constraints and boundary requirements. The relative position ${\mathbf{r}}_0$ and velocity ${\mathit{v}}_0$ of two spacecrafts at initial time $t_0$ are known, and the relative position $r_f$ and velocity $v_f$ of two spacecrafts at terminal time $t_f$ are known too. Assuming that the times of impulses are n, and the impulse imposed at time $t_i$ is $\Delta {\mathbf{v}}_i={[{\begin{array}{ccc}\hfill \Delta v_{\mathrm{xi}}\hfill & \hfill \Delta v_{\mathrm{yi}}\hfill & \hfill \Delta v\hfill \end{array}}_{\mathrm{zi}}]}^T$ ( $i=1,2,\dots ,n$ ), the immediate relative velocity before the burn is ${\mathbf{v}}_i^{-}={[\begin{array}{ccc}\hfill v_{\mathrm{ix}}^{-}\hfill & \hfill v_{\mathrm{iy}}^{-}\hfill & \hfill v_{\mathrm{iz}}^{-}\hfill \end{array}]}^T$ , and the immediate relative velocity after the burn is ${\mathbf{v}}_i^{+}={[\begin{array}{ccc}\hfill v_{\mathrm{ix}}^{+}\hfill & \hfill v_{\mathrm{iy}}^{+}\hfill & \hfill v_{\mathrm{iz}}^{+}\hfill \end{array}]}^T$ . By calculating the relative position and velocity at time $t_{i+1}$ , and then inducing the hypotheses above into equation (6), we obtain the general formula

$${\mathbf{r}}_{i+1}={\mathbf{\Phi }}_{11}(\Delta t){\mathbf{r}}_i+{\mathbf{\Phi }}_{12}(\Delta t){\mathbf{v}}_i^{+}$$

(16)

$${\mathit{v}}_{i+1}^{-}={\mathbf{\Phi }}_{21}(\Delta t){\mathit{r}}_i+{\mathbf{\Phi }}_{22}(\Delta t){\mathit{v}}_i^{+}$$

(17)

Next, we can calculate the impulse using $\Delta v_i=v_i^{+}-v_i^{-}$ from equations (16) and (17).

The problem of fuel-optimal rendezvous based on C–W equations can be described as follows: calculate $(t_i,{\mathbf{X}}_i,\Delta {\mathbf{v}}_i)$ to acquire the minimum optimal index $J=\Delta v={\displaystyle \sum _{i=0}^{N-1}}\left|\Delta {\mathbf{v}}_i\right|$ . The value of $(t_i,{\mathbf{X}}_i,\Delta {\mathbf{v}}_i)$ must satisfy the state constraints and initial and terminal conditions. It can be expressed by equations as

$$\begin{array}{cc}\multicolumn{1}{c}{}& minJ=\Delta v={\displaystyle \sum _{i=0}^{N-1}}\left|\Delta {\mathbf{v}}_i\right|\hfill \\ \multicolumn{1}{c}{}& s.t.\mathbf{X}(t_0)={\mathbf{X}}_0\hfill \\ \multicolumn{1}{c}{}& \mathbf{X}(t_f)={\mathbf{X}}_f\hfill \\ \multicolumn{1}{c}{}& \left|\frac{y}{x}\right|-\tan \frac{{\alpha }_{F0V}}{2}\le 0\hfill \\ \multicolumn{1}{c}{}& t_m-t_{m-1}-\tau >0(m=1,\dots ,N-1)\hfill \\ \multicolumn{1}{c}{}& max(\left|\Delta {\mathbf{v}}_0\right|,\dots ,\left|\Delta {\mathbf{v}}_{N-1}\right|)-\Delta v_{max}\le 0\hfill \\ \multicolumn{1}{c}{}& {\stackrel{·}{x}}_{N-1}^{-}\le \left|{\stackrel{·}{x}}_{P2}\right|\hfill \end{array}$$

(18)

where s.t. represent related constraints, including initial and terminal constraints, the constraints of field view, maximum ability of maneuver, minimum time interval of impulses, and passive safety.

Choose the time of burns $t_i$ as design variables and calculate the impulse $\Delta {\mathit{v}}_i$ by iterative equations (14) and (15). Accordingly, the optimization variables of the mathematical model can be demonstrated as $y={(t_1,\dots ,t_{N-1},\Delta {\mathbf{v}}_0^{\mathrm{T}},\dots ,\Delta {\mathbf{v}}_{N-2}^{\mathrm{T}},\Delta {\mathbf{v}}_{N-1}^{\mathrm{T}})}^{\mathrm{T}}$ .

Hybrid GA

The GA is a random search method that works by imitating the biological law of evolution from nature, ¹⁸ which is “survival of the fittest, extinction of the unfitness.” Advantages of the GA include limitless continuity and derivation of function, intrinsic implicit parallelism and better global optimization ability, the search space of automatic acquisition and guide optimization, and adaptive adjusting of the search direction. However, the GA also entails some disadvantages, such as poor local search ability. The optimal result in the GA search floats within a certain time range. In this article, the hybrid GA is applied to solve that problem and improve the search precision.

Optimization-based, the hybrid GA takes SQP behind the GA to improve the poor local search ability. ¹⁹ The SQP algorithm is based on local search, with the performance of the search depending on neighborhood function and initial solutions. If the neighborhood function is unsuitably designed or the initial solutions chosen are unfit, the algorithm will be trapped into the local minimum. But the advantage of stronger local search ability exists in the SQP algorithm. So that the algorithms can learn from each other; the hybrid GA contains the advantages of both the GA (better global search ability) and the SQP (higher search precision).

1. Step 1. Solve the mathematical model of the nonlinear constraints optimal problem by a GA of continuous coding. Because solving this optimal problem requires extensive analysis as well as calculation of objective function and constraints by the GA, an analytical trajectory model is applied to lower computational cost. This model involves calculating the optimal rendezvous model with the solutions of C–W equations. For N impulses rendezvous, the second, third, and N − 1th impulses are determined after the first impulse, and the time for those impulses is given. Only the first impulse and the time for those impulses are free variables. Hence, we choose $\mathbf{z}={(t_1,\dots ,t_{N-1},\Delta {\mathbf{v}}_0^{\mathrm{T}},\dots ,\Delta {\mathbf{v}}_{N-3}^{\mathrm{T}})}^{\mathrm{T}}$ as the optimization variables of the GA and $J=\Delta v={\displaystyle \sum _{i=0}^{N-1}}\left|\Delta {\mathbf{v}}_i\right|$ as fitness function. The minimum time interval of impulses’ constraint $t_m-t_{m-1}-\tau >0(m=1,\dots ,N-1)$ is a kind of linear constraint, while the maximum ability of maneuver constraint $max(\left|\Delta {\mathbf{v}}_0\right|,\dots ,\left|\Delta {\mathbf{v}}_{N-1}\right|)-\Delta v_{max}\le 0$ , constraint of field view $\left|\frac{y}{x}\right|-\tan \frac{{\alpha }_{F0V}}{2}\le 0$ , and constraint of passive safety $\left|\Delta {\mathbf{v}}_{\mathrm{N}-1}\right|\le \Delta v_{\mathrm{safe}}$ are nonlinear constraints. We will get the optimum individual as a preliminary optimization by the GA function in the library functions of MATLAB.

2. Step 2. Set the solutions of the GA as initial solutions of the SQP. Solve the mathematical model of the nonlinear constraints’ optimal problem again by the SQP. Adopt the numerical integration of the C–W equations and then acquire precise solutions with the SQP. The SQP solver used in this study is MATLAB’s fmincon (Figure 4).

Figure 4.

The solving framework of hybrid genetic algorithm.

The control method of trajectory tracking

By using the hybrid GA above to solve the mathematical model of optimal rendezvous with a series of constraints, we can acquire the optima trajectory with linearized orbit model (equation (2)) without disturbances. And that trajectory is called nominal trajectory. However, during the on-orbit operation, the spacecrafts are influenced by many disturbances, including the linearization error of C–W equations, orbit perturbations, navigation deviation, and control error. The trajectory with optima control under disturbances is called practical trajectory. Therefore, we need to add new maneuvering controls on the original optima impulses of practical trajectory to eliminate the effects of disturbances. By doing so, we ensure that the chaser can satisfy the terminal state requirement of the rendezvous mission. In this article, we employ a tracking control method that corrects the deviations between practical trajectory and nominal trajectory with the velocity impulses offered by the control system. That control method can ensure that the chaser runs along the nominal trajectory.

Control principle

The tracking control method used in the close-range rendezvous phase is to track the nominal trajectory. The trajectories are updated continually based on the practical trajectory at the present maneuver time and the nominal trajectory at the next maneuver time to decrease the deviation between practical trajectory and nominal trajectory. Those trajectories are connected together to form the reference trajectory. As shown in Figure 5, a trajectory spanning from initial point ${\mathbf{r}}_0$ to hovering point ${\mathbf{r}}_f$ is adopted in this article. Where the maneuver time at ${\mathbf{r}}_0$ is $t_0$ , ${\mathbf{r}}_1$ is the position to which corrected impulse is applied on the time, ${\mathbf{c}}_1$ is the nominal position that makes sure that less deviation is obtained between practical and design nominal trajectory, and the maneuver times at ${\mathbf{r}}_1$ and ${\mathbf{c}}_1$ are, respectively, $t_1$ and $t_{c1}$ . At $t_0$ , making the nominal state on $t_1$ the control objective, we calculate the transfer impulse and then generate the reference trajectory based on practical state.

Figure 5.

The schematic diagram of control method for trajectory tracking.

However, because the spacecrafts are influenced by the linearization error of C–W equations, orbit perturbations, navigation deviation, control error, and so on, the practical trajectory does not move along the reference trajectory. This means that the chaser spacecraft will not arrive at position ${\mathbf{r}}_1$ but at position ${\mathbf{r}\mathbf{\text{'}}}_i$ . Using the nominal state at the position ${\mathbf{c}}_1$ as control objective, we calculate the impulse of transfer and then generate the new reference trajectory based on the practical state at the position ${\mathbf{r}\mathbf{\text{'}}}_i$ . In the same way, we can determine the corrected impulse at time $t_{c1}$ .

Control algorithm

In order to analyze the influences of the linearization error of the C–W equations along with those of orbit perturbations, navigation deviation, and control error on motion trajectory, and to realize the trajectory tracking perfectly, we have created the flow diagram shown in Figure 6.

Figure 6.

The flow diagram of control for trajectory tracking.

From the flow diagram above, we can see that we acquire the corrected impulse by navigational position and velocity directly through the tracking control method, disregarding the navigation deviation and control deviation. We can calculate the practical state of the chaser spacecraft by employing accurate relative motion equations with perturbations and then analyzing the influences of navigation deviation and control deviation. In the control algorithm of the closing phrases of rendezvous, we have adopted relative motion equations along with first-order analytical solutions of the C–W equations and the accurate relative motion equations with perturbations.

Simulation and analysis

Results and analysis of three-impulse rendezvous simulation

For our calculations, we assume that the target spacecraft runs on a 400-km-high circular orbit. This example addresses the problem of rendezvous at short distance. The facility of a spacecraft is usually confined by FOV during this stage. In this example, therefore, we study a problem of the optimal rendezvous approach with three impulses in plane. Having adopted the optimization model (18), we consider the constraints of field view, maximum maneuver ability, and minimum time interval of impulses and passive safety. Additionally, we seek the four optimization variables $t_1$ , $t_2$ , $\Delta v_{0x}$ , and $\Delta v_{0y}$ to achieve minimum fuel cost. Where $t_1$ represents the time of the second impulse, $t_2$ represents the time of the third impulse, and $\Delta v_{0x}$ and $\Delta v_{0y}$ represent the components of the first impulse along the x- and y-axes separately.

The first impulse is imposed at the initial time represented as $t_0=0$ . The initial relative position is ${\mathbf{r}}_0=[-2000,-500,0]\mathrm{m}$ ; the initial relative velocity is ${\mathbf{v}}_0=[-0.4,-0.2,0]\mathrm{m}/\mathrm{s}$ ; the terminal position of rendezvous is $r_2=[-200,0,0]\mathrm{m}$ ; the terminal velocity of rendezvous is ${\mathbf{v}}_2=[0,0,0]\mathrm{m}/\mathrm{s}$ ; and the three impulses are imposed at the positions of ${\mathbf{r}}_0$ , ${\mathbf{r}}_1$ , and ${\mathbf{r}}_2$ separately. The approach phase is finished in 4000 s, which means $t_2\le 4000\mathrm{s}$ .

With equation (8) applied in the rendezvous constraint of field view ${\alpha }_{\mathrm{FOV}}/2={40}^{\circ }$ , the minimum time interval of impulses is $t_{min}=100\mathrm{s}$ and the maximum ability of maneuver is $\Delta v_{max}=2\mathrm{m}/\mathrm{s}$ . When $x_P>x_{P_0}$ , the critical trajectory is type 2, based on the task requirements above. When considering equation (13), it is easy to see that the absolute value of velocity $\left|{\stackrel{·}{x}}_P\right|$ is decreased in a nonlinear manner with the increase in the absolute value of position $\left|x_P\right|$ . Taking this into account, it is only when the velocity of the guidance maneuver at the failure point is less than the velocity of critical trajectory at the same point that we can acquire the safety rendezvous approach. By substitution of $x_P=200$ and the length of rectangular flight-forbidden zone $a=20\mathrm{m}$ and $b=10\mathrm{m}$ into equation (15), the absolute value of velocity $\left|{\stackrel{·}{x}}_{P2}\right|$ is calculated to be 3.5777 m/s. This means that $x$ component of relative velocity at terminal time $t_f$ should be less than $3.5777\mathrm{m}/\mathrm{s}$ . In this way, the chaser spacecraft will not fly into the flight-forbidden zone after the failure of last impulse. The rendezvous approach is passive safety. This method ensures a safe rendezvous approach.

The hybrid GA is applied in the problem of optimal rendezvous at a short distance in this example. Where the number of population is 100, the maximum generation is 1000, real code compilation is applied, the crossover is 0.8, and the mutation is 0.1. In order to test the effectiveness and applicability of the hybrid GA, we compare the performance of the GA, the SQP, and the hybrid GA with the random initial state. Table 1 denotes the results of this simulation.

Table 1.

The results of different optimization algorithm.

Algorithm	$t_1(\mathrm{s})$	$t_2(\mathrm{s})$	$J(\mathrm{m}/\mathrm{s})$	$\Delta v_0(\mathrm{m}/\mathrm{s})$	$\Delta v_1(\mathrm{m}/\mathrm{s})$	$\Delta v_2(\mathrm{m}/\mathrm{s})$
GA	254.6	1150.7	3.6827	1.9992	0.1462	1.5373
SQP (random)	DNC a	DNC	DNC	DNC	DNC	DNC
SQP (appropriate)	1572.1	3699.0	1.5554	1.3196	0.1156	0.1203
Hybrid GA	1407.6	3.509.9	1.5510	1.3222	0.0958	0.1329

GA: genetic algorithm; SQP: sequential quadratic programming.

a

DNC: did not converge. The optimization method did not converge to an optimal solution.

An appropriate initial state needs to be selected when the SQP algorithm is applied to the optimization mathematical model. Otherwise, the matrix of Hessian is not positive definite and cannot update anymore. That will result in non-convergence and failure of optimization. A random initial state $(t_1,t_2,\Delta v_{0x},\Delta v_{0y})=(0,0,1,1)$ and an appropriate initial state $(t_1,t_2,\Delta v_{0x},\Delta v_{0y})=(500,1500,1,1)$ are selected separately with SQP algorithm optimization in Table 1. The results of simulation reveal that using the GA to optimize the mathematical model is the only thing that can achieve optimal results rapidly and prove the better global optimization ability. However, by using the comparison of the traditional optimization algorithm, we find that the precision of optimization results is poorer when using the SQP algorithm. In other words, the local search ability of the GA is weaker than that of the SQP algorithm. Using the hybrid GA can combine the advantages of the GA (better global search ability) and the SQP (higher search precision), thereby solving the problem of optimization rendezvous approach.

Figure 7 describes the trajectory of optimization rendezvous in plane with the results of the hybrid GA optimization. From the figure, we can see that the chaser spacecraft starts at the position of $(-2000,-500)\mathrm{m}$ and, through three impulses imposed, achieves the terminal position of $(-200,0)\mathrm{m}$ .

Figure 7.

The chaser approaches the target in plane.

Figure 8 denotes the time history of view angle during rendezvous. From this figure, we can see that the view angle exists within the range of 40° at all times. It corresponds with the field view constraint of the chaser spacecraft.

Figure 8.

Time history of view angle during rendezvous.

To demonstrate the advantage of the hybrid GA, we will compare it with linear quadratic regulator (LQR). The LQR is an optimal control approach valid for a linear state-space model. The LQR quadratic cost function can be found in Zhang et al. ²⁰ The state gain matrix and control gain matrix in the LQR are set as follows to obtain an arrival time equivalent to the hybrid GA

$$\{\begin{array}{c}\hfill \mathit{Q}=\mathrm{diag}([\begin{array}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \end{array}])\hfill \\ \multicolumn{1}{c}{\mathit{R}=\mathrm{diag}([\begin{array}{cc}\hfill {10}^{10}\hfill & \hfill {10}^{10}\hfill \end{array}])}\end{array}$$

(19)

The controller gain matrix $\mathit{K}$ can be obtained by solving the Riccati equation. The optimal feedback control is given by the well-known expression

$$\mathit{U}=-K(X_n-X)$$

(20)

The chaser approaches the target in plane, and time history of cost of the application of the LQR is presented in Figures 9 and 10. As shown in Figure 10, the $\Delta V$ cost for rendezvous with the LQR method is 6.7039 m/s. It is found that the hybrid GA can significantly reduce velocity increment compared with the LQR method.

Figure 9.

The chaser approaches the target in plane of LQR method.

Figure 10.

Time history of cost of LQR method.

Results and analysis of tracking control simulation

Considering the influences of the linearization error of the C–W equations, orbit perturbations, navigation deviation, control error, and so on, the tracking control method described above is employed to correct impulses of optimal trajectory, with the aforementioned simulation results as the control objective. Add the corrected impulses on original impulses of nominal trajectory at times $t_0,t_1,t_2$ , and $t_{c1}$ , respectively, where the choice of ${\mathit{c}}_1$ can follow experience. In this article, we chose $t_{c1}=3000\mathrm{s}$ . The navigation deviations are as follows: the error of relative position is selected as $\Delta x=5\mathrm{m}$ and $\Delta y=3\mathrm{m}$ , and the error of relative velocity is chosen as $\Delta v_e=1\times {10}^{-3}\mathrm{m}/\mathrm{s}$ . Accurate relative motion equations are selected with the accurate mathematical model, including the linearization error of the C–W equations and the error of Earth oblateness $J_2$ . The whole simulation of the closed-loop program is set up with MATLAB 6.5. The results of the simulation are shown in Figures 9 and 10.

Figure 11 shows the difference between practical trajectory and nominal trajectory. It is observed from Figure 9 and the partial enlarged view that the difference between practical trajectory and nominal trajectory is extremely small. The practical trajectory can track the nominal trajectory well.

Figure 11.

The difference between practical trajectory and nominal trajectory.

Figure 12 displays the corrected relative velocity of the chaser spacecraft varying with time. From Table 2, we can see that the overall fuel consumption with velocity impulses correcting is more than that without velocity impulses correcting.

Figure 12.

Time history of the relative velocity after correcting.

Table 2.

The velocity impulses contrast between before and after correcting for the tracking control.

Applied time	Velocity impulse (m/s)	Before correcting	After correcting
$t_0(\mathrm{s})$ 0	$\Delta v_1$	1.3222	1.3234
$t_1(\mathrm{s})$ 1407.6	$\Delta v_2$	0.0958	0.0960
$t_{c1}(\mathrm{s})$ 3000.0 s	$\Delta v_3$	0	0.0200
$t_2(\mathrm{s})$ 3509.9	$\Delta v_4$	0.1329	0.1452
	${\displaystyle \sum _{i-0}^4}\Delta v_i$	1.5510	1.5846

Conclusion

In this article, we have studied the problem of multiple-impulse, fuel-optimal rendezvous with the constraints of field view, maximum ability of maneuver, minimum time interval of impulses, and passive safety. Previous research and analysis on this topic have focused on testing whether the passive safety constraint was satisfied by a design that did not account for safety of the trajectory. In this article, however, we have calculated the maximum value of the last impulse before the trajectory design and used that value as a constraint for this trajectory design. Employing this tactic, we have found that the optimal algorithm can be easily solved.

The mathematical model of optimal planning approaches of multiple impulses with path constraints has been constructed based first on C–W equations. Trajectory planning has then been solved by application of the hybrid GA. We have shown that this algorithm contains both the better global searching ability of the GA and the stronger local searching ability of the SQP.

Finally, considering the perturbations and errors existing in the practical problem, we have applied a kind of tracking control method. We have given simulation examples with three impulses to acquire an optimal rendezvous trajectory with path constraints; the simulation results have proven the effectiveness and applicability of the mathematical model and algorithm. At the same time, our simulation has verified the applicability of the tracking control method. However, the order the errors, including navigation, force model, and/or linearization, cause position difference at the terminal time that needs an additional, non-designed impulse to eliminate is not considered in this study. This will be studied in our future work.

This article has shown that when the tracking control method is used, the errors between the practical trajectory and nominal trajectory are corrected perfectly with an additional impulse, and the chaser spacecraft arrives at terminal position accurately and safely. Overall, this article has comprehensively addressed trajectory planning as well as tracking control of the close-range rendezvous mission.