Pursuer’s Control Strategy for Orbital Pursuit-Evasion-Defense Game with Continuous Low Thrust Propulsion

Abstract

This paper studies the orbital pursuit-evasion-defense problem with the continuous low thrust propulsion. A control strategy for the pursuer is proposed based on the fuzzy comprehensive evaluation and the differential game. First, the system is described by the Lawden’s equations, and simplified by introducing the relative state variables and the zero effort miss (ZEM) variables. Then, the objective function of the pursuer is designed based on the fuzzy comprehensive evaluation, and the analytical necessary conditions for the optimal control strategy are presented. Finally, a hybrid method combining the multi-objective genetic algorithm and the multiple shooting method is proposed to obtain the solution of the orbital pursuit-evasion-defense problem. The simulation results show that the proposed control strategy can handle the orbital pursuit-evasion-defense problem effectively.

1. Introduction

Recently, the orbital pursuit-evasion problem has attracted increasing attention in space research [1,2,3,4]. This problem can be formulated as a differential game [5], which aims to obtain the optimal control strategy of the pursuer and/or the evader in the worst-case scenario, so as to realize the interception of the evader or the evasion from the pursuer.

Wong [6] was regarded as the first person to study the orbital pursuit-evasion problem, he solved the problem of intercepting a maneuverable satellite under the assumption of planar motion and constant gravitational field. Since then, many works have focused on the orbital pursuit-evasion problem. In reference [7], a method based on periodically updating the solution of the two-point boundary value problem (TPBVP) was proposed to generate near optimal feedback controls for the orbital pursuit-evasion problem. However, this method is time-consuming and difficult to be applied in real time. In order to overcome these drawbacks, Anderson [8] used a modified first-order differential dynamic programming algorithm to generate near-optimal feedback controls. References [9,10,11] found the saddle-point equilibrium solutions of the three-dimensional orbital pursuit-evasion game respectively by three different hybrid numerical methods. Hafer et al. [12] applied the sensitivity method to the orbital pursuit-evasion problem, which greatly reduces the computation burden for solving this problem numerically. Widhalm studied the problem of avoiding an interception and proposed two optimal evasive-maneuver strategies with the impulsive thrust [13] and the continuous low thrust [1] respectively. Prussing et al. [14] derived minimum-fuel impulsive strategies for return-on-state maneuvers by applying the primer vector theory. Merz [15] developed the guidance laws for the noisy satellite pursuit-evasion game. Woodbury et al. [16] studied an incomplete, imperfect information game and presented the adaptive strategies for the pursuer and the evader. Ghosh et al. [17] developed a near-optimal feedback controller for the two-player pursuit-evasion games by using a new extremal-field approach. The above works were studied in the two-player pursuit-evasion game framework. However, in this framework, the evader can only perform maneuvers by itself to avoid threats. It is called self-defense, which disturbs the original mission of the evader and requires a large additional amount of fuel.

To overcome this disadvantage, a defender is introduced in [18]. The role of the defender is intercepting the pursuer. In this way, the evader can perform its original mission without being disturbed. A hybrid method combined particle swarm optimization with a Newton-Interpolation algorithm was proposed to solve the orbital defense problem. However, because of the introduction of the defender, the pursuer must avoid the interception by the defender while capturing the evader [19], which makes the design of the pursuer’s control strategy more complicated. In order to develop control strategies for pursuers, Liu et al. [19] proposed a distributed online mission plan algorithm for pursuers to access targets. However, these works on the orbital pursuit-evasion-defense game adopted the impulsive thrust, which suffers the drawback that the interception will fail when the target can perform evasive maneuvers [4].

Compared with the impulse thrust, the continuous low thrust allows players to perform multiple, continuous maneuvers, which meets the requirements of the frequently orbital transfers in the game. When applying the continuous low thrust, the hypothesis about players’ maneuverable is removed. It is closer to the actual situation of the orbital pursuit-evasion-defense game. Therefore, in this paper, the orbital pursuit-evasion-defense game model is constructed based on the continuous low thrust. Different from the model based on impulse thrust, the model based on continuous low thrust cannot adopt the Keplerian dynamics [20]. Its dynamic equations are based on the non-Keplerian motion. Two issues need to be solved in this model: (i) The system has a high dimension, which means that it will suffer from the curse of dimensionality [21] when solving the problem; (ii) two objectives, intercepting the evader and evading the defender, should be considered by the pursuer, and the corresponding weights should be determined according to the current state. For the first issue, as the zero effort miss (ZEM) can be used to simplify the linear system [4], the dimension of the system is reduced by introducing the relative state variables and the ZEM variables [22]. For the second issue, the pursuer’s objective function is designed based on the fuzzy comprehensive evaluation, and the pursuer’s control strategy which is suitable for the orbital pursuit-evasion-defense game is proposed. Based on the above model, the orbital pursuit-evasion-defense game is transformed into a TPBVP by applying the differential game theory. A hybrid method combining the multi-objective genetic algorithm and the multiple shooting method is presented to solve the TPBVP.

2. Mathematical Model of Orbital Pursuit-Evasion-Defense Game

2.1. Relative Orbital Dynamics

The orbital pursuit-evasion-defense game occurs in the final phase of the confrontation when the spacecraft are close enough so that they can identify each other with onboard electronic devices [4]. In this type of situation, the motion between the spacecraft can be expressed as relative motion [23]. As is known, Lawden’s equations [24] and Clohessy–Wiltshire (C–W) equations [25] are two linearized equations used to describe the relative motion between spacecrafts. Unlike the C–W equations, which can only be applied to circular orbits, the Lawden’s equations can describe the relative motion of a spacecraft in elliptical orbits. Same as in [26], the dynamics of the participating spacecraft are described in the local-vertical local-horizontal (LVLH) frame centered at a virtual spacecraft. In addition, Lawden’s equations are adopted as the relative dynamic equations of the spacecraft.

As shown in Figure 1, $P,D,E$ respectively represents the pursuer, the defender, and the evader. We establish an elliptical fictitious spacecraft $O$ which is close to the players. The LVLH coordinate system is centered at the point $O$. $OX$ is pointing outward along the radius of the Earth, $OY$ is perpendicular to $OX$ in the reference orbital plane and points to the front of its flight direction, $OZ$ is perpendicular to the orbital plane and forms a right-handed frame with $OX$ and $OY$.

The Lawden’s equations can be expressed as:

$$\{\begin{array}{l}{\ddot{x}}_i={\omega }^2{x}_i+2\omega {\dot{y}}_i+\dot{\omega }{y}_i+2\frac{\mu x_i}{r_t^3}+T_i{u}_{xi}\\ {\ddot{y}}_i=-2\omega {\dot{x}}_i-\dot{\omega }{x}_i+{\omega }^2{y}_i-\frac{\mu y_i}{r_t^3}+T_i{u}_{yi}\\ {\ddot{z}}_i=-\frac{\mu z_i}{r_t^3}+T_i{u}_{zi}\end{array},i=P,D,E$$

(1)

where $\mu$ is the Earth gravitational constant, $r_t$ is the distance between the origin $O$ and the Earth core, $\omega$ and $\dot{\omega }$ represent the orbital angular velocity and acceleration of the origin $O$, respectively. $x_i$, $y_i$, and $z_i$ represent the position components of the players in the relative coordinate system. $T_i$ represents the maximum thrust. $u_{xi}$, $u_{yi}$ and $u_{zi}$ respectively represent control variables in three directions (i.e., x, y, z axis), ranging from 0 to 1.

The state variables (i.e., position and velocity) of the players are represented by ${\mathit{X}}_i$ as follows:

$${\mathit{X}}_i={[x_i,y_i,z_i,{\dot{x}}_i,{\dot{y}}_i,{\dot{z}}_i]}^{\mathrm{T}},i=P,D,E$$

(2)

Thus, the dynamics equations can be written as:

$${\dot{\mathit{X}}}_i={\mathit{A}}_i{\mathit{X}}_i+T_i{\mathit{U}}_i$$

(3)

where

$${\mathit{A}}_i\left(t\right)=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ {\omega }^2+\frac{2\mu }{r_t^3}& \dot{\omega }& 0& 0& 2\omega & 0\\ -\dot{\omega }& {\omega }^2-\frac{\mu }{r_t^3}& 0& -2\omega & 0& 0\\ 0& 0& -\frac{\mu }{r_t^3}& 0& 0& 0\end{array}\right],i=P,D,E$$

(4)

${\mathit{U}}_i(i=P,D,E)$ is the control variable, which can be represented by

$${\mathit{U}}_i={[0,0,0,u_{xi},u_{yi},u_{zi}]}^{\mathrm{T}},\Vert {\mathit{U}}_i\Vert \le 1$$

(5)

2.2. Dimension-Reduction

According to Equation (3), each player has 6 state variables, so the total number of state variables in the game is 18. In the numerically solving process, co-state variables associated with state variables are introduced, and the problem converts to a 36-dimensional TPBVP. However, solving this high-dimensional TPBVP is quite difficult and possesses high computational demands [27]. In order to improve computational efficiency, the dimension of the system needs to be reduced. This process is performed in two steps. First, the relative state variables between the spacecraft are used to replace the system states. Then, the ZEM variables are applied to further reduce the number of variables and equations in the system.

In the first step, the game can be divided into two parts: One is the game between the pursuer $P$ and the evader $E$; the other is the game between the defender $D$ and the pursuer $P$. The relative state variables in two parts, ${\mathit{X}}_{\mathrm{PE}}$ and ${\mathit{X}}_{\mathrm{DP}}$, can be respectively represented by

$$\{\begin{array}{l}{\mathit{X}}_{\mathrm{PE}}={\mathit{X}}_{\mathrm{P}}-{\mathit{X}}_{\mathrm{E}}\\ {\mathit{X}}_{\mathrm{DP}}={\mathit{X}}_{\mathrm{D}}-{\mathit{X}}_{\mathrm{P}}\end{array}$$

(6)

Substituting Equation (6) into Equation (3), the state equations are converted to:

$$\{\begin{array}{l}{\dot{\mathit{X}}}_{\mathrm{PE}}=\mathit{A}{\mathit{X}}_{\mathrm{PE}}+T_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}-T_{\mathrm{E}}{\mathit{U}}_{\mathrm{E}}\\ {\dot{\mathit{X}}}_{\mathrm{DP}}=\mathit{A}{\mathit{X}}_{\mathrm{DP}}+T_{\mathrm{D}}{\mathit{U}}_{\mathrm{D}}-T_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}\end{array}$$

(7)

where $\mathit{A}={\mathit{A}}_{\mathrm{P}}={\mathit{A}}_{\mathrm{E}}={\mathit{A}}_{\mathrm{D}}$.

In the second step, according to the linear system theory, the zero-input state transfer matrix $\mathit{\Phi }(t_{\mathrm{f}},t)$ of the state equation is defined as:

$$\{\begin{array}{l}\dot{\mathit{\Phi }}(t_{\mathrm{f}},t)=-{\Phi }(t_{\mathrm{f}},t)\mathrm{A}\\ \mathit{\Phi }(t_{\mathrm{f}},t_{\mathrm{f}})={\mathit{I}}_6\end{array}$$

(8)

where $t_{\mathrm{f}}$ is the terminal time and ${\mathit{I}}_6$ is the $6\times 6$ unit matrix.

Although two factors, the relative position and the relative velocity are involved in the game, only the first factor needs to be considered at the end of the game. The ZEM is the miss distance if both players do not apply any control from the current moment to the end of the game. Thus, the ZEM variables are introduced to reduce the dimension of the system and defined as:

$$\{\begin{array}{l}{\mathit{Z}}_{\mathrm{PE}}(t)=\mathit{D}\mathit{\Phi }(t_{\mathrm{f}},t){\mathit{X}}_{\mathrm{PE}}\\ {\mathit{Z}}_{\mathrm{DP}}(t)=\mathit{D}\mathit{\Phi }(t_{\mathrm{f}},t){\mathit{X}}_{\mathrm{DP}}\end{array}$$

(9)

where $\mathit{D}=[I_{3\times 3},0_{3\times 3}]$.

Substituting Equation (9) into Equation (7), the state equations are reduced to:

$$\{\begin{array}{l}{\dot{\mathit{Z}}}_{\mathrm{PE}}=\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{E}}{\mathit{U}}_{\mathrm{E}}\\ {\dot{\mathit{Z}}}_{\mathrm{DP}}=\mathit{D}\mathit{\Phi }{T}_{\mathrm{D}}{\mathit{U}}_{\mathrm{D}}-\mathit{D}\mathsf{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}\end{array}$$

(10)

2.3. Design of Objective Function Based on Fuzzy Comprehensive Evaluation

In the orbital pursuit-evasion-defense game, the pursuer must survive from the defender’s interception before it can successfully access the evader. Therefore, the pursuer-evader game and the defender-pursuer game must be considered and weighed in the objective function of the pursuer. As noted by Liu [19], the fuzzy comprehensive evaluation is an effective way to quantify various factors that are difficult to evaluate. Thus, it is used to obtain the weights corresponding to the two games. Detailed design is shown as follows.

By taking the terminal miss distance as the cost, the objective function of the three players can be defined as:

$$\{\begin{array}{l}{J}_{\mathrm{E}}=-\Vert {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})\Vert \\ J_{\mathrm{D}}=\Vert {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})\Vert \\ J_{\mathrm{P}}=k_1\Vert {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})\Vert -k_2\Vert {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})\Vert \end{array}$$

(11)

where the parameter $k_i,i=1,2$ is the weight factor, which satisfies $k_i\ge 0$. $k_1>k_2$ indicates that the pursuer prefers to reduce the terminal miss distance of the pursuer-evader game, while $k_1<k_2$ indicates that the pursuer prefers to increase the terminal miss distance of the defender-pursuer game. The value of $k_i$ is divided into 11 scales, which are shown in

Table 1. The evaluation scales.

$v_i$1	1	2	3	4	5	6	7	8	9	10	11
$k_1$	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
$k_2$	1	0.9	0.8	0.7	0.6	0.5	0.4	0.3	0.2	0.1	0

¹$v_i,i=1,\cdots ,11$ represents the corresponding scales, respectively.

.

According to the analysis above, two factors need to be evaluated, one is the urgency of intercepting the evader at the very moment, denoted by $u_1$; the other is the urgency of evading the defender at the very moment, denoted by $u_2$. The effect of the factor $u_1$ increases as the ZEM distance of the pursuer-evader game decreases. The effect of the factor $u_2$ increases as the ZEM distance of the defender-pursuer game decreases. According to this rule, ${\mathit{Z}}_{\mathrm{PE}}(t)$ and ${\mathit{Z}}_{DP}(t)$ are used to construct the weights of the two factors, which are given by:

$$\{\begin{array}{l}{a}_1=1-{\left(\frac{\Vert {\mathit{Z}}_{\mathrm{PE}}(t)\Vert }{\Vert {\mathit{Z}}_{\mathrm{PE}}(t)\Vert +\Vert {\mathit{Z}}_{\mathrm{DP}}(t)\Vert }\right)}^3\\ a_2=1-a_1\end{array}$$

(12)

where $a_1$ and $a_2$ represent the weights of the factor $u_1$ and the factor $u_2$ respectively. Then the weight vector is expressed as: $\mathit{A}=[a_1,a_2]$.

In order to establish the relationship between weighting factors and evaluation scales, the membership degree of each factor is calculated by the non-linear membership function which is written as follows:

$$\{\begin{array}{l}{u}_1(x)={\left(k\left(x-1\right)\right)}^3\\ u_2(x)=1-{\left(k\left(x-1\right)\right)}^3\end{array}$$

(13)

where $k=0.1$, $x=1,\cdots ,11$ are the corresponding evaluation scales.

Let $r_{ij}=u_i(j)$, where $i=1,2$,$j=1,\cdots 11$, the fuzzy evaluation matrix can be obtained as:

$$\mathit{R}={[r_{ij}]}_{2\times 11}$$

(14)

The fuzzy comprehensive evaluation result vector is generated by the fuzzy synthetic operation of the weight vector and the fuzzy evaluation matrix. The fuzzy synthetic formula is defined as follows:

$$\mathit{B}=\mathit{A}\circ \mathit{R}=(b_1,b_2,\cdots ,b_{11})$$

(15)

where “$\circ$” is a fuzzy synthetic operator. In this paper, the weighted average fuzzy arithmetic operator is adopted, which can make full use of the information of $\mathit{R}$. It is specifically expressed as:

$$b_j=\min \left\{1,{\displaystyle \sum _{i=1}^2{a}_i·r_{ij}}\right\},j=1,2,\cdots ,11$$

(16)

The comprehensive evaluation value is obtained by analyzing the fuzzy comprehensive evaluation result vector. The analysis is done in the following steps. First, the result vector is normalized:

$$b_j^{\prime }=\frac{b_j}{{\displaystyle \sum _{j=1}^{11}{b}_j}}$$

(17)

Then, the normalized vector: ${\mathit{B}}^{\prime }=(b_1^{\prime },b_2^{\prime },\cdots ,b_n^{\prime })$, which is processed using the weighted average principle. The evaluation value can be obtained as follows:

$$b=\frac{{\displaystyle \sum _{j=1}^{11}{(b_j^{\prime })}^k·j}}{{\displaystyle \sum _{j=1}^{11}{(b_j^{\prime })}^k}}$$

(18)

where $k=10$ is a specific coefficient. The purpose of this coefficient is to control the role played by a larger $b_j^{\prime }(j=1,2,\cdots ,11)$. If its value increases, the role of the largest term in $b_j^{\prime }(j=1,2,\cdots ,11)$ will be more prominent.

Finally, the values of $k_1$ and $k_2$ are obtained by finding the evaluation scale corresponding to the evaluation value $b$.

3. Solution Method for Orbital Pursuit-Evasion-Defense Game

3.1. Necessary Conditions for Optimal Strategies

The orbital pursuit-evasion-defense model given in the second section can be formulated as a non-cooperative N-person differential game. Necessary conditions for optimal strategies in this game are provided by Sarma [28] and applied to the system composed of (7) and (8) to obtain the form of optimal strategies.

The Hamiltonian function is introduced as follows:

$$\{\begin{array}{l}{H}_{\mathrm{E}}={\mathit{\lambda }}_{\mathrm{E}}^{\mathrm{T}}{\dot{\mathit{Z}}}_{\mathrm{PE}}={\mathit{\lambda }}_{\mathrm{E}}^{\mathrm{T}}\left(\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{E}}{\mathit{U}}_{\mathrm{E}}\right)\\ H_{\mathrm{D}}={\mathit{\lambda }}_{\mathrm{D}}^{\mathrm{T}}{\dot{\mathit{Z}}}_{\mathrm{DP}}={\mathit{\lambda }}_{\mathrm{D}}^{\mathrm{T}}\left(\mathit{D}\mathit{\Phi }{T}_{\mathrm{D}}{\mathit{U}}_{\mathrm{D}}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}\right)\\ H_{\mathrm{P}}={\mathit{\lambda }}_{\mathrm{PE}}^{\mathrm{T}}{\dot{\mathit{Z}}}_{\mathrm{PE}}+{\mathit{\lambda }}_{\mathrm{DP}}^{\mathrm{T}}{\dot{\mathit{Z}}}_{\mathrm{DP}}={\mathit{\lambda }}_{\mathrm{PE}}^{\mathrm{T}}\left(\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{E}}{\mathit{U}}_{\mathrm{E}}\right)+{\mathit{\lambda }}_{\mathrm{DP}}^{\mathrm{T}}\left(\mathrm{D}\mathit{\Phi }{T}_{\mathrm{D}}{\mathit{U}}_{\mathrm{D}}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}\right)\end{array}$$

(19)

where ${\mathit{\lambda }}_i(i=P,D,PE,DP)$ are the co-state variables of the system.

According to the necessary conditions, the co-state equations are expressed as follows:

$$\{\begin{array}{l}{\dot{\mathit{\lambda }}}_{\mathrm{PE}}=-{\left(\frac{\partial H_{\mathrm{P}}}{\partial {\mathit{Z}}_{\mathrm{PE}}}\right)}^{\mathrm{T}}=0\\ {\dot{\mathit{\lambda }}}_{\mathrm{DP}}=-{\left(\frac{\partial H_{\mathrm{P}}}{\partial {\mathit{Z}}_{\mathrm{DP}}}\right)}^{\mathrm{T}}=0\\ {\dot{\mathit{\lambda }}}_{\mathbf{E}}=-{\left(\frac{\partial H_{\mathrm{E}}}{\partial {\mathit{Z}}_{\mathrm{PE}}}\right)}^{\mathrm{T}}=0\\ {\dot{\mathit{\lambda }}}_{\mathrm{D}}=-{\left(\frac{\partial H_{\mathrm{D}}}{\partial {\mathit{Z}}_{\mathrm{DP}}}\right)}^{\mathrm{T}}=0\end{array}$$

(20)

and the transversality conditions are written as follows:

$$\{\begin{array}{l}{\mathit{\lambda }}_{\mathrm{PE}}(t_{\mathrm{f}})=\frac{\partial J_{\mathrm{P}}}{\partial {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})}=k_1\frac{{\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})\Vert }\\ {\mathit{\lambda }}_{\mathrm{DP}}(t_{\mathrm{f}})=\frac{\partial J_{\mathrm{P}}}{\partial {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}=-k_2\frac{{\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})\Vert }\\ {\mathit{\lambda }}_{\mathrm{E}}(t_{\mathrm{f}})=\frac{\partial J_{\mathrm{E}}}{\partial {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})}=-\frac{{\mathit{Z}}_{PE}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})\Vert }\\ {\mathit{\lambda }}_{\mathrm{D}}(t_{\mathrm{f}})=\frac{\partial J_{\mathrm{D}}}{\partial {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}=\frac{{\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})\Vert }\end{array}$$

(21)

From Equations (20) and (21), we can find the following relationship:

$$\{\begin{array}{l}{\mathit{\lambda }}_{\mathrm{PE}}(t)=-k_1{\mathit{\lambda }}_{\mathrm{E}}(t)\\ {\mathit{\lambda }}_{\mathrm{DP}}(t)=-k_2{\mathit{\lambda }}_{\mathrm{D}}(t)\end{array}$$

(22)

In addition, the optimal control strategies need to satisfy:

$$\{\begin{array}{l}{\mathit{u}}_{\mathrm{D}}^{*}=\underset{\Vert {\mathit{u}}_{\mathrm{D}}\Vert \le 1}{\arg \min }{H}_{\mathrm{D}}\\ {\mathit{u}}_{\mathrm{E}}^{*}=\underset{\Vert {\mathit{u}}_{\mathrm{E}}\Vert \le 1}{\arg \min }{H}_{\mathrm{E}}\\ {\mathit{u}}_{\mathrm{P}}^{*}=\underset{\Vert {\mathit{u}}_{\mathrm{P}}\Vert \le 1}{\arg \min }{H}_{\mathrm{P}}\end{array}$$

(23)

For the sake of brevity, we define new variables ${\mathit{M}}_i(i=D,E,P)$ as follows:

$$\{\begin{array}{l}{\mathit{M}}_{\mathrm{D}}={\mathit{\lambda }}_{\mathrm{D}}^{\mathrm{T}}\mathit{D}\mathit{\Phi }{T}_{\mathrm{D}}\\ {\mathit{M}}_{\mathrm{E}}=-{\mathit{\lambda }}_{\mathrm{E}}^{\mathrm{T}}\mathit{D}\mathit{\Phi }{T}_{\mathrm{E}}\\ {\mathit{M}}_{\mathrm{P}}=-k_1{\mathit{\lambda }}_{\mathrm{E}}^{\mathrm{T}}\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}+k_2{\mathit{\lambda }}_{\mathrm{D}}^{\mathrm{T}}\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}\end{array}$$

(24)

Combining Equations (19), (22), (23), and (24) yields:

$$\{\begin{array}{l}{\mathit{u}}_{\mathrm{D}}^{*}={[u_{\mathrm{xD}}^{*},u_{\mathrm{yD}}^{*},u_{\mathrm{zD}}^{*}]}^{\mathrm{T}}=-\frac{{[{\mathit{M}}_{\mathrm{D}}(4),{\mathit{M}}_{\mathrm{D}}(5),{\mathit{M}}_{\mathrm{D}}(6)]}^{\mathrm{T}}}{\Vert {[{\mathit{M}}_{\mathrm{D}}(4),{\mathit{M}}_{\mathrm{D}}(5),{\mathit{M}}_{\mathrm{D}}(6)]}^{\mathrm{T}}\Vert }\\ {\mathit{u}}_{\mathrm{E}}^{*}={[u_{\mathrm{xE}}^{*},u_{\mathrm{yE}}^{*},u_{\mathrm{zE}}^{*}]}^{\mathrm{T}}=-\frac{{[{\mathit{M}}_{\mathrm{E}}(4),{\mathit{M}}_{\mathrm{E}}(5),{\mathit{M}}_{\mathrm{E}}(6)]}^{\mathrm{T}}}{\Vert {[{\mathit{M}}_{\mathrm{E}}(4),{\mathit{M}}_{\mathrm{E}}(5),{\mathit{M}}_{\mathrm{E}}(6)]}^{\mathrm{T}}\Vert }\\ {\mathit{u}}_{\mathrm{P}}^{*}={[u_{\mathrm{xP}}^{*},u_{\mathrm{yP}}^{*},u_{\mathrm{zP}}^{*}]}^{\mathrm{T}}=-\frac{{[{\mathit{M}}_{\mathrm{P}}(4),{\mathit{M}}_{\mathrm{P}}(5),{\mathit{M}}_{\mathrm{P}}(6)]}^{\mathrm{T}}}{\Vert {[{\mathit{M}}_{\mathrm{P}}(4),{\mathit{M}}_{\mathrm{P}}(5),{\mathit{M}}_{\mathrm{P}}(6)]}^{\mathrm{T}}\Vert }\end{array}$$

(25)

Combining Equation (25) and the form of control vector, the optimal control variables are expressed as Equation (26), which satisfies Equation (27).

$$\{\begin{array}{l}{\mathit{U}}_{\mathrm{D}}^{*}={[0,0,0,u_{\mathrm{xD}}^{*},u_{\mathrm{yD}}^{*},u_{\mathrm{zD}}^{*}]}^{\mathrm{T}}\\ {\mathit{U}}_{\mathrm{E}}^{*}={[0,0,0,u_{\mathrm{xE}}^{*},u_{\mathrm{yE}}^{*},u_{\mathrm{zE}}^{*}]}^{\mathrm{T}}\\ {\mathit{U}}_{\mathrm{P}}^{*}={[0,0,0,u_{\mathrm{xP}}^{*},u_{\mathrm{yP}}^{*},u_{\mathrm{zP}}^{*}]}^{\mathrm{T}}\end{array}$$

(26)

$$\{\begin{array}{c}{J}_{\mathrm{P}}\left({\mathit{U}}_{\mathrm{P}}^{*},{\mathit{U}}_{\mathrm{E}}^{*},{\mathit{U}}_{\mathrm{D}}^{*}\right)\le J_{\mathrm{P}}\left({\mathit{U}}_{\mathrm{P}},{\mathit{U}}_{\mathrm{E}}^{*},{\mathit{U}}_{\mathrm{D}}^{*}\right)\\ J_{\mathrm{E}}\left({\mathit{U}}_{\mathrm{P}}^{*},{\mathit{U}}_{\mathrm{E}}^{*}\right)\le J_{\mathrm{E}}\left({\mathit{U}}_{\mathrm{P}}^{*},{\mathit{U}}_{\mathrm{E}}\right)\\ J_{\mathrm{D}}\left({\mathit{U}}_{\mathrm{P}}^{*},{\mathit{U}}_{\mathrm{D}}^{*}\right)\le J_{\mathrm{D}}\left({\mathit{U}}_{\mathrm{P}}^{*},{\mathit{U}}_{\mathrm{D}}\right)\end{array}$$

(27)

Equations (10), (20), (21), and (24)–(26) constitute a TPBVP.

3.2. Hybrid Numerical Method

So far, the orbital pursuit-evasion-defense problem has been transformed into a 12-dimensional TPBVP. Generally, this kind of problem cannot be solved analytically, and numerical algorithms must be employed [9]. Numerical algorithms for solving this kind of problems include collocation method [29] and multiple shooting method [30]. The collocation method suffers from poor accuracy and high computational burden, while the multiple shooting method has high accuracy but is very sensitive to the initial guess. As noted by Pontani [9], evolutionary methods constitute an effective statistical search technique for selecting the best parameters. Thus, we apply evolutionary methods to generate the initial guess for the multiple shooting method. A hybrid method combining the multi-objective genetic algorithm and the multiple shooting method is proposed to obtain the solution of the orbital pursuit-evasion-defense game. First, the initial guesses of unknown parameters are obtained by using the multi-objective genetic algorithm. Then the exact solution of the TPBVP is solved by using the multiple shooting method.

For the sake of clarity, the state equations, the co-state equations, the initial conditions, and the terminal conditions are arranged.

Combining Equations (10), (20), and (26), the state equations and the co-state equations can be concluded as follows:

$$\{\begin{array}{l}{\dot{\mathit{Z}}}_{\mathrm{PE}}=\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}^{*}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{E}}{\mathit{U}}_{\mathrm{E}}^{*}\\ {\dot{\mathit{Z}}}_{\mathrm{DP}}=\mathit{D}\mathit{\Phi }{T}_{\mathrm{D}}{\mathit{U}}_{\mathrm{D}}^{*}-\mathit{D}\mathit{\Phi }{T}_{\mathrm{P}}{\mathit{U}}_{\mathrm{P}}^{*}\\ {\dot{\mathit{\lambda }}}_{\mathrm{E}}=0\\ {\dot{\mathit{\lambda }}}_{\mathrm{D}}=0\end{array}$$

(28)

The initial conditions of Equation (28) are expressed as follows:

$$\{\begin{array}{l}{\mathit{Z}}_{\mathrm{PE}}(0)=\mathit{D}\mathit{\Phi }(t_{\mathrm{f}},0){\mathit{X}}_{\mathrm{PE}}(0)\\ {\mathit{Z}}_{\mathrm{DP}}(0)=\mathit{D}\mathit{\Phi }(t_{\mathrm{f}},0){\mathit{X}}_{\mathrm{DP}}(0)\end{array}$$

(29)

where ${\mathit{X}}_{\mathrm{PE}}\left(0\right)={\mathit{X}}_{\mathrm{P}}\left(0\right)-{\mathit{X}}_{\mathrm{E}}(0)$, ${\mathit{X}}_{\mathrm{DP}}\left(0\right)={\mathit{X}}_{\mathrm{D}}\left(0\right)-{\mathit{X}}_{\mathrm{P}}(0)$.

According to Equation (21), the terminal conditions are written as follows:

$$\{\begin{array}{l}{\mathit{\lambda }}_{\mathrm{E}}(t_{\mathrm{f}})=\frac{\partial J_{\mathrm{E}}}{\partial {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})}=-\frac{{\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})\Vert }\\ {\mathit{\lambda }}_{\mathrm{D}}(t_{\mathrm{f}})=\frac{\partial J_{\mathrm{D}}}{\partial {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}=\frac{{\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})\Vert }\end{array}$$

(30)

3.2.1. Multi-Objective Genetic Algorithm

In the multi-objective genetic algorithm preprocessing, the terminal time $t_{\mathrm{f}}$ and the unknown initial co-state variables ${\mathit{\lambda }}_{\mathrm{E}}(0)$ and ${\mathit{\lambda }}_{\mathrm{D}}(0)$ are taken as parameters (individuals). According to the terminal conditions, the objective functions of the multi-objective genetic algorithm are set as follows:

$$\{\begin{array}{c}{J}_1=\Vert {\mathit{\lambda }}_{\mathrm{E}}(t_{\mathrm{f}})+\frac{{\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{PE}}(t_{\mathrm{f}})\Vert }\Vert \\ J_2=\Vert {\mathit{\lambda }}_{\mathrm{D}}(t_{\mathrm{f}})+\frac{{\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})}{\Vert {\mathit{Z}}_{\mathrm{DP}}(t_{\mathrm{f}})\Vert }\Vert \end{array}$$

(31)

The safe distance constraint is applied to ensure that the distance between any two players is greater than the safe distance before the terminal time. The best parameters are obtained by setting the reasonable population size, the appropriate maximum generation, and the suitable operators (i.e., crossover and mutation). The multi-objective genetic algorithm improved by Deb [31] is applied to this problem. This algorithm can reduce the complexity of computation and maintain the diversity of solutions. In this paper, we used the default operators in the toolkit on multi-objective genetic algorithm which is provided by Aravind Seshadri [32]. In addition, the population size and the number of generations are set as 100 and 200 respectively. Because of the use of the multi-objective genetic algorithm, the preprocessing time is relatively long. Thus, this algorithm is suitable for off-line calculation.

3.2.2. Multiple Shooting Method

In order to better illustrate the application of the multiple shooting method in this problem, a new state vector is defined:

$$\mathsf{\Omega }\left(t\right)=\left[{\mathit{Z}}_{\mathrm{PE}}(t),{\mathit{Z}}_{\mathrm{DP}}(t),{\mathit{\lambda }}_E(t),{\mathit{\lambda }}_{\mathrm{D}}(t)\right]$$

(32)

Substituting Equation (32) into Equation (28), the system equations can be expressed as follows:

$$\dot{\mathsf{\Omega }}\left(t\right)=f\left(t,\mathsf{\Omega }\left(t\right)\right)$$

(33)

The multiple shooting method transforms the TPBVP into a series of initial value problems. The specific steps are given as follows:

Step 1.

Divide the time interval $\left[0,t_{\mathrm{f}}\right]$ into m subintervals, and $t_k\left(k=0,\cdots ,\mathrm{m}\right)$ represents the boundary points of subintervals, which satisfy $0=t_0<t_1<\cdots <t_{\mathrm{m}}=t_{\mathrm{f}}$.

Step 2.

For each subinterval $\left[t_i,t_{i+1}\right]\left(i=0,\cdots ,\mathrm{m}-1\right)$, consider the initial value problem: $\dot{\mathsf{\Omega }}\left(t\right)=f\left(t,\mathsf{\Omega }\left(t\right)\right)$, $\mathsf{\Omega }\left(t_i\right)=s_i$, where $s_i$ is the initial value of the problem.

Step 3.

Calculate the initial guess by the multi-objective genetic algorithm.

Step 4.

Solve the initial value problem on each subinterval to obtain the solution $\mathsf{\Omega }\left(t,t_i,s_i\right)$.

Step 5.

Determine whether the condition $\mathsf{\Omega }\left(t_{i+1},t_i,s_i\right)=s_{i+1}$ and boundary conditions (29) and (30) are satisfied. If not, use the Newton method to modify the initial value and return to step 4. If the conditions are satisfied, the solution of the TPBVP is obtained successfully.

We point out that the accuracy of the initial guess value affects the solution obtained by the multiple shooting method. If the accuracy of the initial guess value is not enough, the convergence point may not be the desired solution. Moreover, it may increase the number of iterations and prolong the calculation time.

4. Results and Discussion

In this section, the following four examples are given to verify the effectiveness of the proposed strategy. Among these, Example 1 and Example 2 are taken as one group. Their initial conditions and maneuver parameters are the same. The differences between the two examples are that when performing orbital maneuvers, the pursuer in Example 1 adopts the control strategy based on the fuzzy comprehensive evaluation, while the pursuer in Example 2 does not consider the impact of the defender, that is, the parameters $k_1=1,k_2=0$ in the objective function $J_{\mathrm{P}}$. Example 3 and Example 4 are taken as the other group, with the differences between the two examples being the same as those between Example 1 and Example 2 in the first group. The initial orbital altitude of their reference orbit $h=500 \mathrm{km}$, the acceleration of gravity $g=9.8e-3\mathrm{km}/s^2$, and the radius of the Earth $R=6371.393 \mathrm{km}$. During the game, the safety distance between players is set as 0.5 km.

Example 1. The maximum unit mass thrusts of the pursuer, the evader, and the defender are $T_{\mathrm{P}}=0.09\times \mathrm{g}$, $T_{\mathrm{E}}=0.01\times \mathrm{g}$, and $T_{\mathrm{D}}=0.02\times \mathrm{g}$, respectively, and the game time is 267.4124 s. The initial positions and velocities of the pursuer, the evader, and the defender are shown in

Table 2. Positions and velocities of the initial time.

Parameter	Pursuer	Evader	Defender
$X/\mathrm{km}$	0	12	6
$Y/\mathrm{km}$	0	16	8
$Z/\mathrm{km}$	20	0	10
$V_X/\left(\mathrm{km}·{\mathrm{s}}^{-1}\right)$	0	0	0
$V_Y/\left(\mathrm{km}·{\mathrm{s}}^{-1}\right)$	0	0	0
$V_Z/\left(\mathrm{km}·{\mathrm{s}}^{-1}\right)$	0	0	0

. The pursuer adopts the control strategy based on the fuzzy comprehensive evaluation.

Figure 2 shows the curves of the positions of the three players changing with time in the directions of X, Y, Z. From Figure 2, it can be seen that the pursuer bypasses the interception of the defender and eventually catches up with the evader. From

Table 3. Position of each player at the end of the game.

Parameter	Pursuer	Evader	Defender
$X/\mathrm{km}$	15.24	15.28	10.24
$Y/\mathrm{km}$	17.44	17.68	11.72
$Z/\mathrm{km}$	−2.098	−2.363	4.87

, it can be seen that at the terminal moment, the distance between the pursuer and the evader is 0.3598 km, which is shorter than the safety distance 0.5 km, indicating that at the terminal moment, the pursuer catches up with the evader. Figure 3 shows the distance between the defender and the pursuer during the game. It reaches the shortest distance at 203.5 s. After that, the distance between the pursuer and the defender becomes longer, the shortest distance being 0.5099 km, which is longer than the safety distance 0.5 km, indicating that during the game the pursuer successfully bypasses the defender.

Figure 4 shows the curves of the control variable of each player changing with time in the directions of X, Y, Z. Figure 5 shows the curve of ZEM distance changing with time. From the figures, it can be seen that when the ZEM distance ${\mathit{Z}}_{\mathrm{DP}}(t)$ is close to 0 (i.e., 40 s to 130 s), the pursuer will consider more about evading the defender. So in this phase, the pursuer’s control curve is nearer to the control curve of the defender. During the time when the ZEM variables ${\mathit{Z}}_{\mathrm{DP}}(t)$ are not close to 0, the pursuer almost ignores the impact of the defender, so the control curve of the pursuer at this stage almost superposes with that of the evader. Through this strategy, the pursuer successfully bypasses the defenders during the game and finally captures the evader.

Example 2. The maximum unit mass thrusts of the pursuer, the evader, and the defender are $T_{\mathrm{P}}=0.09\times \mathrm{g}$, $T_{\mathrm{E}}=0.01\times \mathrm{g}$, and $T_{\mathrm{D}}=0.02\times \mathrm{g}$, respectively, and the game time is 159.81193 s. The positions and velocities of the pursuer, the evader and the defender in the initial time are shown in Table 2. The pursuer does not consider the impact of the defender when performing orbital maneuvers.

As shown in Figure 6, the defender successfully intercepts the pursuer at the terminal moment. Figure 7 shows the distance between the defender and the pursuer during the game. According to Figure 7, the distance becomes shorter and shorter in the entire game, which is caused by the pursuer’s not considering the impact of the defender. At the terminal moment, the distance between the defender and the pursuer is 0.4004 km, which is shorter than the safety distance 0.5 km.

Figure 8 shows the control variable of each player changing with time in the game. As shown in the figure, the control curve of the pursuer overlaps with that of the evader in the whole procedure, the reason being that the pursuer only considers the evader when performing orbital maneuvers.

Comparing Example 1 with Example 2, it can be seen that with the control strategy based on the fuzzy comprehensive evaluation, the pursuer can successfully bypass the defender, and finally capture the evader. The pursuer, for not considering the impact of the defender, is eventually intercepted by the defender.

Example 3. The maximum unit mass thrusts of the pursuer, the evader, and the defender are $T_{\mathrm{P}}=0.14\times \mathrm{g}$, $T_{\mathrm{E}}=0.01\times \mathrm{g}$, and $T_{\mathrm{D}}=0.14\times \mathrm{g}$, respectively. The maneuverability of the defender and that of the pursuer are the same, and the game time is 177.87788 s. The positions and velocities of the pursuer, the evader, and the defender in the initial time are shown in

Table 4. Positions and velocities of the initial time.

Parameter	Pursuer	Evader	Defender
$X/\mathrm{km}$	0	8	18
$Y/\mathrm{km}$	0	9	24
$Z/\mathrm{km}$	30	12	0
$V_X/\left(\mathrm{km}·{\mathrm{s}}^{-1}\right)$	0	0	0
$V_Y/\left(\mathrm{km}·{\mathrm{s}}^{-1}\right)$	0	0	0
$V_Z/\left(\mathrm{km}·{\mathrm{s}}^{-1}\right)$	0	0	0

. The pursuer adopts the control strategy based on the fuzzy comprehensive evaluation.

Figure 9 shows the curves of the positions of the three players changing with time in directions of X, Y, Z. As shown in Figure 9, the defender successfully intercepts the pursuer at the end of the game. Figure 10 shows the curves of the control variable of each player changing with time in the directions of X, Y, Z. At 160 s or so, the pursuer starts to change the control strategy to evade the defender. However, because of the same maneuverability of the defender and the pursuer, the pursuer does not successfully bypass the interception of the defender.

Table 5. Position of each player at the end of the game.

Parameter	Pursuer	Evader	Defender
$X/\mathrm{km}$	8.293	9.083	8.638
$Y/\mathrm{km}$	8.804	9.609	9.026
$Z/\mathrm{km}$	11.89	10.51	11.79

shows the position of each player at the terminal moment. From Table 5, it can be seen that at the terminal moment, the distance between the defender and the pursuer is 0.4223 km, which is shorter than the safety distance 0.5 km, and at the terminal moment, the distance between the pursuer and the evader is 1.7823 km, which is longer than the safety distance 0.5 km. All the above show that at the terminal moment, the defender successfully intercepts the pursuer and that the evader successfully evades the capture of the pursuer.

Example 4. The maximum unit mass thrusts of the pursuer, the evader, and the defender are $T_{\mathrm{P}}=0.14\times \mathrm{g}$, $T_{\mathrm{E}}=0.01\times \mathrm{g}$, $T_{\mathrm{D}}=0.14\times \mathrm{g}$, respectively. The maneuverability of the defender and that of the pursuer are the same, and the game time is 177.19431 s. The positions and velocities of the pursuer, the evader, and defender in the initial time are shown in Table 4. The pursuer does not consider the impact of the defender when performing orbital maneuvers.

As shown in Figure 11, the defender intercepts the pursuer at the terminal moment. From

Table 6. Position of each player at the end of the game.

Parameter	Pursuer	Evader	Defender
$X/\mathrm{km}$	8.384	9.069	8.714
$Y/\mathrm{km}$	8.956	9.593	9.132
$Z/\mathrm{km}$	11.81	10.51	11.68

, it can be seen that the distance between the pursuer and the defender in the game is 0.3959 km, which is shorter than the safety distance 0.5 km, and that the distance between the pursuer and the evader is 1.6016 km, which is longer than the safety distance 0.5 km. This shows that the defender intercepts the pursuer successfully at the terminal moment, and the evader evades the capture of the pursuer successfully. Figure 12 shows the curves of the control variable of each player changing with time in the directions of X, Y, Z. From the figure, it can be seen that the control curves of the pursuer remain overlapped with those of the evader.

Comparing Example 3 with Example 4, it can be seen that, because of the different control strategies adopted by the pursuer, the time that the defender takes to intercept the pursuer in Example 3 is longer than that in Example 4. Moreover, at the terminal moment, the distance between the pursuer and the defender in Example 4 is shorter than that in Example 3.

The comparison between Example 1 and Example 2 shows that when the control variable of the pursuer is in a dominant position, the optimal control strategy proposed in this paper makes the pursuer bypass the defender and capture the evader. The comparison between Example 3 and Example 4 shows that when the control variable of the pursuer is not in a dominant position, the optimal control strategy proposed in this paper prolongs the time that the defender takes to intercept the pursuer.

5. Conclusions

The fuzzy comprehensive evaluation and the differential game theory are applied to design the control strategy of the pursuer in the orbital pursuit-evasion-defense problem. The hybrid method combining the multi-objective genetic algorithm and the multiple shooting method is proposed to solve the problem. The simulation results show that when the pursuer control is in a dominant position, the control strategy proposed in this paper can make the pursuer bypass the defender and capture the evader, and that when the pursuer control is not in a dominant position, the control strategy proposed in this paper can prolong the time that the defender takes to intercept the pursuer. The proposed control strategy is applicable to the orbital pursuit-evasion-defense scenario, in which the players adopt the continuous low thrust propulsion. When the ZEM distance between the pursuer and the defender is close to zero, the control strategy can be automatically switched to parallel with the defender’s control strategy, so that the pursuer can effectively avoid the interception of the defender.

However, the limitation of this paper is that the terminal time of the game is given by the genetic algorithm, which is not accurate. Further research will be carried out on the accurate calculation of the terminal time.