2.2 System model

2.2.1 dynamical model

The UAV dynamics are modeled as Equation (1):

(1) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\dot x} { = V\,{\rm{cos\,\Theta\, cos\,\Psi }}}\\[3pt]{\dot y} { = V\,{\rm{sin\,\Theta }}}\\[3pt]{\dot z} { = V\,{\rm{cos\,\Theta \,sin\,\Psi }}}\\[5pt]{\dot V} { = - \dfrac{D}{m} - g\,{\rm{sin\,\Theta }}}\\[12pt]{{\dot{\Theta}}} { = \dfrac{{L\,{\rm{cos}}\,\nu }}{{mV}} - \dfrac{{g\,{\rm{cos\,\Theta }}}}{V}}\\[12pt]{{\dot{\Psi}}} { = \dfrac{{L\,{\rm{sin}}\,\nu }}{{mV\,{\rm{cos\,\Theta }}}}}\end{array}} \right.\end{align}

Where $V$ is the magnitude of velocity, ${\rm{\Theta }}$ is the trajectory angle, ${\rm{\Psi }}$ is the trajectory yaw angle, $x,y,z$ is the coordinates in the North Celestial Inertial System, $D$ is the aerodynamic drag, $L$ is the aerodynamic lift, $\nu $ is the velocity pitch angle, $g$ is the local gravitational acceleration, and $m$ is the mass.

2.2.2 Relative motion model

In the differential game model, we simplify the target coordinate system. First, we switch from the Northeastern system to the one in Fig. 1. It’s split into the firing and horizontal planes, as shown in Fig. 2.

Figure 2. Planar motion decoupling.

Where ${\lambda _{Di}}$ is the line-of-sight elevation angle, ${V_{Dmi}},{V_{Ddi}}$ are the velocitys projection of the firing plane $M$ and ${D_i},\left( {i = 1,2} \right)$ , ${\gamma _{Dmi}},{\gamma _{Ddi}}$ are the firing plane velocity angles, ${\phi _{Dmi}},{\phi _{Ddi}}$ are the front angles of the firing planes of escapee $M$ and pursuer ${D_i}$ . ${\lambda _{Ti}}$ is the line-of-sight azimuth, ${V_{Tmi}},{V_{Tdi}}$ are the velocity projections of M and ${D_i}\!\left( {i = 1,2} \right)$ in the horizontal plane, ${\gamma _{Tmi}},{\gamma _{Tdi}}$ are the horizontal plane velocity angles, and ${\phi _{Tmi}},{\phi _{Tdi}}$ are the front angles of the horizontal plane M and ${D_i}\!\left( {i = 1,2} \right)$ .

Modeling of relative motion in the firing plane is as Equation (2):

(2) \begin{align}\left\{ {\begin{array}{l}{{\phi _{Dmi}}} { = {\gamma _{Dmi}} - {\lambda _{Di}}}\\[4pt]{{\phi _{Ddi}}} { = {\gamma _{Ddi}} - {\lambda _{Di}}}\\[4pt]{{{\dot R}_{Di}}} { = {V_{Dmi}}\,{\rm{cos}}\,{\phi _{Dmi}} - {V_{Ddi}}\,{\rm{cos}}\,{\phi _{Ddi}}}\\[4pt]{{R_{Di}}{{\dot \lambda }_{Di}}} { = {V_{Dmi}}\,{\rm{sin}}\,{\phi _{Dmi}} - {V_{Ddi}}{\phi _{Ddi}}}\end{array}} \right.\end{align}

Where ${R_{Di}}$ is the distance between M and ${D_i}\!\left( {i = 1,2} \right)$ .

Let $\frac{{{V_{Dm1}}}}{{{V_{Dd1}}}} = {\alpha _{D1}},\frac{{{V_{Dm2}}}}{{{V_{Dd2}}}} = {\alpha _{D2}}$ , then let the transformed relative distance be as Equation (3):

(3) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{{\bar R}_{D1}}} { = {\alpha _{D1}}\,{\rm{cos}}\,{\phi _{Dm1}} - \,{\rm{cos}}\,{\phi _{Dd1}}}\\[5pt]{{{\bar R}_{D2}}} { = {\alpha _{D2}}\,{\rm{cos}}\,{\phi _{Dm2}} - \,{\rm{cos}}\,{\phi _{Dd2}}}\end{array}} \right.\end{align}

Let the transformed relative distance ${\bar R_{D1}},{\bar R_{D2}}$ be the state variable ${x_1},{x_2}$ for the differential game model and the angular velocity of the plane ${\dot \gamma _{Dmi}},{\dot \gamma _{Ddi}}$ be the control variable ${u_{Dmi}},{u_{Ddi}}$ , then the system state is as Eqaution 4:

(4) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{{\dot x}_1}} { = {f_{D1}}({x,t}){u_{Dm1}} + {g_{D1}}({x,t}){u_{Dd1}} + {h_{D1}}({x,t})}\\[3pt]{{{\dot x}_2}} { = {f_{D2}}({x,t}){u_{Dm2}} + {g_{D2}}({x,t}){u_{Dd2}} + {h_{D2}}({x,t})}\end{array}} \right.\end{align}

Where ${f_{Di}}({x,t}) = - {\alpha _{Di}}\,{\rm{sin}}\,{\phi _{Dmi}}$ , ${g_{Di}}({x,t}) = \,{\rm{sin}}\,{\phi _{Ddi}}$ , ${h_{Di}}({x,t}) = {\dot \lambda _{Di}}\left( {{\alpha _{Di}}\,{\rm{sin}}\,{\phi _{Dmi}} - \,{\rm{sin}}\,{\phi _{Ddi}}} \right)$ , $\left( {i = 1,2} \right)$ .

${\dot \gamma _{Dmi}}$ are synthesised as the total firing plane control variable as the angular velocity command ${u_{Dm}}$ for the trajectory angle, i.e., ${\dot{\Theta}} = {u_{Dm}}$ . The same is true for the horizontal plane. But, the control variable there is the angular velocity command for the trajectory yaw, so ${\dot{\Psi}} = {u_{Tm1}} = {u_{Tm2}} = {u_{Tm}}$ . The control variable in the horizontal plane is the angular velocity command for the trajectory yaw.

Combining Equation (1) gives the lift as Equation (5):

(5) \begin{align}L = \sqrt {{{(mV{u_{Dm}} + mg\,{\rm{cos\,\Theta }})}^2} + {{(mV\,{\rm{cos\,\Theta }}{u_{Tm}})}^2}} \end{align}

The overload is $N = \frac{L}{{mg}}$ , and after the overload constraint, the real aerodynamic lift $L$ and aerodynamic drag $D$ are obtained using atmospheric interpolation, and the real trajectory angular velocity ${\dot{\Theta}}$ and trajectory yaw angular velocity ${\dot{\Psi}}$ are obtained from Equation (1).

3.1 Performance indicators

The firing plane and the horizontal plane are decoupled. From the perspective of evasion, the horizontal plane is taken as an example, and the design performance index is:

(6) \begin{align}{\min\limits_{{u_{Tm}}} J = } { - \frac{1}{2}{\beta _1}\left( {{t_{f1}}} \right)\bar R_{T1}^2\left( {{t_{f1}}} \right) - \frac{1}{2}{\beta _2}\left( {{t_{f2}}} \right)\bar R_{T2}^2\left( {{t_{f2}}} \right) + \frac{1}{2}{\kappa _T}\int u_{Tm}^2dt}\end{align}

Where ${\beta _i}(t) = \frac{{{V_{Tdi}}(t)}}{{{V_{d1}}(t) + {V_{d2}}(t)}}$ is the weight of the terminal value index, ${t_{fi}}\left( {i = 1,2} \right)$ is the terminal moment and ${\bar R_{T1}},{\bar R_{T2}}$ is the transformed relative distance from the escapee and pursuer on the horizontal plane corresponding to the transformation in (3). ${\kappa _T}$ is the weight of the energy metric.

The performance index aims to increase the distance between the escapee and pursuer. It also aims to reduce the escapee’s energy use. This index shows the evasion and energy. The parameter ${\beta _i}(t)$ is modifiable to dynamically adjust according to the capabilities of the pursuer, thereby augmenting the significance of the threat of a proficient pursuer. This adaptive mechanism serves to mitigate the inherent limitation of relying solely on distance metrics for evaluating the dynamics of the game. Parameter ${\kappa _T}$ based on engineering experience.

To represent the amount of last-minute avoidance, which can’t be estimated or predicted, it is converted to an integral index [Reference Xing30]:

(7) \begin{align}\frac{1}{2}{\beta _i}\bar R_{Ti}^2\left( {{t_{fi}}} \right) & = {\frac{1}{2}{\beta _{i0}}\bar R_{Ti}^2\left( {{t_0}} \right) + \frac{1}{2}\int \nolimits_{{t_0}}^{{t_{fi}}} d\!\left( {{\beta _i}\bar R_{Ti}^2} \right)dt}\nonumber\\[5pt]& = {\frac{1}{2}{\beta _{i0}}\bar R_{Ti}^2\left( {{t_0}} \right) + \int \nolimits_{{t_0}}^{{t_{fi}}} \,\left( {{\beta _i}{x_i}{{\bar R}_{Ti}} + \frac{1}{2}{{\dot \beta }_i}\bar R_{Ti}^2} \right)dt}\end{align}

Thereby the performance index is converted to:

(8) \begin{align}\min\limits_{{u_{Tm}}} J & = - \frac{1}{2}{\beta _{10}}\bar R_{T10}^2 - \frac{1}{2}{\beta _{20}}\bar R_{T20}^2 - \int \nolimits_{{t_0}}^{{t_f}} \left[{l_1}\left( {{\beta _1}{x_1}\bar R_{T1}^2 + \frac{1}{2}{{\dot \beta }_1}\bar R_{T1}^2} \right)\right.\nonumber\\[5pt]& \quad \left.+ {l_2}\left({\beta _2}{x_2}\bar R_{T2}^2 + \frac{1}{2}{{\dot \beta }_2}\bar R_{T2}^2\right)\right]dt + {\kappa _T}\int \nolimits_{{t_0}}^{{t_f}} \frac{1}{2}u_T^2dt\end{align}

where ${\beta _{10}},{\beta _{20}}$ is the value of ${\beta _1},{\beta _2}$ taken at the initial moment and ${\bar R_{T10}},{\bar R_{T20}}$ the initial distance between the escapee M and the pursuer ${D_i}\!\left( {i = 1,2} \right)$ after the transformation. ${t_f} = {\rm{max}}\!\left\{ {{t_{f1}},{t_{f2}}} \right\}.$ ${l_1} = \left\{ \begin{array}{l@{\quad}l}{1,} & {t \lt {t_{f1}}}\\[4pt]{0,} & {t \geqslant {t_{f1}}}\end{array} \right.,{l_2} = \left\{ \begin{array}{l@{\quad}l}{1,} & {t \lt {t_{f2}}}\\[4pt]{0,} & {t \geqslant {t_{f2}}}\end{array} \right.$ are the time operators.

This procedural change streamlines the solution process. It maintains the integrity of the solution and avoids adding errors.

3.2 HJB equations solving

The Hamiltonian function is:

(9) \begin{align}H & = - {l_1}\left( {{\beta _1}{x_1}{{\bar R}_{T1}} + \frac{1}{2}{{\dot \beta }_1}\bar R_{T1}^2} \right) - {l_2}\!\left( {{\beta _2}{x_2}{{\bar R}_{T2}} + \frac{1}{2}{{\dot \beta }_2}\bar R_{T2}^2} \right) + \frac{1}{2}{\kappa _T}u_{Tm}^2\nonumber\\[4pt] & \quad + {\lambda _1}\!\left( {{f_{T1}}({x,t}){u_{Tm}} + {g_{T1}}({x,t}){u_{Td1}} + {h_{T1}}({x,t})} \right)\nonumber\\[4pt] & \quad + {\lambda _2}\!\left( {{f_{T2}}({x,t}){u_{Tm}} + {g_{T2}}{u_{Td2}} + {h_{T2}}({x,t})} \right) \end{align}

According to the reference [Reference Bressan and Shen31], let the

(10) \begin{align}\phi & = - {l_1}\!\left( {{\beta _1}{x_1}\bar R_{T1}^2 + \frac{1}{2}{{\dot \beta }_1}\bar R_{T1}^2} \right) - {l_2}\!\left( {{\beta _2}{x_2}{{\bar R}_{T2}} + \frac{1}{2}{{\dot \beta }_2}\bar R_{T2}^2} \right)\nonumber\\[3pt]& \quad + \frac{1}{2}{\kappa _T}u_{Tm}^2 + {\lambda _1}\!\left( {{f_{T1}}({x,t}){u_{Tm}}} \right) + {\lambda _2}\!\left( {{f_{T2}}({x,t}){u_{Tm}}} \right) \end{align}

Where ${u_{Tm}}$ reaches the optimal solution.

(11) \begin{align}\frac{{\partial \phi }}{{\partial {u_{Tm}}}} & = {\frac{{ - {l_1}\!\left( {{\beta _1}{x_1}\bar R_{T1}^2 + \frac{1}{2}{{\dot \beta }_1}\bar R_{T1}^2} \right) - {l_2}\!\left( {{\beta _2}{x_2}{{\bar R}_{T2}} + \frac{1}{2}{{\dot \beta }_2}\bar R_{T2}^2} \right)}}{{\partial {u_{Tm}}}}} + {\kappa _T}{u_{Tm}} + {\lambda _1}\!\left( {{f_{T1}}({x,t})} \right) + {\lambda _2}\!\left( {{f_{T2}}({x,t})} \right)\end{align}

Taylor expands $\overline {{R_i}}, {x_i}$ to obtain

(12) \begin{align}\left\{\begin{array}{l} {{{\bar R}_i}} = {{\bar R}_{i0}} + {x_i}t + \frac{1}{2}{{\dot x}_i}{t^2} + o({t^3})\\[5pt]{{x_i}} = {x_{i0}} + \dot xt + o({t^2}) \end{array}\right.\end{align}

solving for ${u_{Tm}}$ is

(13) \begin{align}{u_{Tm}} & = \frac{1}{{{\kappa _T}}}\left({\beta _1}{f_{T1}}({x,t}){{\bar R}_1}{l_1}t + \frac{1}{2}{f_{T1}}({x,t})\left( {{{\dot \beta }_1}{{\bar R}_1} + {\beta _1}{x_1}} \right){l_1}{t^2} + {\beta _2}{f_{T2}}({x,t}){{\bar R}_2}{l_2}t\right.\nonumber\\[4pt] &\quad \left.+ \frac{1}{2}{f_{T2}}({x,t})\left( {{{\dot \beta }_2}{{\bar R}_2} + {\beta _2}{x_2}} \right){l_2}{t^2}\right) - \left( {{\lambda _1}{f_{T1}}({x,t}) + {\lambda _2}{f_{T2}}({x,t})} \right)\end{align}

3.3 Solving for covariates using Cauchy’s initial value

3.3.1 Iterative solution method

Using the regular equation, Equation (14) is obtained:

(14) \begin{align}\left\{\begin{array}{l}{{{\dot \lambda }_1}({x,t})} = - \dfrac{{\partial H}}{{\partial {x_1}}} = {\beta _1}{{\bar R}_1}{l_1}\\[9pt]{{{\dot \lambda }_1}({x,t})} = - \dfrac{{\partial H}}{{\partial {x_1}}} = {\beta _1}{{\bar R}_1}{l_1}\end{array}\right.\end{align}

Let ${\dot \lambda _i}\!\left( {x,{t_0}} \right) = {\beta _0}{\bar R_{i0}}$ , from Equation (12), and Equation (14) we know that

(15) \begin{align}{\dot \lambda _i} = {\beta _i}{l_i}\!\left( {{{\bar R}_{i0}} + {x_i}t + \frac{1}{2}{{\dot x}_i}{t^2},{{\ddot \lambda }_i} = {\beta _i}{l_i}\!\left( {{x_i} + {{\dot x}_i}t} \right)} \right)\end{align}

Then, the Cauchy initial value problem is satisfied:

(16) \begin{align}{\ddot \lambda _i} + p(t){\dot \lambda _i} = {q_i}({x,t})\end{align}

Where $p(t) = \frac{2}{t},{q_i}({x,t}) = - {\beta _i}{l_i}\!\left( {{x_i} + \frac{{2{{\bar R}_{i0}}}}{t}} \right)\left( {i = 1,2} \right)$ .

Referring to the reference (Reference Youness, Megahed, Eladdad and Madkour32), it is deduced that:

(17) \begin{align}{{\lambda _i}({x,t}) = } {\frac{1}{{p(t)}}\left( {I\!\left( {{q_i}({x,t})} \right) + I\!\left( {{\lambda _i}({x,t})\dot p(t)} \right) - {{\dot \lambda }_i}({x,t}) + {\beta _{i0}}{{\bar R}_{i0}}} \right)}\end{align}

An iterative approach was used to calculate ${\lambda _i}({x,t})$

(18) \begin{align}\left\{\begin{array}{l}\lambda _i^{(n)}({x,t}) = \dfrac{1}{{p(t)}}\left( {I\!\left( {{q_i}({x,t})} \right) + I\!\left( {\lambda _i^{({n - 1})}({x,t})\dot p(t)} \right) - {{\dot \lambda }_i}({x,t}) + {\beta _{i0}}{{\bar R}_{i0}}} \right)\\[8pt]{\lambda ^{\left( 0 \right)}}({x,t}) = - \dfrac{t}{2}\left( {{{\bar R}_i}{\beta _i} + {\beta _{i0}}{{\bar R}_{i0}}} \right){l_i}\end{array} \right.\end{align}

Taking the value of the 5th iteration, ${\lambda _i}$ is obtained:

(19) \begin{align}{\lambda _i} & = {l_i}\!\left(\frac{1}{{48}}{{\bar R}_{i0}}{\beta _{i0}}t\,{\rm{l}}{{\rm{n}}^4}\,t + \frac{1}{{48}}{{\bar R}_i}{\beta _i}t\,{\rm{l}}{{\rm{n}}^4}\,t - \frac{1}{6}{{\bar R}_i}{\beta _i}t\,{\rm{l}}{{\rm{n}}^3}\,t + \frac{1}{6}{{\bar R}_{i0}}{\beta _{i0}}t\,{\rm{l}}{{\rm{n}}^3}\,t\right.\nonumber\\[3pt]& \quad + \frac{1}{6}{{\bar R}_{i0}}{\beta _i}\,{\rm{l}}{{\rm{n}}^3}\,t + \frac{1}{2}{{\bar R}_i}{\beta _i}t\,{\rm{l}}{{\rm{n}}^2}\,t - \frac{1}{4}{{\bar R}_{i0}}{\beta _{i0}}t\,{\rm{l}}{{\rm{n}}^2}\,t + \frac{1}{2}{{\bar R}_{i0}}{\beta _i}t\,{\rm{l}}{{\rm{n}}^2}\,t\nonumber\\[3pt]& \quad \left.- {{\bar R}_i}{\beta _i}t\,{\rm{ln}}\,t + \frac{1}{2}{{\bar R}_{i0}}{\beta _{i0}}t\,{\rm{ln}}\,t + {{\bar R}_{i0}}{\beta _i}t\,{\rm{ln}}\,t + {{\bar R}_i}{\beta _i}t - \frac{1}{2}{{\bar R}_{i0}}{\beta _{i0}}t\right)\end{align}

3.3.2 Proof of convergence

From Equation (18), we know that

(20) \begin{align}\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t}) = {\frac{1}{{p(t)}}\left( {I\!\left( {{\lambda ^{({n - 1})}}({x,t})\dot p(t)} \right) - I\!\left( {\lambda _i^{\left( {n - 2} \right)}({x,t})\dot p(t)} \right)} \right)}\end{align}

When n = 1, from Equations (17) and (18):

(21) \begin{align}\lambda _i^{(1)}({x,t}) - \lambda _i^{\left( 0 \right)}({x,t}) & = \frac{1}{{p(t)}}[I\!\left( {q({x,t})} \right) + I\!\left( {\lambda _i^{\left( 0 \right)}({x,t})\dot p(t)} \right) - 2{{\bar R}_i}\beta {l_i}\nonumber\\[5pt]& = - \frac{t}{2}\left[ {I\!\left( {q({x,t})} \right) - \int \,\frac{{{{\bar R}_i}{\beta _i} + {{\bar R}_{i0}}{\beta _{i0}}}}{t}{l_i}dt - 2{{\bar R}_i}{\beta _i}{l_i}} \right]\end{align}

Since ${\bar R_i},{\bar R_{i0}},{\beta _i},{\beta _{i0}} \gt 0,{t_0} \gt 0,0 \lt t \lt {t_f}$ and ${x_i},\frac{{2\overline {{R_{i0}}} }}{t}$ are bounded, Equation (22) is deflated to obtain

(22) \begin{align}{\left| {\lambda _i^{(1)}({x,t}) - \lambda _i^{\left( 0 \right)}({x,t})} \right|} { \le \frac{t}{2}\left( {\left| {I\!\left( {q({x,t})} \right)} \right| + \int \frac{{{{\bar R}_i}{\beta _i} + {{\bar R}_{i0}}{\beta _{i0}}}}{t}{l_i}dt + 2{{\bar R}_i}{\beta _i}{l_i}} \right) \le Qt}\end{align}

where Q is the upper bound of $\left( {\left| {I\!\left( {q({x,t})} \right)} \right| + \int \frac{{{{\bar R}_i}{\beta _i} + {{\bar R}_{i0}}{\beta _{i0}}}}{t}{l_i}dt + 2{{\bar R}_i}{\beta _i}{l_i}} \right)$

When n = 2, from Equation (18):

(23) \begin{align}{\left| {\lambda _i^{(2)}({x,t}) - \lambda _i^{(1)}({x,t})} \right|} { = \left| {\frac{1}{{p(t)}}\left( I \right))\lambda _i^{(1)}({x,t}) - \lambda _i^{\left( 0 \right)}({x,t}))\dot p(t))|} \right. \le \left| { - \frac{t}{2} \cdot Qt\int \frac{2}{{{t^2}}}dt} \right| = Qt}\end{align}

When n = 3, from Equations (18) and (23):

(24) \begin{align}{\left| {\lambda _i^{\left( 3 \right)}({x,t}) - \lambda _i^{\left( 2 \right)}({x,t})} \right|} { = \left| {\frac{1}{{p(t)}}\left(I((\lambda _i^{(2)}({x,t}) - \lambda _i^{(1)}({x,t}))\dot p(t))\right)} \right| \le Qt}\end{align}

And so on to get $\left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right| \le Qt$ , so $\left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right|$ bounded.

When $n \geqslant 3$ ,

(25) \begin{align}\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t}) & = {\frac{1}{{p(t)}}\left( {I\!\left( {\lambda _i^{({n - 1})}({x,t})\dot p(t)} \right) - I\!\left( {\lambda _i^{\left( {n - 2} \right)}({x,t})\dot p(t)} \right)} \right)}\nonumber\\& = {\frac{1}{{p(t)}}\left( {p(t)\left( {\lambda _i^{({n - 1})}({x,t}) - \lambda _i^{\left( {n - 2} \right)}({x,t})} \right)} \right)}\nonumber\\& - \frac{1}{{p(t)}}\int {{\left( {\frac{1}{{p(t)}}} \right)}^{\rm{'}}}\dot p(t)\left( {\lambda _i^{\left( {n - 2} \right)}({x,t}) - \lambda _i^{\left( {n - 3} \right)}({x,t})} \right)p(t)dt\nonumber\\& = \left( {\lambda _i^{({n - 1})}({x,t}) - \lambda _i^{\left( {n - 2} \right)}({x,t})} \right) + \frac{t}{2}\int \frac{2}{{{t^3}}}\left( {\lambda _i^{\left( {n - 2} \right)}({x,t}) - \lambda _i^{\left( {n - 3} \right)}({x,t})} \right)dt\end{align}

There exists $\xi \in \left( {{t_0},t} \right)$ , which is obtained by using the median theorem of the integral

(26) \begin{align}{\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} & = \left( {\lambda _i^{({n - 1})}({x,t}) - \lambda _i^{\left( {n - 2} \right)}({x,t})} \right) + \frac{t}{2} \cdot \frac{2}{{t_0^3}}\int \nolimits_{{t_0}}^\xi \left( {\lambda _i^{\left( {n - 2} \right)}({x,t}) - \lambda _i^{\left( {n - 3} \right)}({x,t})} \right)dt\nonumber\\[4pt]& = \left( {\lambda _i^{({n - 1})}({x,t}) - \lambda _i^{\left( {n - 2} \right)}({x,t})} \right) + \frac{t}{{t_0^3}}\int \nolimits_{{t_0}}^\xi \left( {\lambda _i^{\left( {n - 2} \right)}({x,t}) - \lambda _i^{\left( {n - 3} \right)}({x,t})} \right)dt\end{align}

Combining Equation (24) that

(27) \begin{align}\left| {\frac{t}{{t_0^3}}\int \nolimits_{{t_0}}^\xi \left( {\lambda _i^{\left( {n - 2} \right)}({x,t}) - \lambda _i^{\left( {n - 3} \right)}({x,t})} \right)dt} \right| & = \left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t}) - \left( {\lambda _i^{({n - 1})}({x,t}) - \lambda _i^{\left( {n - 2} \right)}({x,t})} \right)} \right|\nonumber\\[4pt]& \le \left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right| + \left| {\lambda _i^{({n - 1})}({x,t}) - \lambda _i^{\left( {n - 2} \right)}({x,t})} \right|\nonumber\\[4pt]& \le 2Qt\end{align}

Therefore, $\left| {{\frac{1}{{{t_0}}}}^3}\int \nolimits_{{t_0}}^\xi {\left( \lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})\right)dt} \right|$ is bounded, and let the upper bound of $\left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right|$ be S. Then

(28) \begin{align}\left| {\frac{1}{{t_0^3}}\int \nolimits_{{t_0}}^\xi \left( {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)dt} \right| \le \left| {\frac{1}{{{t_0}^3}}S\left( {\xi - {t_0}} \right)} \right| = \frac{{S\!\left( {\xi - {t_0}} \right)}}{{t_0^3}}\end{align}

There exists $L \gt 0$ such that $\frac{{S\!\left( {\xi - {t_0}} \right)}}{{t_0^3}} \le L$ , then $S \le \frac{{Lt_0^3}}{{\left( {\xi - {t_0}} \right)}}$ , when $n \to \infty $ , ${t_0} \to 0$ , $S \to 0$ , so $\left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right| \to 0$ , $\lambda _i^{(n)}({x,t})$ converge.

3.3.3 Proof of uniqueness

The following is a proof of uniqueness by contradiction.

Known from 3.3.2, $\lambda _i^{\left( 0 \right)}$ is bounded, $\left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right|$ is bounded and ${\lambda ^{(n)}}({x,t})$ converges, so ${\lambda ^{(n)}}({x,t})$ is bounded.

Assuming that the solution of ${\lambda _i}({x,t})$ is not unique, there exists a bounded function $\mu ({x,t})$ that is a solution of ${\lambda _i}({x,t})$ , and there exists $N \gt 0$ that satisfies

(29) \begin{align}\left\{\begin{array}{l}\mu ({x,t}) = \dfrac{1}{{p(t)}}\left( {I\!\left( {q({x,t})} \right) + I\!\left( {\mu ({x,t})\dot p(t)} \right) - {{\dot \lambda }_t}({x,t}) + {\beta _{i0}}{{\bar R}_{i0}}} \right)\\[11pt]0 \lt \left| {\mu ({x,t}) - \lambda _t^{(n)}({x,t})} \right| \lt N\end{array} \right.\end{align}

Then

(30) \begin{align}\mu ({x,t}) - \lambda _i^{(n)}({x,t}) & = \frac{1}{{p(t)}}\left( {I\!\left( {\mu ({x,t})\dot p(t)} \right) - I\!\left( {\lambda _i^{({n - 1})}({x,t})\dot p(t)} \right)} \right)\nonumber\\[4pt]& = - \frac{t}{2}\int \nolimits_{{t_0}}^t \frac{2}{{{t^2}}}\left( {\mu ({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)dt\end{align}

By the integral median theorem, there exists $\zeta \in \left( {{t_0},t} \right)$ such that

(31) \begin{align}{\mu ({x,t}) - \lambda _i^{(n)}({x,t})} & = - \frac{t}{2}\int \nolimits_{{t_0}}^t \frac{2}{{{t^2}}}\left( {\mu ({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)dt\nonumber\\[4pt]& = - \frac{t}{{t_0^2}}\int\limits^{\!\!\!\zeta} \left( {\mu ({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)dt\end{align}

Equation (30) deflates to give:

(32) \begin{align}N \gt \left| {\mu ({x,t}) - \lambda _i^{(n)}({x,t})} \right| = \left| {\frac{t}{{t_0^2}}\int \nolimits_{{t_0}}^\zeta \left( {\mu ({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)dt} \right|\end{align}

Let $0 \lt L \le N$ , $L$ is the lower bound of $\left( {\mu ({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)$ . So when $n \to \infty $ , there are

(33) \begin{align}{\left| {\mu ({x,t}) - \lambda _i^{(n)}({x,t})} \right|} { = \left| {\frac{t}{{t_0^2}}\int \nolimits_{{t_0}}^\zeta \left( {\mu ({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right)dt} \right| \geqslant \frac{{Lt\!\left( {\zeta - {t_0}} \right)}}{{_0^2}}}\end{align}

Therefore

(34) \begin{align}L \le \frac{{Nt_0^2}}{{t\!\left( {\zeta - {t_0}} \right)}}\end{align}

When $n \to \infty, {t_0} \to 0,{t_0} \lt t \lt {t_{fj}}$ , then $L \to 0$ , which contradicts $0 \lt L \le N$ . Therefore, the assumption is not valid and the solution of ${\lambda _i}({x,t})$ is unique.

Table 1. Simulation parameterst (normalised)

Figure 4. Comparison of experimental trajectories with method 1.

This section details the development of an optimal avoidance strategy for UAVs based on a nonlinear dynamical model. It demonstrates the convergence and uniqueness of the key variables, known as covariates, within this strategy. Consequently, the proposed approach offers a rigorously validated solution to the nonlinear optimal strategy, addressing the shortcomings of traditional analytical methods. This advancement holds significant implications for improving UAV evasion tactics.

4.1 Simulation condition

The simulation experiment uses the non-cooperative differential game given in this paper for a game scenario with two pursuers D1, D2 and one escapee M. The simulation stops when the escapee evades pursuit and reaches the destination. Or, it stops if the escapee is not pursued. The shooting plane energy metrics are weighted ${\kappa _D} = 10$ , and the horizontal plane energy metrics are weighted ${\kappa _T} = 100$ . The destination is set at the origin (0, 0, 0), and c is a constant. The simulation parameters are set as in Table 1.

Figure 5. Comparison of experimental trajectories with method 2.

Figure 6. Comparison of experimental trajectories with method 3.

4.2 Simulation results

The escapee M uses the strategy designed in this paper, and the pursuers D1 and D2 use proportional guidance. Under Table 1, the evasion is successful, M evades 112.633c against D1, and 2.533c against D2. The final distance between M and T is 0.133c, the final velocity of M is 0.192, and the time when M reaches T is 44.180.

According to the Apollo circle method described in reference (Reference Makkapati, Sun and Tsiotras33), as method 1, the comparison experiment involving escaping with optimal avoidance strategy is conducted. The simulation yielded successful results, reaching the destination and M escaping against D1 49.366c and D2 1.067c. The final distance between M and T was 0.2376c, while M’s final velocity was 0.190, and it took 73.600s for M to reach T.

According to the method using lift coefficients as control quantities, as method 2, described in reference (Reference Liang, Li, Wu, Zheng, Chu and Wang6), the comparison experiment involving escaping with optimal avoidance strategy is conducted. M only succeeds in avoidance but not reaching T, and M escaping against D1 226.342c and D2 765.333c. The final distance between M and T was 30c.

According to the method based on a dynamic multi-objective algorithm, as method 3, described in reference (Reference Zheng-Ping, Hui-Cai, Zhi-Qiang, Kai-Feng and Ze-Xuan34), the comparison experiment involving escaping with optimal avoidance strategy is conducted. The simulation yielded successful results, reaching the destination and M escaping against D1 1.027c and D2 1.816c. The final distance between M and T was 0.333c, while M’s final velocity was 0.160, and it took 200.600s for M to reach T.

The trajectories are shown in Figs 4, 5 and 6, where M is the escapee’s trajectory and D1 and D2 are the pursuers’ trajectories. The dashed line shows the control experiment, contrast-M is the escapee’s trajectory, contrast-D1 and contrast-D2 are the pursuers’ trajectories.

Through the analysis depicted in Figs 4–6, it is evident that the optimal avoidance strategy delineated in this study surpasses the methodology expounded in the existing methods. This superiority is demonstrated by the heightened efficacy in enabling entity M to more adeptly evade the pursuers D1 and D2. Subsequent to the successful evasion, M systematically employs the optimal avoidance strategy to navigate towards the target point T. Using this method greatly reduces the final distance between M and T. It also cuts the task completion time. The final speed greatly rises with the best avoidance strategy. These results highlight the strategy’s high effectiveness and superiority.

Figures 7–10 serve as visual representations of the simulated outcomes arising from the evasion process for both the proposed optimal avoidance strategy and the control group experiment, denoted by the dashed line. The simulation results distinctly illustrate that the optimal avoidance strategy delineated in this study empowers entity M to execute manoeuvers of heightened versatility, characterised by a smaller angle-of-attack, able to efficiently complete the escape to reach T.

Figure 7. Overload variation curve.

Figure 8. Track angle variation curve.

Table 2. Intrusion process interference settings

Figure 9. Track yaw angle variation curve.

Figure 10. Angle-of-attack variation curve.

In order to test the effectiveness of the optimal avoidance strategy proposed in this paper, we conducted 1,000 anti-jamming simulation experiments [Reference Zhou Hongyu and Xiaogang35]. We establishes a Gaussian distribution of the interference factors and recorded the results in Table 2. During the experiments, we observes that M’s evasion against D1 ranged from 128.133c to 90.233c, while M’s evasion against D2 ranged from 2.699c to 1.911c. These distances are greater than the capture radius of the pursuer, indicating that M was able to evade successfully in all 1,000 experiments.

During the experiments, we observes that M’s evasion against D1 ranged from 128.133c to 90.233c, while M’s evasion against D2 ranged from 2.699c to 1.911c. These distances are greater than the capture radius of the pursuer, indicating that M was able to evade successfully in all 1,000 experiments.

The results of the experiment are depicted in Figs 11–13. These figures demonstrate that the non-cooperative differential game optimal avoidance strategy is robust enough to meet different constraints, even in the presence of disruptive factors.

Figure 11. Angle-of-attack variation curve.

Figure 12. Track angle variation curve.

Figure 13. Track yaw angle variation curve.

In comparison, methods 1 and 2 have poor robustness, while method 3 has a 30% probability of successful escape and has no comparative value, so no graphs are shown.

The simulation comparison indicates that for the ‘two pursuing one’ UAV pursuit-evasion scenario, the proposed method in this paper excels in evading pursuit and conserving energy. It also demonstrates strong robustness under various interferences. In three-dimensional simulations, the method shows significant potential for engineering applications.