Nomenclature

Abbreviations

AU

Unit acceleration

CW

Clohessy-Wiltshire

DU

Unit distance

GA

Genetic algorithm

GA-PSO

Genetic algorithm-particle swarm optimisation

NLP

Nonlinear programming

PSO

Particle swarm optimisation

SA

Simulated annealing

SQP

Sequential quadratic programming

TU

Unit time

Symbols

$\left( {{r_i},{\xi _i},{\phi _i},{\gamma _i},{v_i},{\zeta _i}} \right)$

state variables of spacecraft

$\left( {{a_i},{e_i},{i_i},{\Omega _i},{\omega _i},{f_i}} \right)$

six orbital elements

${O_e} - {x_e}{y_e}{z_e}$

inertial coordinate system

${\boldsymbol{R}_{{\zeta _i}}}{\boldsymbol{R}_{{\phi _i}}}{\boldsymbol{R}_{{\xi _i}}}\ {\boldsymbol{R}_{{\theta _i}}}{\boldsymbol{R}_{{i_i}}}{\boldsymbol{R}_{{\Omega _i}}}$

matrix of coordinate transformation

$O - {x_O}{y_O}{z_O}$

orbit coordinate system of spacecraft

$\mu $

gravitational coefficient of the earth

${\hat{\boldsymbol{c}}_1}$ , ${\hat{\boldsymbol{c}}_2}$ , ${\hat{\boldsymbol{c}}_3}$

coordinate base vector of ${O_e} - {x_e}{y_e}{z_e}$

${J_2}$

perturbation coefficient

$\boldsymbol{r}$ , $\hat{\boldsymbol{\theta }}$ , $\boldsymbol{h}$

coordinate base vector of $O - {x_O}{y_O}{z_O}$

${g_0}$

standard gravitational acceleration

$T$

value of thrust on spacecraft

${I_s}$

specific impulses of thrust

$\dot m$

propellant mass flow rate

$\beta $

angle between thrust vector $\boldsymbol{T}$ and its projection in $O - {x_O}{y_O}$ plane

$\alpha $

angle between the projection of thrust vector $\boldsymbol{T}$ in $O - {x_O}{y_O}$ plane and velocity vector $\boldsymbol{v}$

$\zeta $

angle between the projection of velocity vector $\boldsymbol{v}$ in plane $O - {y_L}{z_L}$ and axis $O{y_L}$

$\gamma $

angle between the velocity vector $\boldsymbol{v}$ and its projection in plane $O - {y_L}{z_L}$

${\boldsymbol{u}_E}$

control actions of the pursuer and the target

${\boldsymbol{u}_P}$

control actions of the pursuer

${\boldsymbol{x}_E}$

state vector of the target

${\boldsymbol{x}_P}$

state vector of the pursuer

${\dot{\boldsymbol{x}}_E}$

state differential vector of the target

${\dot{\boldsymbol{x}}_P}$

state differential vector of the pursuer

$J$

performance index function

${\boldsymbol{\psi }}( \cdot )$

set of boundary constraints

$H$

Hamiltonian function

$\left( {\boldsymbol{u}_P^*,\boldsymbol{u}_E^*} \right)$

saddle point solution

${\boldsymbol{f}_P}()$ ${\boldsymbol{f}_E}()$

state equations of pursuer and target

${\rm{\Phi }}$

terminal conditions of the system

${{\boldsymbol{\lambda }}_E}$

covariates corresponding to ${\boldsymbol{f}_E}()$

${{\boldsymbol{\lambda }}_P}$

covariates corresponding to ${\boldsymbol{f}_P}()$

${k_{{m_P}}}$

weight coefficient of the fuel consumption of the pursuer

${k_t}$

weight coefficient of interception time

${t_f}$

interception time

${t_0}$

initial time of manoeuvring

${T_E}$

value of thrust on target

${T_P}$

value of thrust on pursuer

${m_E}$

mass of target

${m_P}$

mass of pursuer

${\boldsymbol{\upsilon }}$

covariates corresponding to the boundary constraint conditions

${I_{{s_P}}}$

specific impulses of pursuer

${u_1}\sim {u_4}$

components of ${\boldsymbol{u}_P}$ and ${\boldsymbol{u}_E}$

${\lambda _1}\sim {\lambda _{12}}$

components of ${{\boldsymbol{\lambda }}_P}$ and ${{\boldsymbol{\lambda }}_E}$

$u_1^{\left( A \right)}u_1^{\left( B \right)}\sim u_4^{\left( A \right)}u_4^{\left( B \right)}$

the possible values of ${u_1}\sim {u_4}$

${x_1}\sim {x_{12}}$

components of ${\boldsymbol{x}_P}$ and ${\boldsymbol{x}_E}$

${\eta _i}$

decision variables

${k_i}$

weight value of $J$

${h_{safe}}$

safe altitude

$\varepsilon $

specified small value

${R_E}$

Earth radius

$g$

Gravitational acceleration

$\mu $

Gravitational constant

${c_i}$

distance from the set of functions composed of Equations (26), (30) and (34)

2.0 Transformation between state variables and six elements of orbits

For constant continuous small thrust engine, the control variable is the azimuth angle of thruster nozzle, in order to solve this problem conveniently, in the orbit coordinate system, a new set of state variables $\left( {{r_i},{\xi _i},{\phi _i},{\gamma _i},{v_i},{\zeta _i}} \right)$ is defined to describe the motion state of spacecraft, ${r_i}$ , ${\gamma _i}$ , ${v_i}$ , ${\zeta _i}$ are geocentric distance, flight path angle, flight speed, azimuth of spacecraft, respectively, ${\xi _i}$ and ${\phi _i}$ are geographical longitude and latitude of spacecraft, respectively.

The purpose of defining the state variable is to obtain the real-time change of the control angle of the engine thruster relative to its own spacecraft, which is more intuitive. If the state variable in inertial coordinate system is used, only the inertial component of thrust vector can be obtained. The geocentric equatorial inertial coordinate system ${O_e} - {x_e}{y_e}{z_e}$ is established with the geocentric as the coordinate origin ${O_e}$ . The right-hand system is composed of ${x_e}$ , ${y_e}$ and ${z_e}$ . ${O_e}{x_e}$ represents the direction which the North Pole pointing to the earth’s rotation axis, and ${O_e}{z_e}$ represents the direction pointing to the vernal equinox. The orbit coordinate system of spacecraft (represented by $O - {x_O}{y_O}{z_O}$ ) is established with centre of mass as the coordinate origin (represented by $O$ ), in which the direction of the radial vector $\boldsymbol{r}$ is the axis $o{x_O}$ , the axis $o{z_O}$ is the orbital angular momentum direction, the axis $o{y_O}$ is the direction of $\hat{\boldsymbol{\theta }}$ , and the three axis constitute the right-hand system. Another coordinate system adopted in our study is the local horizontal coordinate system $O - {x_L}{y_L}{z_L}$ , and the coordinate origin of $O - {x_L}{y_L}{z_L}$ is the same with the one of $O - {x_O}{y_O}{z_O}$ ; the ${y_L}$ axis points east in the local horizontal plane; the radial vector direction of ${O_e}$ pointing to $O$ is the ${x_L}$ -axis; the ${z_L}$ -axis points north in the local horizontal plane.

Figures 4 and 5 are the orbit coordinate system of spacecraft, which respectively described by state variables $\left( {{r_i},{\xi _i},{\phi _i},{\gamma _i},{v_i},{\zeta _i}} \right)$ and six orbital elements $\left( {{a_i},{e_i},{i_i},{\Omega _i},{\omega _i},{f_i}} \right)$ in the geocentric equatorial inertial coordinate system ${O_e} - {x_e}{y_e}{z_e}$ . Where, ${\hat{\boldsymbol{c}}_1}$ , ${\hat{\boldsymbol{c}}_2}$ and ${\hat{\boldsymbol{c}}_3}$ are the coordinate base vector of ${O_e} - {x_e}{y_e}{z_e}$ , and $\boldsymbol{r}$ , $\hat{\boldsymbol{\theta }}$ and $\boldsymbol{h}$ are the coordinate base vector of $O - {x_O}{y_O}{z_O}$ , and $N$ is the intersection line between the orbital plane and the equatorial plane.

Figure 4. The orbit coordinate system represented by state variables.

Figure 5. The orbit coordinate system represented by orbital elements.

According to Fig. 4, the coordinate base vector of $O - {x_O}{y_O}{z_O}$ can be obtained by three times coordinate transformations, which is given as follows:

(1) \begin{align}\left[ \begin{array}{l}\boldsymbol{r}\\{\hat{\boldsymbol{\theta }}}\\\boldsymbol{h}\end{array} \right] & = {\boldsymbol{r}_{{\zeta _i}}}{\boldsymbol{R}_{{\phi _i}}}{\boldsymbol{R}_{{\xi _i}}}\left[ \begin{array}{l}{{\hat{\boldsymbol{c}}}_1}\\{{\hat{\boldsymbol{c}}}_2}\\{{\hat{\boldsymbol{c}}}_3}\end{array} \right] \nonumber\\& = \left[ {\begin{array}{ccccccc}{\cos {\phi _i}\cos {\xi _i}} & & & {\cos {\phi _i}\sin {\xi _i}} & & & {\sin {\phi _i}}\\{ - \sin {\zeta _i}\sin {\phi _i}\cos {\xi _i} - \cos {\zeta _i}\sin {\xi _i}}& & & { - \sin {\zeta _i}\sin {\phi _i}\sin {\xi _i} + \cos {\zeta _i}\cos {\xi _i}}& & & {\sin {\zeta _i}\cos {\phi _i}}\\{ - \cos {\zeta _i}\sin {\phi _i}\cos {\xi _i} + \sin {\zeta _i}\sin {\xi _i}} & && { - \cos {\zeta _i}\sin {\phi _i}\sin {\xi _i} - \sin {\zeta _i}\cos {\xi _i}} & && {\cos {\zeta _i}\cos {\phi _i}}\end{array}} \right]\left[ \begin{array}{l}{{\hat{\boldsymbol{c}}}_1}\\{{\hat{\boldsymbol{c}}}_2}\\{{\hat{\boldsymbol{c}}}_3}\end{array} \right]\end{align}

Similarly, according to Fig. 5, the vector ${\left[ {\begin{array}{*{20}{c}}\boldsymbol{r} & {\hat{\boldsymbol{\theta }}} & \boldsymbol{h}\end{array}} \right]^T}$ can also be obtained by three times coordinate transformations as follows:

(2) \begin{align} & \left[ \begin{array}{l} \boldsymbol{r}\\{\hat{\boldsymbol{\theta }}}\\\boldsymbol{h}\end{array} \right] = {\boldsymbol{R}_{{\theta _i}}}{\boldsymbol{R}_{{i_i}}}{\boldsymbol{R}_{{\Omega _i}}}\left[ \begin{array}{l}{{\hat{\boldsymbol{c}}}_1}\\{{\hat{\boldsymbol{c}}}_2}\\{{\hat{\boldsymbol{c}}}_3}\end{array} \right]\nonumber\\& = \left[ {\begin{array}{ccccccc}{\cos {\theta _i}\cos {\Omega _i} - \sin {\theta _i}\cos {i_i}\sin {\Omega _i}} &&& {\cos {\theta _i}\sin {\Omega _i} + \sin {\theta _i}\cos {i_i}\cos {\Omega _i}} &&& {\sin {\theta _i}\sin {i_i}}\\{ - \sin {\theta _i}\cos {\Omega _i} - \cos {\theta _i}\cos {i_i}\sin {\Omega _i}} & &&{ - \sin {\theta _i}\sin {\Omega _i} + \cos {\theta _i}\cos {i_i}\cos {\Omega _i}} & &&{\cos {\theta _i}\sin {i_i}}\\{\sin {i_i}\sin {\Omega _i}} & &&{ - \sin {i_i}\cos {\Omega _i}} & &&{\cos {i_i}}\end{array}} \right]\left[ \begin{array}{l}{{\hat{\boldsymbol{c}}}_1}\\{{\hat{\boldsymbol{c}}}_2}\\{{\hat{\boldsymbol{c}}}_3}\end{array} \right].\end{align}

Then, ${\boldsymbol{R}_{{\zeta _i}}}{\boldsymbol{R}_{{\phi _i}}}{\boldsymbol{R}_{{\xi _i}}} = {\boldsymbol{R}_{{\theta _i}}}{\boldsymbol{R}_{{i_i}}}{\boldsymbol{R}_{{\Omega _i}}}$ , and the relation of relative angles is given as follows:

(3) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\cos {\phi _i}\cos {\xi _i} = \cos {\theta _i}\cos {\Omega _i} - \sin {\theta _i}\cos {i_i}\sin {\Omega _i}}\\{\cos {\phi _i}\sin {\xi _i} = \cos {\theta _i}\sin {\Omega _i} + \sin {\theta _i}\cos {i_i}\cos {\Omega _i}}\\{\sin {\phi _i} = \sin {\theta _i}\sin {i_i}{\rm{ }}}\\{\sin {\zeta _i}\cos \phi = \cos {\theta _i}\sin {i_i}{\rm{ }}}\\{\cos {\zeta _i}\cos {\phi _i} = \cos {i_i}{\rm{ }}}\end{array}} \right.\end{align}

The state variables ${\phi _i}$ , ${\xi _i}$ and ${\zeta _i}$ can be calculated by Equation (3). For ${r_i}$ and ${v_i}$ , it can be directly solved by the six elements of orbit as follows:

(4) \begin{align}{r_i} = \frac{{{a_i}(1 - {e_i}^2)}}{{1 + {e_i}\cos {\theta _i}}},\ v_i = \sqrt{\mu \left( \frac{2}{r_i} - \frac{1}{a_i} \right) } \end{align}

And the flight path angle can be obtained from the radial velocity component of the spacecraft:

(5) \begin{align}{v_{ri}} = \frac{\mu }{{{h_i}}}{e_i}\sin {f_i} = \sqrt {\frac{\mu }{{{p_i}}}} {e_i}\sin {f_i} = \sqrt {\frac{\mu }{{{a_i}\left( {1 - e_i^2} \right)}}} {e_i}\sin {f_i} = {\dot r_i} = {v_i}\sin {\gamma _i}\end{align}

At this point, when the number of six orbits of the spacecraft is given, all the elements in the corresponding state variables can be obtained according to the above transformation relationship. If the state variables are known, the number of six orbitals can also be obtained by the transformation method mentioned above.

Thus, the orbit motion equation of spacecraft with J2 term can be established as follows: [Reference Liu24]

(6) \begin{align} {\begin{cases}{\dot \gamma = \dfrac{{T\cos \beta \sin \alpha }}{{mv}} + \dfrac{{v\cos \gamma }}{r} - \dfrac{{\mu \cos \gamma }}{{v{r^2}}} - \dfrac{{3{G_{{J_2}}}}}{{4v{r^4}}}\left( {\left( {3\cos \left( {2\phi } \right) - 1} \right)\cos \gamma - 2\sin \left( {2\phi } \right)\cos \zeta \sin \gamma } \right)}\\[5pt]{\dot r = v\sin \gamma }\\[5pt]{\dot v = \dfrac{{T\cos \alpha \cos \beta }}{m} - \dfrac{{\mu \sin \gamma }}{{{r^2}}} - \dfrac{{3{G_{{J_2}}}}}{{4{r^4}}}\left( {\left( {3\cos \left( {2\phi } \right) - 1} \right)\sin \gamma + 2\sin \left( {2\phi } \right)\cos \zeta \cos \gamma } \right)}\\[8pt]{\dot \phi = \dfrac{{v\cos \gamma \sin \zeta }}{r}}\\[8pt]{\dot \xi = \dfrac{{v\cos \gamma \cos \zeta }}{{r\cos \phi }}}\\[10pt]{\dot \zeta = \dfrac{{T\sin \beta }}{{mv\cos \gamma }} - \dfrac{{v\tan \phi \cos \gamma \cos \zeta }}{r} + \dfrac{{3{G_{{J_2}}}}}{{2v{r^4}\cos \gamma }}\sin \left( {2\phi } \right)\sin \zeta }\\[10pt]{\dot m = - \dfrac{T}{{{g_0}{I_s}}}}\end{cases}} \end{align}

where ${G_{{J_2}}} = {J_2}\mu R_e^2$ , $\mu = 3.986 \times {10^5}{\rm{k}}{{\rm{m}}^{\rm{3}}}{\rm{/}}{{\rm{s}}^{\rm{2}}}$ is the gravitational coefficient of the earth; ${J_2} = 1.018 \times {10^{ - 3}}$ is perturbation coefficient; ${R_e} = 6378.165\,\textrm{km}$ is the average radius of the earth; $\dot m$ is the propellant mass flow rate; ${g_0} = 9.80665m/{s^2}$ is the standard gravitational acceleration; $T$ is the thrust value; ${I_s}$ is the according specific impulse; the description and definition of $\alpha $ , $\beta $ , $\gamma $ , $\zeta $ are given in Figs 6–8.

Figure 6. Illustration of $O - {x_L}{y_L}{z_L}$ .

Figure 7. Relationship between ${O_e} - {x_e}{y_e}{z_e}$ and $O - {x_L}{y_L}{z_L}$ .

Figure 8. Illustration of thrust under $O - {x_O}{y_O}{z_O}$ .

3.0 Zero-sum differential games

Assuming that both the pursuer and the target can fully know the information of dynamics and thrust acceleration of the other side, and both sides can update and optimise their own control output according to the other party’s information, then the interception of manoeuvring target can be regarded as a zero-sum differential game problem.

Assuming that the pursuer (interceptor) is a player $P$ and the target is another one (denoted by $E$ ), the control actions of the pursuer and the target are ${\boldsymbol{u}_P}$ and ${\boldsymbol{u}_E}$ , respectively. Since the two spacecraft are independent, then the dynamic model can be divided, so the state equations of the two spacecrafts are as follows:

(7) \begin{align}\begin{array}{*{20}{c}}{ \dot{\boldsymbol{x}}_P = {\boldsymbol{f}_P}\left( {{\boldsymbol{x}_P},{\boldsymbol{u}_P},t} \right)} \qquad {{\dot{\boldsymbol{x}}_E} = {\boldsymbol{f}_E}\left( {{\boldsymbol{x}_E},{\boldsymbol{u}_E},t} \right)}\end{array}\end{align}

where $t \in \left( {{t_0},{t_f}} \right)$ , and the initial time ${t_0}$ and terminal time ${t_f}$ can be fixed or free.

The set of boundary constraints is as follows:

(8) \begin{align}{\boldsymbol{\psi }}\left( {{\boldsymbol{x}_{P0}},{\boldsymbol{x}_{E0}},{\boldsymbol{x}_{Pf}},{\boldsymbol{x}_{Pf}},{t_0}{\rm{,}}{t_f}} \right) = {0}\end{align}

The performance index function is as follows:

(9) \begin{align}J = \varphi \left( {{\boldsymbol{x}_{P0}},{\boldsymbol{x}_{E0}},{\boldsymbol{x}_{Pf}},{\boldsymbol{x}_{Pf}},{t_0}{\rm{,}}{t_f}} \right)\end{align}

The pursuer wants to intercept the target as soon as possible with the least fuel consumption, that is, the minimum hope; the target wants to extend the intercepted time infinitely, have enough time to escape the pursuit, or hope that the fuel of the pursuer is exhausted in the interception process, that is, the hope is the most. The above problems can be described by the following inequalities:

(10) \begin{align}J\left( {\boldsymbol{u}_P^*,{\boldsymbol{u}_E}} \right) \le J\left( {\boldsymbol{u}_P^*,\boldsymbol{u}_E^*} \right) \le J\left( {{\boldsymbol{u}_P},\boldsymbol{u}_E^*} \right),\,\forall {\boldsymbol{u}_P},{\boldsymbol{u}_E}\end{align}

The problem belongs to the traditional zero-sum differential game, and the saddle point solution is $\left( {\boldsymbol{u}_P^*,\boldsymbol{u}_E^*} \right)$ . It can also be seen from the above equation. That if the pursuer or target adopts non optimal control strategy, the function value of performance index will deteriorate relative to the target value of saddle point solution no matter what strategy is adopted. If the target cannot be intercepted successfully when the saddle point solution is obtained, it must not be intercepted.

The corresponding pursuit trajectory can be obtained by integrating the saddle point solution into the system state equation. In order to obtain the necessary conditions of existence of $\left( {\boldsymbol{u}_P^*,\boldsymbol{u}_E^*} \right)$ , the Hamiltonian function $H$ and terminal conditions of the system ${\rm{\Phi }}$ are defined as follows:

(11) \begin{align}\begin{array}{*{20}{c}}{H = {\boldsymbol{\lambda }}_P^{\boldsymbol{T}}\,{\boldsymbol{f}_P} + {\boldsymbol{\lambda }}_E^{\boldsymbol{T}}}\,{\boldsymbol{f}_E} = {H_P} + {H_E} \qquad {{\rm{\Phi }} = \varphi + {{\boldsymbol{\upsilon }}^T}{\boldsymbol{\psi }}}\end{array}\end{align}

Where, ${{\boldsymbol{\lambda }}_P}$ and ${{\boldsymbol{\lambda }}_E}$ are the covariates corresponding to the state equation of the pursuer and the target, respectively, and ${\boldsymbol{\upsilon }}$ represents the covariates corresponding to the boundary constraint conditions. The necessary condition of existence is:

(11) \begin{align}\begin{array}{*{20}{c}}{{{\dot{\boldsymbol{\lambda }}}_P} = - {{\left[ {\dfrac{{\partial H}}{{\partial {\boldsymbol{x}_P}}}} \right]}^T} = - {{\left[ {\dfrac{{\partial {f_P}}}{{\partial {\boldsymbol{x}_P}}}} \right]}^T}{{\boldsymbol{\lambda }}_P}} \qquad {{{\dot{\boldsymbol{\lambda }}}_E} = - {{\left[ {\dfrac{{\partial H}}{{\partial {\boldsymbol{x}_E}}}} \right]}^T} = - {{\left[ {\dfrac{{\partial {f_E}}}{{\partial {\boldsymbol{x}_E}}}} \right]}^T}{{\boldsymbol{\lambda }}_E}}\end{array}\end{align}

The corresponding boundary conditions to the covariate variables ${{\boldsymbol{\lambda }}_P}$ and ${{\boldsymbol{\lambda }}_E}$ are:

(12) \begin{align}\begin{array}{l}\begin{array}{*{20}{c}}{{\lambda _{Pk}}\left( {{t_0}} \right) + \dfrac{{\partial {\rm{\Phi }}}}{{\partial {x_{Pk}}\left( {{t_0}} \right)}} = 0} & {{x_{Pk}}\left( {{t_0}} \right)free} &\quad {k = 1,2, \ldots ,{n_P}}\end{array}\\[10pt]\begin{array}{*{20}{c}}{{\lambda _{Ek}}\left( {{t_0}} \right) + \dfrac{{\partial {\rm{\Phi }}}}{{\partial {x_{Ek}}\left( {{t_0}} \right)}} = 0} & {{x_{Ek}}\left( {{t_0}} \right)free} &\quad {k = 1,2, \ldots ,{n_E}}\end{array}\end{array}\end{align}

(13) \begin{align}\begin{array}{l}\begin{array}{*{20}{c}}{{\lambda _{Pk}}\left( {{t_f}} \right) - \dfrac{{\partial {\rm{\Phi }}}}{{\partial {x_{Pk}}\left( {{t_f}} \right)}} = 0} & {{x_{Pk}}\left( {{t_f}} \right)free} &\quad {k = 1,2, \ldots ,{n_P}}\end{array}\\[10pt]\begin{array}{*{20}{c}}{{\lambda _{Ek}}\left( {{t_f}} \right) - \dfrac{{\partial {\rm{\Phi }}}}{{\partial {x_{Ek}}\left( {{t_f}} \right)}} = 0} & {{x_{Ek}}\left( {{t_f}} \right) free} &\quad {k = 1,2, \ldots ,{n_E}}\end{array}\end{array}\end{align}

The cross-sectional conditions corresponding to the two-point boundary values are as follows:

(14) \begin{align}\begin{array}{*{20}{c}}{\dfrac{{\partial \Phi }}{{\partial {t_0}}} - H\left( {{t_0}} \right) = 0} & {{t_0}\begin{array}{*{20}{c}}{} &\quad {is}\ \end{array}free}\\[10pt]{\dfrac{{\partial \Phi }}{{\partial {t_f}}} + H\left( {{t_f}} \right) = 0} & {{t_f}\begin{array}{*{20}{c}}{} &\quad {is}\ \end{array}free}\end{array}\end{align}

For the sum of control variables ${\boldsymbol{u}_P}$ and ${\boldsymbol{u}_E}$ , the Hamiltonian function is separable, so:

(15) \begin{align}\mathop {\min }\limits_{{\boldsymbol{u}_P}} \mathop {\max }\limits_{{\boldsymbol{u}_E}} H = \mathop {\max }\limits_{{\boldsymbol{u}_E}} \mathop {\min }\limits_{{\boldsymbol{u}_P}} H = \min {\boldsymbol{\lambda }}_P^T\,{\boldsymbol{f}_P} + \max {\boldsymbol{\lambda }}_E^T\,{\boldsymbol{f}_E}\end{align}

The saddle point solution $\left( {\boldsymbol{u}_P^*,\boldsymbol{u}_E^*} \right)$ must satisfy the following two equations:

(16) \begin{align}\begin{array}{*{20}{c}}{{\boldsymbol{u}_P} = \arg \mathop {\min }\limits_{{\boldsymbol{u}_P}} H = \arg \mathop {\min }\limits_{{\boldsymbol{u}_P}} \left( {{\boldsymbol{\lambda }}_P^T\,{\boldsymbol{f}_P}} \right)} & {{\boldsymbol{u}_E} = \arg \mathop {\max }\limits_{{\boldsymbol{u}_E}} H = \arg \mathop {\max }\limits_{{\boldsymbol{u}_E}} \left( {{\boldsymbol{\lambda }}_E^T\,{\boldsymbol{f}_E}} \right)}\end{array}\end{align}

When the control variables are not constrained, the Equation (13) can be replaced by the following equation:

(17) \begin{align}\begin{array}{*{20}{c}}{{{\left[ {\dfrac{{\partial H}}{{\partial {\boldsymbol{u}_P}}}} \right]}^T} = {{\left[ {\dfrac{{\partial {\boldsymbol{f}_P}}}{{\partial {\boldsymbol{u}_P}}}} \right]}^T}{{\boldsymbol{\lambda }}_P} = 0} & {{{\left[ {\dfrac{{\partial H}}{{\partial {\boldsymbol{u}_E}}}} \right]}^T} = {{\left[ {\dfrac{{\partial {\boldsymbol{f}_E}}}{{\partial {\boldsymbol{u}_E}}}} \right]}^T}{{\boldsymbol{\lambda }}_E} = 0}\end{array}\end{align}

The above equation guarantees ${\boldsymbol{u}_P}$ and ${\boldsymbol{u}_E}$ is the extreme point. In order to solve the saddle point solution, ${\boldsymbol{u}_P}$ and ${\boldsymbol{u}_E}$ should also satisfy the following second order conditions:

(18) \begin{align}{H_{{\boldsymbol{u}_P}{\boldsymbol{u}_P}}} = \frac{{{\partial ^2}H}}{{\partial \boldsymbol{u}_P^2}} \geq 0,{H_{{\boldsymbol{u}_E}{\boldsymbol{u}_E}}} = \frac{{{\partial ^2}H}}{{\partial \boldsymbol{u}_E^2}} \le 0\end{align}

That is to say, the Hessian matrix of Hamiltonian function with respect to control variable ${\boldsymbol{u}_P}$ should be semi positive definite, and the Hessian matrix with respect to ${\boldsymbol{u}_E}$ should be semi negative definite. The control variables of ${\boldsymbol{u}_P}$ and ${\boldsymbol{u}_E}$ corresponding to the minimax of Hamiltonian function respectively.

The first-order and second-order necessary conditions of optimal control can be obtained by the above derivation. The satellite interception problem is transformed into a two-point boundary value problem. The two-point boundary value problem is composed of Equations (7), (8), (11)–(14), (17) and (18). The key to solve the two-point boundary value problem is to find initial value of the covariate which satisfy terminal constraints.

4.0 Saddle point solutions of differential games

The interception time and the fuel consumption of the pursuer in the interception process are selected as the performance indexes of the pursuit and escape problem:

(19) \begin{align}J = \int_{{t_0}}^{{t_f}} {\left[ {{k_t} + {k_{{m_P}}}\left| {{{\dot m}_P}} \right|} \right]dt} = \int_{{t_0}}^{{t_f}} {\left[ {{k_t} + {k_{{m_P}}}\left| { - \frac{{{T_P}}}{{{g_0}{I_{{s_P}}}}}} \right|} \right]dt} = \left( {{k_t} + \frac{{{k_{{m_P}}}{T_P}}}{{{g_0}{I_{{s_P}}}}}} \right)\left( {{t_f} - {t_0}} \right)\end{align}

Where, ${k_t}$ is the weight coefficient of interception time and ${k_{{m_P}}}$ is the weight coefficient of the fuel consumption of the pursuer. If the pursuer and the target manoeuvre at the same time, then ${t_0} = 0$ . Given the initial states of the two sides, assuming that the manoeuvrability of the pursuer is greater than that of the target spacecraft, that is ${T_P}/{m_P} \gt {T_E}/{m_E}$ , the pursuer must be able to catch up with the target spacecraft. Then $J$ can be transformed into:

(20) \begin{align}J = \left( {{k_t} + \frac{{{k_{{m_P}}}{T_P}}}{{{g_0}{I_{{s_P}}}}}} \right){t_f}\end{align}

Since the quantities in brackets above are constant, the performance index can also be simplified as the optimisation of interception time:

(21) \begin{align}J = {t_f}\end{align}

Let the state variables of pursuer and target and the corresponding covariates be as follows:

(22) \begin{align}{x_{P\left( E \right)}} & = {\left[ {\begin{array}{*{20}{l}}{{x_{1\left( 7 \right)}}} & {{x_{2\left( 8 \right)}}} & {{x_{3\left( 9 \right)}}} & {{x_{4\left( {10} \right)}}} & {{x_{5\left( {11} \right)}}} & {{x_{6\left( {12} \right)}}}\end{array}} \right]^T}\nonumber\\& = {\left[ {\begin{array}{*{20}{l}}{{r_{P\left( E \right)}}} & {{v_{P\left( E \right)}}} & {{\gamma _{P\left( E \right)}}} & {{\xi _{P\left( E \right)}}} & {{\phi _{P\left( E \right)}}} & {{\zeta _{P\left( E \right)}}}\end{array}} \right]^T}\end{align}

(23) \begin{align}\begin{array}{*{20}{c}}{{{\boldsymbol{\lambda }}_P} = {{[{\lambda _1},{\lambda _2},{\lambda _3},{\lambda _4},{\lambda _5},{\lambda _6}]}^T}} \qquad {{{\boldsymbol{\lambda }}_E} = {{[{\lambda _7},{\lambda _8},{\lambda _9},{\lambda _{10}},{\lambda _{11}},{\lambda _{12}}]}^T}}\end{array}\end{align}

By changing the angle of the thruster $\alpha $ and $\beta $ , the pursuer and the target can adjust the pursuit trajectory, so ${u_P}$ and ${u_E}$ are represented respectively as follows:

(24) \begin{align}\begin{array}{*{20}{c}}{{u_P} = {{\left[ {\begin{array}{ccc}{{\alpha _P}} & & {{\beta _P}}\end{array}} \right]}^T} = {{\left[ {{u_1}\quad {u_2}} \right]}^T}} \qquad {{u_E} = {{\left[ {\begin{array}{*{20}{c}}{{\alpha _E}} & &{{\beta _E}}\end{array}} \right]}^T} = {{\left[ {{u_3} \quad {u_4}} \right]}^T}}\end{array}\end{align}

For the problem of pursuiting and escaping with given initial boundary conditions, the terminal boundary conditions can be obtained according to the end conditions of interception. Assuming that the fuel of the two spacecraft is sufficient, the following conditions should be met by the pursuer and the target at the interception time ${t_f}$ , that is, ${\left[ {\begin{array}{*{20}{c}}r & \xi & \phi \end{array}} \right]^T}$ of the two spacecraft is the same:

(25) \begin{align}{r_P} = {r_E}\left| {_{t = {t_f}}} \right.,{\xi _P} = {\xi _E}\left| {_{t = {t_f}}} \right.,{\phi _P} = {\phi _E}\left| {_{t = {t_f}}} \right.\end{align}

The set of terminal boundary constraints can be obtained as follows:

(26) \begin{align}{\boldsymbol{\psi }} = \left[ \begin{array}{l}{x_{1f}} - {x_{7f}}\\{x_{4f}} - {x_{10f}}\\{x_{5f}} - {x_{11f}}\end{array} \right] = {\textit{0}}\end{align}

The ccovariates corresponding to the boundary set are:

\begin{align*}{\boldsymbol{\upsilon }} = {[{\upsilon _1},{\upsilon _2},{\upsilon _3}]^T}\end{align*}

For $H$ , the necessary and sufficient conditions of minimax algorithm are applied to obtain equation of the covariate as follows:

(27) \begin{align}\begin{array}{*{20}{c}}{{{\dot{\boldsymbol{\lambda }}}_P} = - {{\left[ {\dfrac{{\partial H}}{{\partial {\boldsymbol{x}_P}}}} \right]}^T} = - {{\left[ {\dfrac{{\partial {f_P}}}{{\partial {\boldsymbol{x}_P}}}} \right]}^T}{{\boldsymbol{\lambda }}_P}} \qquad {{{\dot{\boldsymbol{\lambda }}}_E} = - {{\left[ {\dfrac{{\partial H}}{{\partial {\boldsymbol{x}_E}}}} \right]}^T} = - {{\left[ {\dfrac{{\partial {f_E}}}{{\partial {\boldsymbol{x}_E}}}} \right]}^T}{{\boldsymbol{\lambda }}_E}}\end{array}\end{align}

(28) \begin{align}\left[ {\begin{array}{*{20}{l}}{{{\dot \lambda }_1}}\\{{{\dot \lambda }_2}}\\ \vdots \\{{{\dot \lambda }_{11}}}\\{{{\dot \lambda }_{12}}}\end{array}} \right] = - {\left[ {\begin{array}{*{20}{l}}{\dfrac{{\partial {{\dot x}_1}}}{{\partial {x_1}}}} & {\dfrac{{\partial {{\dot x}_1}}}{{\partial {x_2}}}} & \ldots & {\dfrac{{\partial {{\dot x}_1}}}{{\partial {x_{11}}}}} & {\dfrac{{\partial {{\dot x}_1}}}{{\partial {x_{12}}}}}\\[10pt]{\dfrac{{\partial {{\dot x}_2}}}{{\partial {x_1}}}} & {\dfrac{{\partial {{\dot x}_2}}}{{\partial {x_2}}}} & \ldots & {\dfrac{{\partial {{\dot x}_2}}}{{\partial {x_{11}}}}} & {\dfrac{{\partial {{\dot x}_2}}}{{\partial {x_{12}}}}}\\[2pt] \vdots & \vdots & \ddots & \vdots & \vdots \\[2pt] {\dfrac{{\partial {{\dot x}_{11}}}}{{\partial {x_1}}}} & {\dfrac{{\partial {{\dot x}_{11}}}}{{\partial {x_2}}}} & \ldots & {\dfrac{{\partial {{\dot x}_{11}}}}{{\partial {x_{11}}}}} & {\dfrac{{\partial {{\dot x}_{11}}}}{{\partial {x_{12}}}}}\\[10pt] {\dfrac{{\partial {{\dot x}_{12}}}}{{\partial {x_1}}}} & {\dfrac{{\partial {{\dot x}_{12}}}}{{\partial {x_2}}}} & \ldots & {\dfrac{{\partial {{\dot x}_{12}}}}{{\partial {x_{11}}}}} & {\dfrac{{\partial {{\dot x}_{12}}}}{{\partial {x_{12}}}}}\end{array}} \right]^T}\left[ {\begin{array}{*{20}{l}}{{\lambda _1}}\\ {{\lambda _2}}\\ \vdots \\{{\lambda _{11}}}\\{{\lambda _{12}}}\end{array}} \right]\end{align}

The specific representation of each component is as follows:

(29) \begin{align} {\begin{cases}\begin{array}{l}{{\dot \lambda }_1} = \left( {{x_2}{\lambda _3}\cos {x_3} + {\lambda _4}{x_2}\cos {x_3}\dfrac{{\cos {x_6}}}{{\cos {x_5}}} + {x_2}{\lambda _5}\cos {x_3}\sin {x_6}} \right.\\[10pt]\qquad \left. { - {x_2}{\lambda _6}\cos {x_3}\sin {x_5}\dfrac{{\cos {x_6}}}{{\cos {x_5}}}} \right)\dfrac{1}{{x_1^2}} - \left( {2{\mu _E}{\lambda _2}\sin {x_3} + 2{\mu _E}{\lambda _3}\dfrac{{\cos {x_3}}}{{{x_2}}}} \right)\dfrac{1}{{x_1^3}}\end{array}\\[15pt]\begin{array}{l}{{\dot \lambda }_2} = - {\lambda _3}\left[ {\cos {x_3}\left( {\dfrac{1}{{{x_1}}} + \dfrac{{{\mu _E}}}{{x_1^2x_2^2}}} \right) - \dfrac{{{T_P}}}{{{m_P}x_2^2}}\sin {u_1}\cos {u_2}} \right] - {\lambda _4}\cos {x_3}\dfrac{{\cos {x_6}}}{{{x_1}\cos {x_5}}}\\[10pt]\qquad - {\lambda _5}\cos {x_3}\dfrac{{\sin {x_6}}}{{{x_1}}} + {\lambda _6}\left( {\dfrac{{{T_P}\sin {u_2}}}{{{m_P}x_2^2\cos {x_3}}} + \dfrac{{\cos {x_3}\sin {x_5}\cos {x_6}}}{{{x_1}\cos {x_5}}}} \right) - {\lambda _1}\sin {x_3}\end{array}\\[15pt]\begin{array}{l}{{\dot \lambda }_3} = - {x_2}{\lambda _1}\cos {x_3} + {\mu _E}{\lambda _2}\dfrac{{\cos {x_3}}}{{x_1^2}} + {\lambda _3}\sin {x_3}\left( {\dfrac{{{x_2}}}{{{x_1}}} - \dfrac{{{\mu _E}}}{{x_1^2{x_2}}}} \right) + {x_2}{\lambda _4}\sin {x_3}\dfrac{{\cos {x_6}}}{{{x_1}\cos {x_5}}}\\[10pt]\qquad + {x_2}{\lambda _5}\sin {x_3}\dfrac{{\sin {x_6}}}{{{x_1}}} - {\lambda _6}\left[ {\dfrac{{{T_P}\sin {u_2}\sin {x_3}}}{{{m_P}{x_2}{{\left( {\cos {x_3}} \right)}^2}}} + {x_2}\dfrac{{\sin {x_3}\sin {x_5}\cos {x_6}}}{{{x_1}\cos {x_5}}}} \right]\end{array}\\[7pt] {{{\dot \lambda }_4} = 0}\\[3pt] {{{\dot \lambda }_5} = \dfrac{{{x_2}\cos {x_3}\cos {x_6}}}{{{x_1}{{\left( {\cos {x_5}} \right)}^2}}}\left( { - {\lambda _4}\sin {x_5} + {\lambda _6}} \right)}\\[10pt] {{{\dot \lambda }_6} = \dfrac{{{x_2}\cos {x_3}}}{{{x_1}\cos {x_5}}}\left( {{\lambda _4}\sin {x_6} - {\lambda _5}\cos {x_5}\cos {x_6} - {\lambda _6}\sin {x_5}\sin {x_6}} \right)}\end{cases}} \end{align}

Replace the subscript of Equation (29) to get the component expression of the target’s co-state equation.

According to the boundary conditions of the covariates in Equations (12) and (13), the results are as follows:

\begin{align*}pursuer\left\{ {\begin{array}{*{20}{l}}{{\lambda _1}\left( {{t_f}} \right) = {\upsilon _1}}\\[6pt]{{\lambda _4}\left( {{t_f}} \right) = {\lambda _4} = {\upsilon _2}}\\[6pt]{{\lambda _5}\left( {{t_f}} \right) = {\upsilon _3}}\\[6pt]{{\lambda _2}\left( {{t_f}} \right) = {\lambda _3}\left( {{t_f}} \right) = {\lambda _6}\left( {{t_f}} \right) = 0}\end{array}} \right.target\left\{ {\begin{array}{*{20}{l}}{{\lambda _7}\left( {{t_f}} \right) = - {\upsilon _1}}\\[6pt]{{\lambda _{10}}\left( {{t_f}} \right) = {\lambda _{10}} = - {\upsilon _2}}\\[6pt]{{\lambda _{11}}\left( {{t_f}} \right) = - {\upsilon _3}}\\[6pt]{{\lambda _8}\left( {{t_f}} \right) = {\lambda _9}\left( {{t_f}} \right) = {\lambda _{12}}\left( {{t_f}} \right) = 0}\end{array}} \right.\end{align*}

The conditions that the covariates satisfy are obtained:

(30) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{\lambda _1}\left( {{t_f}} \right) + {\lambda _7}\left( {{t_f}} \right) = 0}\\[6pt]{{\lambda _4}\left( {{t_f}} \right) + {\lambda _{10}}\left( {{t_f}} \right) = {\lambda _4} + {\lambda _{10}} = 0}\\[6pt]{{\lambda _5}\left( {{t_f}} \right) + {\lambda _{11}}\left( {{t_f}} \right) = 0}\\[6pt]{{\lambda _2}\left( {{t_f}} \right) = {\lambda _3}\left( {{t_f}} \right) = {\lambda _6}\left( {{t_f}} \right) = {\lambda _8}\left( {{t_f}} \right) = {\lambda _9}\left( {{t_f}} \right) = {\lambda _{12}}\left( {{t_f}} \right) = 0}\end{array}} \right.\end{align}

For the control variables, according to the first-order conditions of Equation (17), we can get the following results:

(31) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\dfrac{{\partial H}}{{\partial {u_1}}} = \dfrac{{{T_P}}}{{{m_P}{x_2}}}\cos {u_2}\left( {{\lambda _3}\cos {u_1} - {x_2}{\lambda _2}\sin {u_1}} \right) = 0}\\[8pt]{\dfrac{{\partial H}}{{\partial {u_2}}} = - \dfrac{{{T_P}}}{{{m_P}}}{\lambda _2}\sin {u_2}\cos {u_1} + \dfrac{{{T_P}}}{{{m_P}{x_2}}}\left( {{\lambda _6}\dfrac{{\cos {u_2}}}{{\cos {x_3}}} - {\lambda _3}\sin {u_1}\sin {u_2}} \right) = 0}\\[8pt]{\dfrac{{\partial H}}{{\partial {u_3}}} = \dfrac{{{T_E}}}{{{m_E}{x_8}}}\cos {u_4}\left( {{\lambda _9}\cos {u_3} - {x_8}{\lambda _8}\sin {u_3}} \right) = 0}\\[8pt]{\dfrac{{\partial H}}{{\partial {u_4}}} = - \dfrac{{{T_E}}}{{{m_E}}}{\lambda _8}\sin {u_4}\cos {u_3} + \dfrac{{{T_E}}}{{{m_E}{x_8}}}\left( {{\lambda _{12}}\dfrac{{\cos {u_4}}}{{\cos {x_9}}} - {\lambda _9}\sin {u_3}\sin {u_4}} \right) = 0}\end{array}} \right.\end{align}

Assume that $\cos {u_2} \ne 0$ , $\cos {u_4} \ne 0$ , then the possible values of ${u_1}$ , ${u_2}$ , ${u_3}$ and ${u_4}$ can be obtained as follows:

(32) \begin{align}\begin{array}{l}\begin{array}{*{20}{c}}{u_1^{\left( A \right)} = \arctan \left( {\dfrac{{{\lambda _3}}}{{{x_2}{\lambda _2}}}} \right)} & {u_1^{\left( B \right)} = \arctan \left( {\dfrac{{{\lambda _3}}}{{{x_2}{\lambda _2}}}} \right) + \pi }\end{array}\\[9pt]u_2^{\left( A \right)} = \arctan \left( {\dfrac{{{\lambda _6}}}{{\cos {x_3}\left( {{x_2}{\lambda _2}\cos {u_1} + {\lambda _3}\sin {u_1}} \right)}}} \right)\\[9pt]u_2^{\left( B \right)} = \arctan \left( {\dfrac{{{\lambda _6}}}{{\cos {x_3}\left( {{x_2}{\lambda _2}\cos {u_1} + {\lambda _3}\sin {u_1}} \right)}}} \right) + \pi \\[9pt]\begin{array}{*{20}{c}}{u_3^{\left( A \right)} = \arctan \left( {\dfrac{{{\lambda _9}}}{{{x_8}{\lambda _8}}}} \right)} & {u_3^{\left( B \right)} = \arctan \left( {\dfrac{{{\lambda _9}}}{{{x_8}{\lambda _8}}}} \right) + \pi }\end{array}\\[9pt]u_4^{\left( A \right)} = \arctan \left( {\dfrac{{{\lambda _{12}}}}{{\cos {x_9}\left( {{x_8}{\lambda _8}\cos {u_3} + {\lambda _9}\sin {u_3}} \right)}}} \right)\\[9pt]u_4^{\left( B \right)} = \arctan \left( {\dfrac{{{\lambda _{12}}}}{{\cos {x_9}\left( {{x_8}{\lambda _8}\cos {u_3} + {\lambda _9}\sin {u_3}} \right)}}} \right) + \pi \end{array}\end{align}

Therefore, for each participant, there are four possible control modes satisfying the first-order condition.

There are four control modes for the pursuer and the target spacecraft respectively:

(33) \begin{align}\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{pursuer} & {\begin{array}{*{20}{c}}{\left\{ {u_1^{\left( A \right)},u_2^{\left( A \right)}} \right\}} & {\left\{ {u_1^{\left( A \right)},u_2^{\left( B \right)}} \right\}}\\[10pt]{\left\{ {u_1^{\left( B \right)},u_2^{\left( A \right)}} \right\}} & {\left\{ {u_1^{\left( B \right)},u_2^{\left( B \right)}} \right\}}\end{array}}\end{array}} & {\begin{array}{*{20}{c}}{target} & {\begin{array}{*{20}{c}}{\left\{ {u_3^{\left( A \right)},u_4^{\left( A \right)}} \right\}} & {\left\{ {u_3^{\left( A \right)},u_4^{\left( B \right)}} \right\}}\\[10pt]{\left\{ {u_3^{\left( B \right)},u_4^{\left( A \right)}} \right\}} & {\left\{ {u_3^{\left( B \right)},u_4^{\left( B \right)}} \right\}}\end{array}}\end{array}}\end{array}\end{align}

The correct and appropriate control strategy can be selected by considering the second-order condition in Equation (18), that is, the following relationship exists:

\begin{align*}\begin{array}{l}\dfrac{{{\partial ^2}H}}{{\partial \boldsymbol{u}_P^2}} = \left[ {\begin{array}{ccc}{H_{11}^{\left( P \right)}} & & {H_{12}^{\left( P \right)}}\\[3pt]{H_{21}^{\left( P \right)}} && {H_{22}^{\left( P \right)}}\end{array}} \right] \geq 0\\[15pt]\dfrac{{{\partial ^2}H}}{{\partial \boldsymbol{u}_E^2}} = \left[ {\begin{array}{ccc}{H_{11}^{\left( E \right)}} && {H_{12}^{\left( E \right)}}\\[3pt]{H_{21}^{\left( E \right)}} && {H_{22}^{\left( E \right)}}\end{array}} \right] \le 0\end{array}\end{align*}

The sub matrices in the above equation are as follows:

\begin{align*}\begin{array}{*{20}{c}}{H_{11}^{\left( P \right)} = - \dfrac{{{T_P}\cos {u_2}}}{{{m_P}{x_2}}}\left( {{x_2}{\lambda _2}\cos {u_1} + {\lambda _3}\sin {u_1}} \right)}\\[9pt]{H_{12}^{\left( P \right)} = H_{21}^{\left( P \right)} = \dfrac{{{T_P}\sin {u_2}}}{{{m_P}{x_2}}}\left( {{x_2}{\lambda _2}\sin {u_1} - {\lambda _3}\cos {u_1}} \right)}\\[9pt]\begin{array}{l}H_{22}^{\left( P \right)} = - \dfrac{{{T_P}}}{{{m_P}}}\left[ {\dfrac{1}{{{x_2}}}\left( {{\lambda _3}\sin {u_1}\cos {u_2} + {\lambda _6}\dfrac{{\sin {u_2}}}{{\cos {x_3}}}} \right) + {\lambda _2}\cos {u_1}\cos {u_2}} \right]\\[9pt]\begin{array}{*{20}{c}}{H_{11}^{\left( E \right)} = - \dfrac{{{T_E}\cos {u_4}}}{{{m_E}{x_8}}}\left( {{x_8}{\lambda _8}\cos {u_3} + {\lambda _9}\sin {u_3}} \right)}\\[9pt]{H_{12}^{\left( E \right)} = H_{21}^{\left( E \right)} = \dfrac{{{T_E}\sin {u_4}}}{{{m_E}{x_8}}}\left( {{x_8}{\lambda _8}\sin {u_3} - {\lambda _9}\cos {u_3}} \right)}\\[9pt]{H_{22}^{\left( E \right)} = - \dfrac{{{T_E}}}{{{m_E}}}\left[ {\dfrac{1}{{{x_8}}}\left( {{\lambda _9}\sin {u_3}\cos {u_4} + {\lambda _{12}}\dfrac{{\sin {u_4}}}{{\cos {x_9}}}} \right) + {\lambda _8}\cos {u_3}\cos {u_4}} \right]}\end{array}\end{array}\end{array}\end{align*}

The above two second-order conditions are also equivalent to the following two inequalities:

\begin{align*}\begin{array}{*{20}{c}}{H_{11}^{\left( P \right)} + H_{22}^{\left( P \right)} \geq 0} & {H_{11}^{\left( P \right)}H_{22}^{\left( P \right)} - H_{12}^{\left( P \right)}H_{21}^{\left( P \right)} \geq 0}\\[6pt]{H_{11}^{\left( E \right)} + H_{22}^{\left( E \right)} \le 0} & {H_{11}^{\left( E \right)}H_{22}^{\left( E \right)} - H_{12}^{\left( E \right)}H_{21}^{\left( E \right)} \geq 0}\end{array}\end{align*}

The appropriate control rate at any time can be selected from the four control schemes in Equation (33) by the above two inequalities.

When the terminal interception time ${t_f}$ is free, the terminal condition for the Hamiltonian function of the system is:

(34) \begin{align}\frac{{\partial \Phi }}{{\partial {t_f}}} + H\left( {{t_f}} \right) = \frac{{\partial \left( {\varphi + {{\boldsymbol{\upsilon }}^T}{\boldsymbol{\psi }}} \right)}}{{\partial {t_f}}} + H\left( {{t_f}} \right) = \frac{{\partial {t_f}}}{{\partial {t_f}}} + H\left( {{t_f}} \right) = 1 + {H_f} = 0\end{align}

Since a is time independent, its value is constant, i.e. $\forall t,H = - 1$ .

For bilateral optimal control, according to the above derivation, the control variables represented by covariates and state variables are obtained, and the optimal control in continuous space is transformed into Two Point Boundary Value Problem, that is to solve a series of nonlinear differential equations.

5.0 Hybrid optimisation method for optimal interception orbit

Because the dynamic equation of the two-point boundary value problem is a nonlinear equation, it is difficult to obtain its analytical solution. Generally, numerical solutions are used to solve the problem, such as the near extreme value method, difference method or shooting method. The exact solution must have a good initial value to get the global optimal solution, otherwise it may fall into the local optimal region. In order to solve this problem, this paper will use the hybrid optimisation algorithm of global search to get the suboptimal solution and local optimisation to obtain the optimal value. That is, firstly, the sub optimal solution of the optimisation problem is obtained by the global search of genetic algorithm, and then it is regarded as the sequence quadratic program. The initial value of SQP algorithm can overcome the disadvantage that SQP algorithm is sensitive to the initial parameters, it can meet the requirements of convergence and accuracy as well. The flow chart of the proposed hybrid optimisation method is given as follows in Fig. 9.

Figure 9. Flow chart of hybrid optimisation method.

The design variable is defined as $\boldsymbol{x} = \left[ {{\lambda _{1,0}},{\lambda _{2,0}},{\lambda _{3,0}},{\lambda _4},{\lambda _{5,0}},{\lambda _{6,0}},{\lambda _{7,0}},{\lambda _{8,0,}}{\lambda _{9,0}},{\lambda _{10}},{\lambda _{11,0}},{\lambda _{12,0}},{t_f}} \right]$ , where ${{\boldsymbol{\lambda }}_0}$ is the initial value of the covariate variables and ${t_f}$ is the end time of the game.

In order to transform the two-point boundary value problem into an optimisation problem without constraints, the objective function of design optimisation is the square of the distance from the set of functions composed of Equations (26), (30) and (34), is:

(35) \begin{align}J = \sum\limits_{i = 1}^{13} {{k_i}{\eta _i}c_i^2} \end{align}

Where ${k_i}$ is the weight, ${\eta _i} = 0/1$ , when ${\eta _i} = 0$ , the distance is less than $\varepsilon $ , and $\varepsilon $ is a specified small value. In the process of solving, the equation condition is gradually satisfied with the distance approaching 0.

After getting the initial values of the co-state variables and interception time by genetic algorithm, the SQP algorithm is used for accurate optimisation. In order to solve the above problems, in addition to the design variables and objective functions unchanged, the orbital heights of the pursuer and the target spacecraft should be added as inequality constraints, that is, the orbital altitudes ${r_P}$ and ${r_E}$ of the two spacecraft at any time should be greater than a certain safe altitude:

\begin{align*}\begin{array}{l}\left( {{R_e} + {h_{safe}}} \right) - {r_P} \le 0\\\left( {{R_e} + {h_{safe}}} \right) - {r_E} \le 0\end{array}\end{align*}

Where ${R_e}$ is the earth radius, ${h_{safe}}$ is the safe altitude. The above two inequalities require two spacecraft will not enter the earth’s atmosphere during flight. Moreover, it must be satisfied at any time in the integration process.

6.0 Simulation results and analysis

In order to facilitate the calculation, normalisation is carried out. Set unit distance ${\rm{DU}} = 1{R_E} = 6378.165\left( {{\rm{km}}} \right)$ , unit time ${\rm{TU}} = \sqrt {R_E^3/\mu } = 806.8\left( {\rm{s}} \right)$ , unit acceleration ${\rm{AU = DU/T}}{{\rm{U}}^2} = g = 9.798 \times {10^{ - 3}}\left( {{\rm{km/}}{{\rm{s}}^2}} \right)$ , gravitational constant $\mu = 1$ . Set ${T_P}/{m_P} = 0.1{\rm{AU}}$ , ${T_E}/{m_E} = 0.05{\rm{AU}}$ , ${h_{safe}} = 100$ km, the orbital parameters of the pursuer and the target satellite are shown in Table 1.

Table 1. Six elements of the track of the pursuer and the target

It can be seen from Table 1 that the pursuer operates on an elliptical orbit with an eccentricity of 0.044, and the target runs on a circular orbit with an orbit height of 1,000 km. The corresponding state variables can be obtained through the transformation relationship in Section 2, as shown in Table 2.

Table 2. Initial values of state variables for pursuer and target

Firstly, 12 covariates are estimated initially, and the upper and lower limits of the value range are given. Then, the genetic algorithm is used to solve the initial problem, where the population size is set to 50, the mating probability and mutation probability are the default values of the algorithm, which are 0.8 and 0.2, respectively. Then the optimal fitness function and average fitness function of each generation in the global optimisation process can be obtained, and the sub-optimal solutions of covariate and time are shown in Table 3.

Table 3. Covariate variables and sub optimal solution of interception end time

Moreover, the average convergence performances of two algorithms (GA-SQP based on differential games: GA-SQP-DG, UOC: unilateral optimal control method) are given in this section. ${J_{initial\_g}}$ , ${J_{generation(i)}}$ are assumed as the average cost in the initial generation and in $i$ th generation, and the ratio indicated as ${J_{generation(i)}}/{J_{initial\_g}}$ is adopted to display the result.

The comparison result is shown in Fig. 10.

Figure 10. Convergence rate comparison.

Figure 11. Intercepting trajectory of pursuer and target.

As shown in Fig. 10, the vertical axis designates the average cost value of 50 runs at each generation normalised by the average cost of the initial generation. The horizontal axis designates the generation number. It is noticeable that the convergence rate of GA-SQP-DG is faster than that of UOC.

Taking the results obtained by genetic algorithm in Table 3 as the initial value of SQP algorithm, the accurate solutions of covariate variables and intercept end time are obtained by optimisation.

Figures 11 and 12 show the transfer orbit and the relative distance change between the two spacecraft under the optimal strategy. As can be seen from the figures, the pursuer can successfully intercept the target in a limited time, and the intercepting orbit is in the same direction as the target’s orbit.

Figure 12. The relative distance between the pursuer and the target.

Figures 13–15 show the time-dependent thrust angles of the pursuer and the target satellite during interception.

Figure 13. Angle-of-attack of pursuer.

Figure 14. Angle of sideslip of pursuer.

Figure 15. Thrust angle of target satellite.

Figures 15–21 show the change of the state variables of the pursuer and the target satellite with time during the interception process. As can be seen from the figures, only the final position is required to be the same at the end of interception, that is, $\left( {{r_i},{\xi _i},{\phi _i}} \right)$ is consistent. The size of $\left( {{\gamma _i},{v_i},{\zeta _i}} \right)$ is not required.

Figure 16. Geocentric distance of pursuer and target.

Figure 17. Longitude of pursuer and target.

Figure 18. Latitude of pursuer and target.

Figure 19. Flight path angle of pursuer and target.

Figure 20. Velocity of pursuer and target.

Figure 21. Azimuth angle of pursuer and target.