Optimal avoidance strategy based on nonlinear approximate analytic solution of non-cooperative differential game

Abstract

This study examines the pursuit-evasion game involving unmanned aerial vehicles (UAVs), with a specific focus on the scenario of N-pursuers-one-escapee. The primary objective is to develop an optimal strategy for the escapee when the pursuers possess superior capabilities. To obtain this objective, we conduct the following study. Firstly, to enhance realism, a non-cooperative differential game model is formulated, incorporating multiple motion characteristics, including aerodynamics, overloading, and imposed constraints. Secondly, the end-value performance index is subsequently converted to an integral one, simplifying the solution process of the Hamilton-Jacobi-Bellman (HJB) equation. An iterative method is utilised to determine the covariates using the Cauchy initial value problem, and its convergence and uniqueness are established. The optimal avoidance strategy is subsequently derived from the covariates. Finally, the superiority of the proposed strategy is validated through simulation experiments and compared to three advanced optimal avoidance strategies. A total of 1,000 anti-jamming simulation experiments are conducted to verify the robustness of the proposed strategy.

Nomenclature

$x,y,z$

the coordinates in the North Celestial Inertial System

$V$

the magnitude of velocity

$D$

the aerodynamic drag

$L$

the aerodynamic lift

$m{\rm{\;}}$

the mass

$g$

the local gravitational acceleration

${D_i}$ , $\left( {i = 1,2} \right)$

the i-th pursuer

$M$

the escapee

${V_{Dmi}}$

the velocity projection of $M{D_i}$ the firing plane of M

${V_{Ddi}}$

the velocity projection of $M{D_i}$ the firing plane of D_i

${V_{Tmi}}$

the component of velocity of $M$ in the horizontal plane

${V_{Tdi}}$

the component of velocity of ${D_i}$ in the horizontal plane

${R_{Di}}$

the relative distance of the firing plane

${\bar R_{Di}}$

the transformed relative distance of the firing plane

${\bar R_{Ti}}$

the transformed relative distance of the horizontal plane

${x_i}$

the relative velocity

${f_{Di}}$

the $M{D_i}$ of the firing plane system $M$ control volume coefficient

${g_{Di}}$

the $M{D_i}$ of the firing plane system $D$ control volume coefficient

${h_{Di}}$

the $M{D_i}$ of the firing plane system coefficient

${f_{Ti}}$

the $M{D_i}$ of the horizontal plane system $M$ control volume coefficient

${g_{Ti}}$

the $M{D_i}$ of the horizontal plane system $D$ control volume coefficient

${h_{Ti}}$

the $M{D_i}$ of the horizontal plane system coefficient

${u_{Dm}}$

the firing plane control variable

${u_{Tm}}$

the horizontal plane control variable

$J$

the perferomance index

${t_{fi}}$

the terminal moment of $M{D_i}$

${t_0}$

the initial time

${l_i}$

the time operator

$H$

the Hamiltonian function

$p(t)$

the left term of the Cauchy’s initial value equation

$q(t)$

the right term of the Cauchy’s initial value equation

$Q$

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)$

$S$

the upper bound of $\left| {\lambda _i^{(n)}({x,t}) - \lambda _i^{({n - 1})}({x,t})} \right|$

$L$

A sufficiently large constant (used in the Section 3.3.2)

$N$

A sufficiently large constant (used in the Section 3.3.3)

$n$

number of iterations

Greek symbol

${\rm{\Theta }}$

the trajectory angle

${\rm{\Psi }}$

the trajectory yaw angle

$\nu $

the velocity pitch angle

$\alpha $

the angle-of-attack

${\lambda _{Ti}}$

the line-of-sight azimuth

${\lambda _{Di}}$

the line-of-sight elevation angle

${\phi _{Tmi}}$

the front angle of the horizontal plane of $M$

${\phi _{Dmi}}$

the front angle of the $M{D_i}$ firing plane of $M$

${\phi _{Tdi}}$

the front angle of the horizontal plane of ${D_i}$

${\phi _{Ddi}}$

the front angles of the $M{D_i}$ firing plane of ${D_i}$

${\gamma _{Tmi}}$

the horizontal plane velocity angle of $M$

${\gamma _{Dmi}}$

the firing plane velocity angle of $M$

${\gamma _{Tdi}}$

the horizontal plane velocity angle of ${D_i}$

${\gamma _{Ddi}}$

the firing plane velocity angle of ${D_i}$

${\alpha _{Di}}$

conversion factor

${\beta _i}$

the weight of the terminal value index

${\beta _{i0}}$

the initial weight of the terminal value index

${\kappa _T}$

the weight of the energy metric

$\phi $

Hamiltonian variation

${\lambda _i}$

the covariant variable

$\zeta $

the moment in $\left( {t,{t_0}} \right)$ in proof of convergence

$\xi $

the moment in $\left( {t,{t_0}} \right)$ in proof of uniqueness

1.0 Introduction

In a pursuit-evasion game, the escapee avoids the pursuer to reach a goal. Meanwhile, the pursuer aims to catch the escapee before it arrives. The escapee’s situation changes over time, demanding a dynamic avoidance strategy. Differential games are often applied here. Problems vary, based on the number of players on each side. For example, there are games with one pursuer and one escapee, or with N pursuers and one escapee, or with N pursuers and N escapees. In a game with multiple players, the one-pursuer-multi-escapees scenario is the subject of concern. An example is the ‘two cutters and a ship’ game. Here, two fast pursuers team up to catch a slower escapee quickly. This game has sparked much research. Differential games have applications in economics, sports, robotics and air combat [1].

This study focuses on the best strategy for UAVs in pursuit-evasion games. It is important to note that the applications of these mathematical tools may not be suitable for the scenario. However, the motion of the vehicle is non-linear and complex. Moreover, directly solving the joint multiple equations will lead to a dimensional disaster. Therefore, it highlights the growing in the optimal avoidance problem with nonlinear differential games for multiple vehicles in recent years [2, 3].

Existing studies often use linearisation and geometric methods to solve the optimal avoidance strategy of an escapee.

(1) Linearisation method: Solving the UAV’s non-linear dynamic equations with Linear Quadratic Differential Game (LQDG) makes it simple [4–9]. However, the variable zero effort miss (ZEM), an important factor, only works for slower, closer scenarios [10–13]. LQDG uses a transfer matrix to solve ZEM, leading to large errors. Additionally, simplifying the evasion problem with many pursuers into a two-player game and merging strategies overlooks multiple interceptors’ impact [14]. This, in turn, weakens the evasion effect.

(2) Geometric approach: The geometric approach uses Apollo circles [15], Tyson polygons [16] and the Voronoi diagram [17, 18] to define escapee and pursuer dominance zones. It also sets the boundary where they meet [19–22]. Geometric methods avoid nonlinear solution problems but assume uniform motion, a big limitation. Real-world relationships are complex and change. This flawlessly simplifies spatial relationships. However, it impacts the vehicle’s evasion effectiveness.

Scholars have solved differential game problems using various methods. These include the viscous solution [23–25] and weak Kolmogorov-Arnold-Moser theory (KAM) theory [26], etc. S. Hamadene [27] proposed the solution of the mixed zero and stochastic differential game by applying the oriented stochastic differential equations with local solutions. Shipei [28] established a nonlinear parabolic variational inequality for the HJB equation under sufficiently smooth conditions, proving that the value function is the unique viscosity solution of its inequality with applications. However, the above mathematical derivation is valid under limited conditions, such as the time periodicity of the viscous solution [29], the control variable tends to 0, and so on. Therefore, a solution method with high applicability and low error is needed for the flying vehicle game problem

This paper suggests an optimal strategy for the ‘N-pursuers-1-escapee’ game. It is adapted to UAV dynamics. To obtain this strategy, we conduct the following study. First, a differential game model is constructed from the UAV’s position-velocity relationship. Secondly, the final avoidance amount is taken as the performance index and converted to integral type to construct the Hamiltonian function to reduce the number of joint equations and avoid falling into dimensional catastrophe. Then, the covariates that satisfy the Cauchy initial value problem are constructed, and its convergence and uniqueness are proved. So the optimal avoidance strategy is obtained. Compared to three advanced optimal avoidance strategies, the simulation results show strong feasibility and robustness of the proposed strategy in this paper.

The contribution of this paper to the problem of pursuit-evasion game by UAVs can be summarised as follows:

1. (1) This method simplifies complex dynamics and finds a solution for vehicle games.

2. (2) It solves HJB system covariates. It works with the Cauchy problem, proving their convergence and uniqueness.

2.0 Problem description and modelling

2.1 Description of problems

This paper’s research objects are both the escapee and pursuers. It focuses especially on the escapee. Its task is to avoid pursuit and reach its destination.

The confrontation scenario is shown in Fig. 1: M is the escapee, D1 and D2 are pursuers and T is M’s destination.

Figure 1. Confrontation scenario.

Where ${\lambda _{D1}},{\lambda _{D2}}$ are the elevation angles of the line of sight, and ${\lambda _{T1}},{\lambda _{T2}}$ is the azimuth angle of the line of sight. Pursuers D1 and D2 pursue escapee M at different heights and directions, respectively. S1 and S2 are the target coordinate systems of D1 and D2.

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.0 Non-cooperative differential game optimal avoidance strategy design

This section solves the problem of avoiding two pursuers in the UAV pursuit-evasion scenario. The solution revolves around a function and key variables. The study introduces a new method for solving complex differential games. It directly uses three-degree-of-freedom equations. First, it sets up a performance index to manage the game dynamics. Then, it creates a Hamiltonian function. This function is derived from the initial index. It determines the best control variable with the covariate. The study then introduces a fresh approach. It integrates the Cauchy initial value problem. This integration helps to systematically derive solutions for the variables. The paper also proves the convergence and uniqueness of these solutions. This work improves the practical and theoretical strategies for UAV evasion.

The process is illustrated in Fig. 3.

Figure 3. Process.

Where $um$ is the optimal avoidance strategy for UAV.

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

(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 [31], 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 (32), 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.0 Simulation results and comparison

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 (33), 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 (6), 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 (34), 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 [35]. 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.

5.0 Conclusion

In response to the challenge posed by UAVs pursuit-evasion scenarios, this paper presents an optimal avoidance strategy based on non-cooperative differential game theory. Embracing the intricate dynamics inherent to UAVs, this strategy is tailored for practical multi-vehicle game scenarios, leveraging a three-degree-of-freedom nonlinear dynamic model. A pivotal aspect of this investigation entails the analytical resolution of nonlinear HJB equations. This analytical solution is attained via the Cauchy initial value problem, with its convergence and uniqueness meticulously proofed. Comparative analysis juxtaposes the simulation outcomes of the proposed optimal evasion strategy against those of three traditional advanced game strategies. The results unequivocally highlight the superior efficacy of the proposed strategy in evasion manoeuvers, accompanied by notable energy conservation. Subsequently, systematic anti-jamming simulation experiments are conducted to validate the viability and robustness of the proposed evasion strategy under diverse uncertainties. Demonstrating applicability and effectiveness, in scenarios where the pursuer possesses superior capabilities to the evading party, this strategy holds promise in addressing UAV pursuit-evasion dynamics. The insights gleaned from this study enrich the discourse surrounding the avoidance of moving obstacles and adept pursuers by UAVs, pertinent to dynamics of real-world UAV operations.