A hybrid optimisation method for intercepting satellite trajectory based on differential game

Abstract

This study addresses orbit design and optimisation for the situation of satellite interception in which the target spacecraft is capable of manoeuvring using continuous magnitude restricted thrust. For the purpose of designing a long-range continuous thrust interception orbit, the orbit motion equations of two satellites with J2 perturbation are constructed. This problem is assumed to be a typical pursuit-evasion problem in differential game theory; using boundary constraint conditions and a performance index function that includes time and fuel consumption, the saddle point solution corresponding to the bilateral optimal is derived, and then this pursuit-evasion problem is transformed into a two-point boundary value problem. A hybrid optimisation method using a genetic algorithm (GA) and sequential quadratic programming (SQP) is derived to obtain the optimal control strategy. The proposed model and algorithm are proved to be feasible for the given simulation cases.

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)

1.0 Introduction

Interception of satellite signals has become a contentious topic in aeronautical engineering. Satellite interception can be split into direct ascent and common orbit types based on the mechanism of launch (Fig. 1). In terms of energy consumption, the two methods are vastly different: The former requires a launch vehicle to launch an intercepting weapon near the target satellite, which results in a high level of energy consumption; the latter, on the other hand, uses satellites, which are light in weight, inexpensive to manufacture, and have high agility and strong autonomous mobility. Satellite interception technology relies on efficient interception; thus, the orbit design and optimisation method of intercepting satellites is becoming an essential issue in order to overcome the bottleneck.

Figure 1. Two typical intercepting satellite modes.

Control systems for intercepting satellites with little mobility have been extensively examined in a large number of studies. It was found that orbital planar motion and motion outside of the orbital plane can be used to decompose satellite relative motion using the well-known Clohessy-Wiltshire (CW) equation. With the fuel expense in consideration, they devised a relative orbital control technique that included three in-plane and one out-of-plane impulsive manoeuvres. An impulsive approach to a near-orbit rendezvous or interception is simple to conceive or implement. However, as impulsive guidance is an open-loop control system, its precision generally falls short of the mission’s specifications. Another more probable real situation is that the interceptor and target satellites can both perform orbital manoeuvres with thrusters onboard, hence it is commonly believed that both spacecraft utilise continuous thrust manoeuvres. Therefore, this topic revolves around finding the best trajectory for a continuous thrust interception with the target’s manoeuvring.

As shown in Fig. 2 is the overview of literatures on satellite or spacecraft interception orbital optimisation. According to the literature, the three most significant categories for this orbital trajectory optimisation are Nonlinear programming (NLP), Metaheuristics and Hybrid algorithms.

Figure 2. Overview of literatures on satellite/spacecraft interception orbital optimisation.

When it comes to the optimisation of space trajectories, gradient-based approaches like NLP appears to be the most prevalent ways for solving problems about trajectory optimisation. As a result, NLP is frequently capable of achieving rapid convergence and obtaining extremely precise outcomes. Because of this, there have been a lot of individual software programs since 1970 (SNOPT, [1] rSQP++, [2] KNITRO [3]). Gradient-based approaches, on the other hand, require an initial guess at all the system’s parameters. Analysers must be familiar with the ideal trajectory because all the nodal states and control values are based on parameters. Poor initial guessing frequently results in the cost function failing to converge a non-global optimum solution. The gradient-based method’s first guess is a difficult problem to solve in spacecraft trajectory optimisation. This issue was originally addressed in 1995, [4] which led to the introduction of metaheuristics, an alternate methodology to address the trajectory optimisation problems.

Significant progress has been made in the application of getting approximate solutions to spacecraft trajectory optimisation problems via metaheuristics in recent years. When it comes to solving spacecraft trajectory optimisation challenges, many publications have employed evolutionary algorithms, which are the most common instances. Such as GA, Particle Swarm Optimisation (PSO) and Genetic Algorithm-Particle Swarm Optimisation (GA-PSO), etc. are widely used. Shafieenejad [5] proposes an approach for optimising continuous-thrust orbital transfer trajectories. Abdelkhalik [6] puts forward a new approach for calculating the ideal orbital transfer trajectory based on a genetic algorithm. In Ref. (7) PSO is adopted in solving difficulties like finding periodic orbits within the confines of the circular restricted three-body problem and optimising orbital transfers. Imperialistic Competition Algorithm is utilised to solve low-thrust orbital transfer problems, which included nonlinear dynamic equations in Ref. (8). Detailed models, aims, methodologies, and solutions for spacecraft trajectory optimisation are reviewed by Shirazi [9]. Even though evolutionary algorithms appear to be the primary choice for most spacecraft trajectory optimisation problems because of their availability and ease of implementation. Metaheuristics are more commonly applied in issues with impulsive models rather than continuous models, according to Ref. (10).

Hybrid optimisation strategies have become increasingly popular in recent years. Techniques like this can often yield the best outcomes, especially in real-world scenarios. The use of both metaheuristics and other techniques in concert is becoming more commonplace. GA and SQP are employed to solve the general multiple-impulse rendezvous problem in Ref. (11). Vasile [12] offers a novel multi-objective agent-based approach that combines a number of heuristics to address three distinct issues in space trajectory design. A hybrid metaheuristic based on GA and simulated annealing (SA) is employed and outperforms traditional GA and SA in the optimisation of orbital manoeuvres for satellites in Ref. (13). Sentinella [14] designs and uses a hybrid evolutionary method for spaceship trajectory optimisation that synergistically uses differential evolution, genetic algorithms and particle swarm optimisation. A new method based on differential games theory has recently been devised, which naturally treats the interception problem as a game of pursuit-evasion in which both sides have movement capabilities. To solve this challenge in particular, Isaacs [15] first focuses on it and comes up with the saddle solution. A modified first-order differential dynamic programming technique was used to construct a near-optimal feedback control for minimax-range pursuit-evasion problems between two constant-thrust spacecraft in Ref. (16). This paper, on the other hand, assumes that the satellites have a high degree of maneuverability, which is simply not true. A numerical technique for solving the open-loop trajectory of a three-dimensional satellite pursuit-evasion interception is presented in Refs (17) and (18), in which each spacecraft has a limited ability to move. The pursuer’s goal in the interception was to hit the target satellite as quickly as feasible, while the evader’s goal was to delay that moment as long as possible. Genetic algorithms were used to initially generate a pre-solution to the saddle-point equilibrium problem. The initial guess was then used in the semi-analytic procedure to determine the precise trajectory of the pursuit-evasion, and a global solution for the nonlinear pursuit-evasion games can be found using an intensive random search and collocation method, but it uses a lot of computing resources. Closed-loop control is more complex to derive than open-loop control. To synthesise almost optimum feedback controllers for nonlinear two player pursuit-evasion games, Ghosh [19] uses extremal fields. A surrogate model of a feedback controller was constructed using universal Kriging to generate sub-optimal control based on current states. The real-time feedback map was constructed by interpolating controls of these open-loop extremals. Stupik [20] researches satellite interception using the linearised CW equation. The typical solutions interpolate the sub-optimal feedback solution with various initial conditions. The number of conjugates to solve drops from 12 to 3. Despite the method is for the open solution, the open loop extremal’s offline solving ability accordingly improves. Although the authors used Kriging to build a real-time feedback control, they cannot guarantee optimality intercept trajectory. Gutman [21] proposes a nonlinear vector guidance law for exo-atmospheric interception using steering jets as the only viable way. Time-to-go is presented as the solution of a quartic polynomial equation for both ideal and non-ideal interceptors, and the proposed optimal guidance law could optimise both parties’ capture time. Zhou [22] designs the pursuer’s control strategy based on fuzzy comprehensive evaluation and differential game. This is done by introducing relative state variables and zero effort miss (ZEM) variables into the Lawden equations, and the orbital pursuit-evasion-defense problem is solved using a hybrid method integrating multi-objective genetic algorithms, but the satellites are assumed to be highly maneuverable. For two separate spacecraft tracking instances, Wang [23] develops two reinforcement learning parameter-self-tuning controllers: tracking a disabled target and a mobile target, and the reinforcement learning is verified to be effective for achieving desired purposes, but one of the flaws of learning-based algorithms is their reliance on large samples of data.

However, most aforementioned works assume that either the interceptor always has adequate fuel or maneuverability during the game or that the target follows a single strategy over the course of the game’s action. In a real interception procedure, as the target is unwilling to cooperate, the interceptor is unable to accurately detect the target’s information, and the maneuverability of interceptor should be strictly regulated. To fill the gap between applications and methodologies, an investigation of the interception of a target satellite with a continuous magnitude restricted thrust is conducted in this study.

The contributions of this paper can be summarised as follows:

1. (1) The J2 perturbation is considered in the orbit motion equation of the continuous thrust interception orbit design, and the mass change of the spacecraft is included. The performance index function of the zero-sum differential game problem considers that the interception time and fuel consumption are both optimal.

2. (2) Based on differential game theory, a hybrid optimisation algorithm is designed with full consideration of target maneuverability, which can effectively improve the efficiency of problem solving.

The rest of the paper is organised as follows. In the forthcoming section, the orbit motion equation of two satellites with J2 term perturbation is established by defining state variables. Next, based on the differential game theory, the corresponding saddle point solution conditions are derived according to the boundary constraint conditions and the performance index function including time and fuel consumption. Then we introduce the zero miss variables and propose a hybrid algorithm combines GA with SQP to solve the problem, and the optimal control strategy is obtained consequently. Several numerical simulations are used to analyse the proposed algorithm. Figure 3 describes the main content and the main solution process of this paper.

Figure 3. Main content and the main solution process of this paper.

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: [24]

(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.

7.0 Conclusions

A detailed study of intercepting satellite trajectory optimisation method based on differential game is achieved with J2 perturbation.

We established a hybrid method which combines SQP with GA to manage the given pursue-evasion problem. The corresponding saddle point solutions for the bilateral optimisation are derived based on differential game theory, then the sub-optimal solution is obtained by GA, set both interception time and fuel consumption as the performance index, then SQP is adopted to transform the given solution into the change of thrust angle of the two spacecraft in the pursuit process, which is the optimal control strategy of both sides.

Based on simulation studies, the effectiveness and practicality the proposed hybrid algorithm are both verified.

In our future research, the influence of other disturbance items to the relative kinematics model will be considered.