Stochastic Potential Game-Based Target Tracking and Encirclement Approach for Multiple Unmanned Aerial Vehicles System

Abstract

Utilizing fully distributed intelligent control algorithms has enabled the gradual adoption of the multiple unmanned aerial vehicles system for executing Target Tracking and Encirclement missions in industrial and civil applications. Restricted by the evasion behavior of the target, current studies focus on constructing zero-sum game settings, and existing strategy solvers that accommodate continuous state-action spaces have exhibited only modest performance. To tackle the challenges mentioned above, we devise a Stochastic Potential Game framework to model the mission scenario while considering the environment’s limited observability. Furthermore, a multi-agent reinforcement learning method is proposed to estimate the near Nash Equilibrium strategy in the above game scenario, which utilizes time-serial relative kinematic information and obstacle observation. In addition, considering collision avoidance and cooperative tracking, several techniques, such as novel reward functions and recurrent network structures, are presented to optimize the training process. The results of numerical simulations demonstrate that the proposed method exhibits superior search capability for Nash strategies. Moreover, through dynamic virtual experiments conducted with speed and attitude controllers, it has been shown that well-trained actors can effectively act as practical navigators for the real-time swarm control.

1. Introduction

The Multiple Unmanned Aerial Vehicles System (MUAVS) plays an important role in both industrial application and civil field with the rapid evolving speed of intelligent distributed control algorithms. Deploying MUAVS in disaster relief efforts [1] or infrastructure checking [2] for remote regions is beneficial not only for saving the expenses of manpower and material resources but also for improving the efficiency of mission completion. In these real-world applications, system must possess both obstacle avoidance and Target Tracking and Encirclement (TTE) capabilities [3] for non-cooperative objects.

TTE problem is a well-established topic in the field of aerial vehicle guidance and control, primarily focused on the estimation of target state and data association [4]. To achieve the goal, classic measures such as the Kalman Filter [5,6,7,8] and Symmetric Measurement Equations [9,10,11] are widely deployed with the automatic vehicle controllers that have practical value. In the controller part, some studies [12,13,14,15] pay attention to Model Predictive Control for achieving tracking. As the multi-agent system has gradually developed in the past decades, approaches for navigating collaborative UAV swarms have been studied. Several methods developed from graph theory are mentioned in [16,17], revealing the ideal results in kinematic simulation. In [18,19,20], the Actor–Critic (AC)-based MARL frameworks with different constraints are proposed to solve single- or multiple-target tracking, trajectory optimization, and collision avoidance problems, whose main innovations focus on unique reward functions. When the number of tracking members is insufficient, a novel Voronoi diagram with the Nash Q-learning method is introduced in [21] to produce UAV control input. Many of these solutions provide optimal results in simulation but are not suitable for taking account of the evasion tendencies of the target in real-world scenarios.

The field of game theory is experiencing rapid growth due to the rising interest in model-free multi-agent reinforcement learning (MARL) techniques [22,23,24,25]. Recent MARL studies focus on improving how agents effectively utilize observation information under partially observable conditions to enhance the efficiency of collective reward acquisition. In [26], the graph-based MARL algorithm is proposed by utilizing directed graph between agents and multi-heads attention mechanism to guide Q network iteration. After that, the branch of research direction based on directed graph gradually shows the top performance of multi-agents application. In [27], a novel information exchange mechanism between different agents is employed to update the critic based on the AC framework. This paper considers the interdependencies among agents and achieves better performance in standard test scenarios. In order to personalize peer-to-peer communication topology during task execution, a self-adapting graph formation is presented in [28]. This method does not restrict the communication topology between agents and can be dynamically adjusted according to the task execution conditions.

On the other hand, the advancement of modern computing technologies has led to significant progress in the numerical solution of general-sum games. A general-sum game is a game where the sum of all players’ payoffs is not necessarily zero, meaning players’ interests can be both competitive and cooperative. After the Nash Q-learning [29,30] created a precedent that RL methods could deal with game problems effectively, many model-free MARL methods as well as NE strategy solvers have been proposed in recent years. In [23], a general MARL algorithm called Policy-Space Response Oracles is introduced to obtain effective meta-strategies in partial observable environments, which attains the highest score in the Poker Game benchmark. The Nash Equilibrium (NE) [31] convergence of AC approaches is proven in [32] by Two-Timescales stochastic approximation, that is, breaking the hamper existing between continuous general-sum game and popular RL algorithms. Meanwhile, there have been many efforts in previous studies to apply game theories to the TTE problem. The Hamilton–Jacobi Equation for the Pursuer–Evader differential game based on the zero-sum formulation is discussed in [33], which provides a new perspective for solving general class target-hitting games. Furthermore, the bounded rational level-k Hamilton–Jacobi–Issca Equation in [34] is used to express the Pursuer–Evader game scenario for different levels of intelligence. While zero-sum games provide a straightforward description of the interaction between the trackers and target, other crucial factors such as collision avoidance and controllability are not accounted for.

Among many environment modeling methods is a game formulation known as the Potential Game (PG) [35]. Agents who have the unified potential function in PG are not acting like a team or opponents strictly but operating for multiple goals simultaneously and distributedly. Such characteristics are suitable for describing the aforementioned TTE scenario of MUAVS that needs to survive in the complex environment. In this condition, system members play as self-interesting agents trying to maximize their rewards affected by relative kinematic state and unknown environment in continuous action and state space. As a traditional problem of PG field, searching the NE solution is a common challenge for all related research [22,32,36]. The NE of a game represents a joint mixed strategy profile in which any agent cannot obtain higher payoff by modifying their own strategy profile. To solve this Polynomial Parity Arguments on Directed (PPAD) graphs difficulty issue [37], recent research has proposed the Stochastic Potential Game (SPG) concept, combining the Markov Decision Process (MDP) with PG to simplify the difficulty of solving [38,39], which adopts MARL approaches based on the Bellman operator as the solver. Although these methods could easily create state-of-the-art performance in discrete system with multiple agents [40], the effectiveness in continuous action-state space, such as MUAVS deployed in the real world, remains to be improved. Meanwhile, most research typically focuses on kinematic simulation without considering the dynamic properties of multi-agent systems; achieving intelligent auto-control in such systems remains an ongoing challenge.

Regarding the TTE problem, when targets possess intelligence comparable to that of the trackers, conventional target tracking algorithms face significant challenges in effectively coordinating UAV swarms for collaborative tracking missions. In real-world scenarios, it is challenging for UAVs to obtain the target’s path or predict its trajectory in advance, and the target may employ evasive maneuvers or countermeasures. Under such circumstances, it becomes essential for the UAV trackers to develop capabilities to address these complex situations effectively. Therefore, in this paper, we address the above-mentioned challenge for multi-agent systems operating in an unknown environment based on the SPG theory. While selecting quadrotors as system agents, a MARL approach is proposed for approximating the NE of the continuous SPG, which considers kinematic and dynamic constraints of both the trackers and target. The primary contribution of this paper is outlined below:

• The SPG-based TTE scenario is constructed within a finite continuous action-state domain that involves adjusting the composite tasks of an MUAVS in various environments to achieve the TTE objective. All trackers are dealing with the uncertainty of the target’s motion pattern, as the target possesses equivalent intelligent decision-making capabilities within the scenario.
• To search for local NE strategies, a Time-Series Multi-Agent Soft Actor–Critic (TMSAC) approach is proposed. It leverages sequential observations from agents and is particularly effective in determining the optimal strategy for MUAVS, especially when involving multiple agents. Furthermore, novel reward functions that account for the SPG condition are designed. These functions are integrated with the TMSAC framework to enhance performance for both the tracker and the target. The convergence of the algorithm is also discussed to justify its search capability for the local NE.
• Considering the dynamical characteristics of agents, a guidance loop based on the actor training from TMSAC is combined with velocity and attitude controllers, deployed in the visualized physical simulation environment to show the effectiveness and the success rate of the proposed method.

The remainder of this article is organized as follows. Before analyzing environment properties, Section 2 presents the dynamic and relative kinematic models used in further research. After that, the construction of SPG and MARL method details are introduced in Section 3 and Section 4. Finally, multiple experiments and visualized simulation for performance tests are discussed in Section 5.

2. System Modeling

In this section, we analyze the relevant quadrotor dynamics and multi-agent kinematics before introducing the game environment and tracking formulation, performing as the fundamental background for further research.

2.1. Quadrotor Dynamics

We select a small-size quadrotor to play as the trackers and target in our research. As a significant factor for the tracking navigator design, it is necessary to notice the following assumptions before constructing a dynamic model.

Assumption 1. 

The UAV model can be regarded as a rigid body, and the Little Disturbance Approach is used to linearize the dynamic formulation on the equilibrium points.

Assumption 2. 

The mass distribution of the quadrotor is symmetrically arranged along all axes of the Body Reference Frame.

As in Figure 1, we present a universal quadrotor model that employs two coordinate systems, the Body Reference Frame (BRF) and the Earth Reference Frame (ERF), which conform to the Right-hand Rule. The BRF origin point $O_b$ is equal to the aircraft barycenter, while the $O_b{Z}_b$ axis is perpendicular to the fuselage plane, with the positive direction orienting to ground. The $O_b{X}_b$ and $O_b{Y}_b$ are perpendicular to each other in the fuselage plane. The ERF origin point $O_g$ is fixed on the ground plane, with the $O_g{X}_g$ heading to geographic north and the $O_g{Y}_g$ orienting to geographic east. The $O_g{Z}_g$ orients to the geocenter.

The velocity $\mathit{\nu }={\left[u,v,w\right]}^{\mathrm{T}}$ and angular velocity $\mathit{\omega }={\left[p,q,r\right]}^{\mathrm{T}}$ are defined based on the BRF. Meanwhile, coordinates vector $\mathit{P}={\left[x,y,z\right]}^{\mathrm{T}}$ and Euler angle $\mathit{\Theta }={\left[\varphi ,\theta ,\psi \right]}^{\mathrm{T}}$ are used to describe the quadrotor position and posture in the ERF. Thus, the transformation between the BRF and the ERF can be described as Equation (1) [41]:

$$\left[\begin{array}{c}\dot{\mathit{P}}\\ \dot{\mathit{\Theta }}\end{array}\right]=\left[\begin{array}{cc}{\mathit{R}}_b^e& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathit{R}}_E\end{array}\right]\left[\begin{array}{c}\mathit{\nu }\\ \mathit{\omega }\end{array}\right],$$

(1)

where ${\mathit{R}}_b^e$ is the coordinate transfer matrix defined as Equation (2), and ${\mathbf{0}}_{3\times 3}$ is a 2-D $3\times 3$ zero matrix. For concision, $sin(·)$ and $cos(·)$ are represented by s and c separately,

$${\mathit{R}}_b^e=\left[\begin{array}{ccc}c\theta ·c\psi & s\varphi ·s\theta ·c\psi -c\varphi ·s\psi & c\varphi ·s\theta ·c\psi +s\theta ·s\psi \\ c\theta ·s\psi & s\varphi ·s\theta ·s\psi +c\varphi ·c\psi & c\varphi ·s\theta ·s\psi -s\theta ·c\psi \\ -s\theta & s\varphi ·c\theta & c\varphi ·c\theta \end{array}\right].$$

(2)

The ${\mathit{R}}_E$ is the Euler rotation matrix defined as

$${\mathit{R}}_E=\left[\begin{array}{ccc}1& sin\varphi ·tan\theta & cos\varphi ·tan\theta \\ 0& cos\varphi & -sin\varphi \\ 0& sin\varphi ·sec\theta & cos\varphi ·sec\theta \end{array}\right].$$

(3)

The mass of quadrotor is presented as m. The principal axes of inertia are completely aligned with the axes of the BRF. According to Assumption 2, the product of inertia $\left[I_{xy},I_{xz},I_{yz}\right]=0$, and moment of inertia are shown as $\left[I_x,I_y,I_z\right]$. Since the output force vector from propellers $\mathit{F}=F_1+F_2+F_3+F_4$ and the moment $\mathit{M}$ made by $\mathit{F}$ are applied in the BRF, the dynamic model of quadrotor is given by Equations (4) and (5) [42]:

$$\dot{\mathit{\nu }}=\left[\begin{array}{c}rv-qw\\ pw-ru\\ qu-pv\end{array}\right]+{\displaystyle \frac{1}{m}}\left[\begin{array}{c}{\mathit{F}}_x\\ {\mathit{F}}_y\\ mg-{\mathit{F}}_z\end{array}\right],$$

(4)

$$\dot{\mathit{\omega }}=\left[\begin{array}{c}\left[\left(I_y-I_z\right)qr+{\mathit{M}}_x\right]/I_x\\ \left[\left(I_z-I_x\right)pr+{\mathit{M}}_y\right]/I_y\\ \left[\left(I_x-I_y\right)pq+{\mathit{M}}_z\right]/I_z\end{array}\right].$$

(5)

2.2. Multi-Agent Kinematics

Multi-agent kinematics represent the relative motion of the MUAVS members and the target, performing as the basis of SPG design. Taking the dynamic properties of quadrotors into account, it is essential to propose the following assumptions before defining the kinematics, which models the game environment reasonably.

Assumption 3. 

In the multi-agent kinematics modeling, both MUAVS and the target are moving in the same two-dimensional plane, and no motion in the ERF Z-axis needs to be considered.

Assumption 4. 

All quadrotors can obtain the relative position information of other agents through the airborne sensors.

Assumption 5. 

Giving by Assumption 3, the ground velocity $v^i$ synthesized by the $\dot{x}$ and $\dot{y}$ of the trackers are restricted in $v_{max}\ge v^i$. The yaw rate $\dot{{\psi }^i}$ satisfies ${\dot{\psi }}_{max}\ge \left|\dot{{\psi }^i}\right|$.

For an MUAVS that has N units, each quadrotor $\mathbf{i}$ in the system satisfies the following formula to indicate the motion state:

$$\begin{array}{cc}{\dot{\kappa }}_u^i& ={\left[\begin{array}{ccc}{v}^{i}sin{\psi }^i& v^{i}cos{\psi }^i& \dot{{\psi }^i}\end{array}\right]}^{\mathrm{T}},\end{array}$$

(6)

where $[x_i,y_i]$ is the position of quadrotor i in the ERF. In the meantime, the target has the same form of the kinematic:

$$\begin{array}{cc}{\dot{\kappa }}_t& ={\left[\begin{array}{ccc}{v}^{t}sin{\psi }^t& v^{t}cos{\psi }^t& \dot{{\psi }^t}\end{array}\right]}^{\mathrm{T}}.\end{array}$$

(7)

Therefore, the relative position based on the Polar Coordinates for tracker–target and tracker–tracker can be given by

$${\kappa }_{u-t}^i={\left[\begin{array}{cc}{d}_t^i& {\phi }_{u-t}^i\end{array}\right]}^{\mathrm{T}},$$

(8)

$$\begin{array}{ccc}\hfill {\kappa }_{u-u}^i=& [d_u^{i,1}\cdots d_u^{i,i-1}\hfill & \hfill d_u^{i,i+1}\cdots d_u^{i,N},\\ \hfill & {\vartheta }_u^{i,1}\cdots {\vartheta }_u^{i,i-1}\hfill & \hfill {\vartheta }_u^{i,i+1}\cdots {\vartheta }_u^{i,N}& {]}^{\mathrm{T}}\hfill \end{array},$$

(9)

where the $d_t^i=\sqrt{{(x_u^i-x_t)}^2+{(y_u^i-y_t)}^2}$ represents the distance between tracker $\mathbf{i}$ and the target. In particular, ${\phi }_{u-t}^i$ is the relative yaw angle indicated as ${\phi }_{u-t}^i={\psi }_i-{\xi }_t^i$, where ${\xi }_t^i=arctan{\displaystyle \frac{x_t-x^i}{y_t-y^i}}$. The $d_u^{i,j}$ is the distance between tracker i and j, while the ${\vartheta }_u^{i,j}=|{\xi }_t^i-{\xi }_t^j|$ is the relative encirclement angle. Based on the aforementioned definitions, the system’s relative kinematics can be expressed as follows:

$${\dot{\kappa }}_{\mathrm{u}-\mathrm{t}}^i=\left[\begin{array}{c}{\dot{d}}_t^i\\ {\dot{\phi }}_{u-t}^i\end{array}\right]=\left[\begin{array}{c}{v}_u^{i}cos{\phi }_{u-t}^i-v_{t}cos({\psi }_t-{\xi }_t^i)\\ \dot{{\psi }^i}-\dot{{\xi }_t^i}\end{array}\right],$$

(10)

$${\dot{\kappa }}_{\mathrm{u}-\mathrm{u}}^i=\left[\begin{array}{c}{\dot{d}}_u^{i,1}\\ ⋮\\ {\dot{d}}_u^{i,i-1}\\ {\dot{d}}_u^{i,i+1}\\ ⋮\\ {\dot{d}}_u^{i,N}\\ {\dot{\vartheta }}_u^{i,1}\\ ⋮\\ {\dot{\vartheta }}_u^{i,i-1}\\ {\dot{\vartheta }}_u^{i,i+1}\\ ⋮\\ {\dot{\vartheta }}_u^{i,N}\end{array}\right]=\left[\begin{array}{c}{v}_{i}cos{\phi }_{u-u}^{i,1}-v_{1}cos\left({\psi }^1-{\xi }_u^{i,1}\right)\\ ⋮\\ v_{i}cos{\phi }_{u-u}^{i,i-1}-v_{i-1}cos\left({\psi }^{i-1}-{\xi }_u^{i,i-1}\right)\\ v_{i}cos{\phi }_{u-u}^{i,i+1}-v_{i+1}cos\left({\psi }^{i+1}-{\xi }_u^{i,i+1}\right)\\ ⋮\\ v_{i}cos{\phi }_{u-u}^{i,N}-v_{N}cos\left({\psi }^N-{\xi }_u^{i,N}\right)\\ \nabla arccos\left(\left({d_t^1}^2+{d_t^i}^2-{d_u^{1,i}}^2\right)/2d_t^1{d}_t^i\right)\\ ⋮\\ \nabla arccos\left(\left({d_t^i}^2+{d_t^{i-1}}^2-{d_u^{i,i-1}}^2\right)/2d_t^i{d}_t^{i-1}\right)\\ \nabla arccos\left(\left({d_t^i}^2+{d_t^{i+1}}^2-{d_u^{i,i+1}}^2\right)/2d_t^i{d}_t^{i+1}\right)\\ ⋮\\ \nabla arccos\left(\left({d_t^i}^2+{d_t^N}^2-{d_u^{i,N}}^2\right)/2d_t^i{d}_t^N\right)\end{array}\right].$$

(11)

3. Stochastic Potential Game

In this section, the SPG formulation is used to establish the game relations between MUAVS and the target, simultaneously considering objectives such as tracking, encirclement, and obstacle avoidance. In the practical scenario where MUAVS executes the tracking and encirclement mission, the system performance is indicated by not only relative kinematics but also the survivability in an unknown environment. It means that the traditional zero-sum Pursuer–Evader game is not suitable for formulating the above situation. Motivated by such a reason, we adopt the SPG theory derived from the general-sum PG and MDP to describe the game status of the system.

The SPG environment is shown in Figure 2. We represent the continuous SPG as the tuple $\mathcal{G}:=\left\{\mathcal{N},\mathcal{S},\mathbf{\Pi },\mathcal{R},\mathcal{P}\right\}$, where $\mathcal{N}:=\{1,2,\dots ,N,t\}$ denotes the set of agents that has both the trackers and target. The remaining elements of the tuple $\mathcal{G}$ are described as follows:

• $\mathcal{S}:=\{s^1,s^2,\dots ,s^N,s^t\}\forall i\in \mathcal{N}$ indicates the set of agents state constructed by relative kinematic information and obstacle observation.
• $\mathbf{\Pi }:={\times }_i^N{\pi }^i,\forall {\pi }^i\in {\mathcal{A}}^i$ is the joint Markov strategy set where ${\pi }^i$ is the mix strategy gathered from continuous action space ${\mathcal{A}}^i$.
• $\mathcal{P}:=\{P^1,P^2,\dots ,P^N,P^t\}\forall P^i:s^i\times {\pi }^i\times s^i\in [0,1]$ is the state transfer probability function expressing the system response characteristic.
• $\mathcal{R}:=\{R^1,R^2,\dots ,R^N,R^t\}\forall i\in \mathcal{N}$ is the rewards function for each agent in $\mathcal{N}$.

In the traditional research field of stochastic game (SG), an augmented MDP involves agents from $\mathcal{N}$ who seek to maximize their rewards. The SPG concept is a subset of SG that combines PG properties, where a global potential function quantifies agents’ behaviors without relying on the characteristic of individual agents. In each step, every agent i takes an action $a^i\in {\pi }^i$, leading to a change in the potential function, while state $\mathcal{S}$ transfers to a new point based on $\mathcal{P}$. The essential task in constructing the SPG is to demonstrate the existence of a potential function that has a unified formulation to represent the rewards change of agents. In the case we focus on, there is a global potential function $\varphi :\mathcal{S}\times \mathbf{\Pi }\to \mathbb{R}$ that for any agent should satisfy Definition 1.

Definition 1. 

An SG is an SPG if a global potential function can describe the differential of all agents’ rewards (cost matrices). For any $i\in \mathcal{N}$,

$$R^i\left[s,(a^i,a^{-i})\right]-R^i\left[\widehat{s},({\widehat{a}}^i,a^{-i})\right]=\varphi \left[s,(a^i,a^{-i})\right]-\varphi \left[\widehat{s},({\widehat{a}}^i,a^{-i})\right],$$

(12)

where $a^{-i}:=\left\{a^1,a^2,\dots ,a^{i-1},a^{i+1},\dots ,a^N\right\},\forall i\in \mathcal{N}$ is the joint action except $a^i$. ${\widehat{a}}^i$ is the next step of agent i, and ${\widehat{s}}^i$ represents the new state caused by ${\widehat{a}}^i$.

If all agents have the same form of reward function, then Definition 1 is naturally satisfied. However, in circumstances where reward functions are not specific or not in a unified format, finding a potential function in a definite form becomes difficult. Therefore, demonstrating the existence of a potential function is an alternative solution to prove that a SG satisfies the conditions of an SPG. Based on Definition 1, the variant equation could be established as

$${\displaystyle \frac{\partial R^i\left[s,(a^i,a^{-i})\right]}{\partial a^i}}={\displaystyle \frac{\partial \varphi \left[s,(a^i,a^{-i})\right]}{\partial a^i}}.$$

(13)

Given the condition stated in Equation (13), Lemma 1 provides another way to prove that the described scenario satisfies the definition of an SPG. In Section 4, it serves as a precondition for deploying MARL methods to search for the NE strategy.

Lemma 1. 

For any player i and j in the game environment, the game $\mathcal{G}$ is a continuous SPG if and only if their reward functions (cost matrices) satisfy

$${\displaystyle \frac{{\partial }^2{R}^i\left[s,(a^i,a^{-i})\right]}{\partial a^i\partial a^j}}={\displaystyle \frac{{\partial }^2{R}^j\left[s,(a^j,a^{-j})\right]}{\partial a^i\partial a^j}}.$$

(14)

In the case of MUAVS versus the target, the definition of reward functions is critical for building the game environment. To describe $R^i$, the following equation is defined to model the payoff during the training process of the algorithm:

$$R^i\left[s,(a^i,a^{-i})\right]=R_s^i(s,a^i)+R_c^i\left[s,(a^i,a^{-i})\right],$$

(15)

where we have the following:

• $R_s^i=R_{\mathrm{col}}^i\left(s^i\right)+R_{\mathrm{act}}^i\left(a^i\right)$ is called self-influence rewards, affected by the agent i obstacle avoidance and the smoothness of the action.
• $R_c^i=R_{\mathrm{trac}}^i(s^i,a^i,a^{-i})+R_{\mathrm{sur}}^i(s^i,a^i,a^{-i})$ represents the inter-influence rewards that include the tracking and encirclement state.

Specifically, depending on whether agent i is a tracker or the target, $R_{\mathrm{trac}}^i$ and $R_{\mathrm{sur}}^i$ have different expressions as shown in Equations (16) and (17):

$$R_{\mathrm{trac}}^i=\left\{\begin{array}{cc}{\displaystyle \frac{1}{N}}{f}_{trac}\left(d_t^i,a^i-a^j\right),\hfill & if i is a Tracker \hfill \\ -{\displaystyle \frac{1}{N}}\sum _{j=0}^N{f}_{trac}\left(d_t^j,a^j-a^i\right),\hfill & if i is the Target\hfill \end{array}\right.,$$

(16)

$$R_{\mathrm{sur}}^i=\left\{\begin{array}{cc}{f}_{sur}\left({\vartheta }_u^{i,-i},a^i,a^{-i}\right),\hfill & if i is a Tracker (Target is not in -i)\hfill \\ 0,\hfill & if i is the Target\hfill \end{array}\right..$$

(17)

Note that both $f_{trac}(·)$ and $f_{sur}(·)$ are defined as continuously differentiable functions for all $a^i\in {\pi }^i$. Additionally, $f_{trac}(·)$ is linear with respect to $a^i-a^j$. Detailed descriptions of these functions are provided in Section 4.2. The term $d_t^i$ represents the distance between tracker i (or j) and target j (or i), which is one element of the state vector s.

Based on the fact that all agents in $\mathcal{N}$ have the same form of rewards function shown in Equation (15), it is feasible to demonstrate the existence of the potential function using Lemma 1. There are two scenarios for Equation (16): either both agents i and j are members of the MUAVS, or one of them is the target. Given that $R_c^i$ is not affected by $a_j$, its partial derivative with respect to $a_j$ is zero as shown in Equation (18). Consequently, the partial derivative of $R^i$ can be transformed as in Equation (19):

$${\displaystyle \frac{{\partial }^2{R}_s^i\left(s,a^i\right)}{\partial a^i\partial a^j}}=0,$$

(18)

$${\displaystyle \frac{{\partial }^2{R}^i\left[s,(a^i,a^{-i})\right]}{\partial a^i\partial a^j}}=\left\{\begin{array}{cc}-{\displaystyle \frac{1}{N}}{\displaystyle \frac{{\partial }^2{f}_{trac}\left(d_t^i,a^i-a^j\right)}{\partial a^i\partial a^j}},\hfill & if j is the Target\hfill \\ {\displaystyle \frac{{\partial }^2{f}_{sur}\left({\vartheta }_u^{i,j},a^i,a^j\right)}{\partial a^i\partial a^j}},\hfill & if i and j are Trackers\hfill \end{array}\right..$$

(19)

When agent j represents the target, according to Equation (20), the partial derivative of $R_j$ can be expressed as shown in Equation (21):

$${\displaystyle \frac{{\partial }^2{f}_{trac}\left(d_t^k,a^k-a^j\right)}{\partial a^i\partial a^j}}=0,\forall k\in \mathcal{N} and k\ne i,$$

(20)

$${\displaystyle \frac{{\partial }^2{R}^j\left[s,(a^j,a^{-j})\right]}{\partial a^i\partial a^j}}=\left\{\begin{array}{cc}{\displaystyle \frac{1}{N}}{\displaystyle \frac{{\partial }^2{f}_{trac}\left(d_t^i,a^i-a^j\right)}{\partial a^i\partial a^j}},\hfill & if j is the Target\hfill \\ {\displaystyle \frac{{\partial }^2{f}_{sur}\left({\vartheta }_u^{j,i},a^j,a^i\right)}{\partial a^i\partial a^j}},\hfill & if i and j are Trackers\hfill \end{array}\right..$$

(21)

There are several points that should be noted:

• Due to the relativity of the Polar Coordinates used in Section 2, ${\vartheta }_u^{i,j}={\vartheta }_u^{j,i}$ when both agent i and j are members of MUAVS.
• According to the linearity of $f_{trac}(·)$, Equation (22) holds when either i or j represents the target in the pair $(i,j)$:

$${\displaystyle \frac{{\partial }^2{f}_{trac}\left(d_t^i,a^i-a^j\right)}{\partial a^i\partial a^j}}=0.$$

(22)

Based on the preceding discussion, both of the aforementioned scenarios satisfy the conditions stipulated in Lemma 1. Therefore, the existence of the potential function of the proposed cost metrics is proven.

Since [43] has proposed that efficiently searching for the NE in SG is PPAD-hard, some researchers have focused on particular scenarios to simplify the difficulty of finding the NE. In [38], it has been proven that solving an MDP can lead to finding the local NE of the SPG. The Markov Perfect Equilibrium (MPE), as presented in Definition 2, is closely related to the local NE in SPG. As an NE variant, it is an appropriate concept to link the MDP with game theories. Therefore, it is feasible to use stochastic dynamic optimization techniques like MARL to search for the MPE strategy in SPG.

Definition 2. 

An action strategy $\overline{\mathit{\pi }}:=\left\{{\overline{\pi }}^1,{\overline{\pi }}^2,\dots ,{\overline{\pi }}^N\right\}\in \mathbf{\Pi }$ is an MPE of game $\mathcal{G}$ if:

$$V_{\overline{\mathit{\pi }}}^i\left(s\right)\ge V_{\mathit{\pi }}^i\left(s\right),\forall \mathit{\pi }\in \mathbf{\Pi },\forall s\in \mathit{S},$$

(23)

where $V^i$ is the value function of MDP.

4. Time-Series Multi-Agent Soft Actor–Critic Approach

In this section, the TMSAC approach is introduced in detail to solve the SPG formulation constructed in Section 3. Since the Multi-Agent Deep Deterministic Policy Gradient (MADDPG) [24] became the common baseline of deep MARL algorithms, a large number of practical methods, such as Counterfactual Multi-Agent (COMA) [44] and Factored Multi-Agent Centralized (FACMAC) policy gradients [45], are consistently refreshing the ranking of the benchmark test results. However, designing MARL algorithms that are suitable for continuous action-state spaces and are robust in real-time control tasks is challenging and progressing slowly. To enhance the efficacy of algorithms, it is reasonable to build a framework that leverages the sequential nature of agents state information in continuous MARL methods, given the time continuity of tasks. Therefore, combined with the centralized training and distributed execution strategy indicated by MADDPG, we propose a new construction of AC that adopts the Gated Recurrent Unit (GRU) structure to handle continuously sequential state information.

4.1. State and Action Embedding

4.1.1. State Vector

To realize the end-to-end manipulation, the state vector $s^i$ of agent i is defined to incorporate observations from airborne sensors, including sequential relative kinematics information and obstacle observations. The state vector can be represented as follows:

$$s^i=\left\{\begin{array}{cc}\left\{\begin{array}{cc}\left[{\kappa }_{\mathrm{u}-\mathrm{t}}^1,\dots ,{\kappa }_{\mathrm{u}-\mathrm{t}}^j,\dots ,{\kappa }_{\mathrm{u}-\mathrm{t}}^N\right]\times L,& O_{obs}^i\end{array}\right\},\hfill & if i is the Target\hfill \\ \left\{\begin{array}{cc}\left[{\kappa }_{\mathrm{u}-\mathrm{t}}^i,{\kappa }_{\mathrm{u}-\mathrm{u}}^i\right]\times L,& O_{obs}^i\end{array}\right\},\hfill & if i is a Tracker\hfill \end{array}\right..$$

(24)

For the target, $\left[{\kappa }_{\mathrm{u}-\mathrm{t}}^1,\dots ,{\kappa }_{\mathrm{u}-\mathrm{t}}^N\right]$ represents the tracker–target relative kinematics as shown in Equation (10). Meanwhile, $\left[{\kappa }_{\mathrm{u}-\mathrm{t}}^i,{\kappa }_{\mathrm{u}-\mathrm{u}}^i\right]$ represents both the relative kinematics between tracker i and the target, and the relative kinematics among trackers within the group. As shown in Figure 3, the obstacle observation $O_{obs}^i$ can be obtained from the laser radar ranging system that is simplified to a $N_{obs}$-dimension distances vector isometrically sampled from the detection circle centered on the agent. Obtaining the state vectors during the tracking process is facilitated by the advanced capabilities of modern airborne sensors. The variable L denotes the observation duration that specifies the duration for which the agents will retain and utilize the historical information.

4.1.2. Action Vector

Considering that the action outputs are inputs of the UAV speed and attitude controller, the ground velocity $v^i$ and yaw rate $\dot{{\psi }^i}$ have to be manipulated simultaneously, performing as a practical UAV navigation system. Given by the kinematics, the actions of agents $a^i$ are designed in two dimensions:

$$a^i=f_{{\pi }^i}\left(s^i\right)={\left[\begin{array}{cc}{v}^i& {\dot{\psi }}^i\end{array}\right]}^{\mathrm{T}},$$

(25)

where $v^i$ and $\dot{{\psi }^i}$ are limited by Assumption 5.

4.2. Composite Rewards Function

In the SPG problem, a key challenge is designing the instructive reward functions that need to be optimized by the MARL algorithm. As shown in Equation (15), the total rewards of each agent are separated into two parts, self-influence rewards, and inter-influence rewards, whose essence is to measure the individual performance and the contribution of trackers simultaneously.

4.2.1. Self-Influence Rewards

Self-influence rewards serve as crucial indicators of an agent’s survival capability and comprise two key components. The first component focuses on control system stability—to ensure reliable and stable operation of the control system, the action sequences output by the algorithm must maintain smoothness and continuity while avoiding dramatic fluctuations. The second component concentrates on evaluating collision avoidance capabilities. In real-world scenarios, agents may encounter various types of obstacles, such as static buildings or dynamic obstacles. Based on these considerations, we decompose self-influence rewards into two specific reward functions: action smoothness reward $R_{act}$ and collision avoidance reward $R_{\mathrm{col}}$. These two rewards are quantitatively calculated through the following mathematical expressions to guide agents in learning safer and more stable behavioral strategies:

$$R_{\mathrm{act}}^i\left(a^i\right)=-{\beta }_{\mathrm{act}}\sum _{t=0}^{N_{\mathrm{act}}}\left|{\displaystyle \frac{\mathrm{d}{a}^i}{\mathrm{d}t}}\right|,$$

(26)

$$R_{\mathrm{col}}^i\left(s^i\right)=R_{\mathrm{obs}}^i\left(O_{\mathrm{obs}}^i\right)+R_{\mathrm{int}}^i\left({\kappa }_{\mathrm{u}-\mathrm{t}}^N,{\kappa }_{\mathrm{u}-\mathrm{t}}^i\right).$$

(27)

where ${\beta }_{\mathrm{act}}$ is the gain factor of the action change rate, and $N_{\mathrm{act}}$ means the size of the window used to sum. As indicated in Equation (26), collision avoidance rewards consist of obstacle-influence part $R_{\mathrm{obs}}^i$ and agent-influence part $R_{\mathrm{int}}^i$, described as follows:

$$R_{\mathrm{obs}}^i=\left\{\begin{array}{cc}{\beta }_{\mathrm{obs}}\sum _{n=0}^{N_{\mathrm{obs}}}\left(l_{\mathrm{n}}-l_{obs}\right),\hfill & if \forall l_{obs}\ge l_n>0\hfill \\ -{\beta }_{\mathrm{obs}}{N}_{\mathrm{obs}}{l}_{\mathrm{obs}},\hfill & if \exists l_n=0\hfill \end{array}\right.,$$

(28)

$$R_{\mathrm{int}}^i=\left\{\begin{array}{cc}{\beta }_{\mathrm{int}}\sum _{j=0}^{N_{\mathrm{int}}}\left(d^j-2l_{\mathrm{int}}\right),\hfill & if \forall 2l_{\mathrm{int}}\ge d^j>0\hfill \\ -2{\beta }_{\mathrm{int}}{N}_{\mathrm{int}}{l}_{\mathrm{int}},\hfill & if \exists d^j=0\hfill \\ 0,\hfill & if \forall d^j>2l_{\mathrm{int}}\hfill \end{array}\right..$$

(29)

where ${\beta }_{\mathrm{obs}}$ and ${\beta }_{\mathrm{int}}$ are respective gain factors. For $R_{\mathrm{obs}}^i$, $l_{\mathrm{n}}$ indicates the detected distance of each unit in $N_{\mathrm{obs}}$ of agent i. When any $l_{\mathrm{n}}=0$, it means a collision has occurred, and in this case, $R_{\mathrm{obs}}$ should be set to a negative value to represent a penalty. Meanwhile, in $R_{\mathrm{int}}^i$, $d^j$ is the abbreviation of $d_u^{i,j}$ and $d_t^i$, and $N_{\mathrm{int}}$ is the collection of agents whose $2l_{\mathrm{int}}\ge d^j$.

4.2.2. Inter-Influence Rewards

In this study, inter-influence rewards consist of two components, tracking utility $R_{\mathrm{trac}}^i$ and encirclement utility $R_{\mathrm{sur}}^i$, which jointly measure the performance of the mission. In detail, the ${\dot{d}}_t^i$ plays the main character in $R_{\mathrm{trac}}^i$ due to its expression of the kinematic changing rate. And ${\vartheta }_u^{i,j}$ is counted into $R_{\mathrm{sur}}^i$ to evaluate the encirclement status. Therefore, the aforementioned reward functions are denoted as

$$R_{\mathrm{trac}}^i=\left\{\begin{array}{cc}{\displaystyle \frac{{\beta }_{\mathrm{trac}}}{N}}\left(v_{\max }^t+v_{\max }-{\dot{d}}_t^i\right),\hfill & if i is a Tracker\hfill \\ -{\displaystyle \frac{{\beta }_{\mathrm{trac}}}{N}}\sum _{j=0}^N\left(v_{\max }^t+v_{\max }-{\dot{d}}_t^i\right),\hfill & if i is the Target\hfill \end{array}\right.,$$

(30)

$$R_{\mathrm{sur}}^i=\left\{\begin{array}{cc}-{\beta }_{\mathrm{sur}}\left({\left({\vartheta }_u^{i,i+1}-{\displaystyle \frac{2\pi }{N_{\mathrm{uav}}}}\right)}^2+{\left({\vartheta }_u^{i-1,i}-{\displaystyle \frac{2\pi }{N_{\mathrm{uav}}}}\right)}^2\right),\hfill & if i is a Tracker\hfill \\ 0,\hfill & if i is the Target\hfill \end{array}\right..$$

(31)

As the identical definition in (1), ${\beta }_{\mathrm{trac}}$ and ${\beta }_{\mathrm{sur}}$ are gain factors. Note that the inverse proportional form of $R_{trac}^i$ is designed to guide the agents’ relative distance rate close to their maximum speed, which means encouraging to select the fastest tracking strategy. Furthermore, it is crucial for each UAV agent to maintain an encirclement angle of ${\vartheta }_u^{i-1,i}={\vartheta }_u^{i,i+1}={\displaystyle \frac{2\pi }{N_{\mathrm{uav}}}}$.

4.3. Algorithm

In this subsection, we will discuss the algorithm framework with structures of the actor and critic. Inspired by the Soft Actor–Critic (SAC) [46] that revealed excellent performance in continuous tasks, the soft Q-value and strategy update methods are adopted in our research. Furthermore, to leverage the time-serial data processing capability of GRU, we propose an AC structure that combines Multilayer Perceptrons (MLPs) and a GRU to extract valuable information from the state vector.

4.3.1. Learning Methods

Based on the recursive Bellman operator, the soft Q-value for each agent i is defined as

$$Q_s^i\left(s_t^i,a_t^i\right)=R^i\left(s_t^i,a_t^i\right)+\gamma V^i\left(s_t^i\right),$$

(32)

$$V^i\left(s_t^i\right)={\mathbb{E}}_{(s_{t+1}^i,a_{t+1}^i)}\left[Q_s^i\left(s_{t+1}^i,a_{t+1}^i\right)-{\alpha }^{i}ln\left({\pi }^i\left(a_t^i|s_t^i\right)\right)\right],$$

(33)

where ${\alpha }^i$ plays as the temperature coefficient to control the percentage of information entropy, keeping a balance between the exploration and exploitation behaviors of the algorithm. And $\gamma$ is the discount factor to influence the estimation length of the Q-value expectation. Given by the soft policy iteration, the optimal strategy is expressed as

$${\pi }^{i\ast }=\arg \max \sum _t^T\mathbb{E}\left[R^i\left(s_{\mathrm{t}},a_{\mathrm{t}}\right)+{\alpha }^i\mathcal{H}\left({\pi }^i\right)\right],$$

(34)

where $\mathcal{H}\left({\pi }^i\right)=-\int {\pi }^i\left(a_{\mathrm{t}}^i|s_{\mathrm{t}}^i\right)ln\left({\pi }^i\left(a_{\mathrm{t}}^i|s_{\mathrm{t}}^i\right)\right)d{a}_{\mathrm{t}}^i$ represents the information entropy. The above equation indicates that the optimal strategy should not only maximize rewards but also take into account entropy to promote exploration during the initial stages of training. There are two main streams in continuous MARL action sample methods, DDPG-like deterministic action output and Gaussian distribution action sample implemented in SAC. Compared with the hard output, probabilistic sampling brings more stochasticity to avoid losing effective actions.

As is typical of most MARL methods, TMSAC employs neural networks to parameterize the actor and the critic, allowing for gradient backpropagation during training. The parameters of AC are denoted as ${\theta }_a^i$ and ${\theta }_c^i$; also, ${\overline{\theta }}_c^i$ for target Critic prepared for updating. Meanwhile, a technique called Double-Critic is implemented, where each agent has two independent critic networks to calculate the soft Q-value. This approach helps to stabilize the training process. According to the Bellman residual definition, the loss function used in the stochastic gradient descent of critic has the following form:

$$\begin{array}{cc}\hfill J_{Q_s}^i\left({\theta }_c^i\right)=& {\displaystyle \frac{1}{2}}{\mathbb{E}}_{(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i)}\left[\left(Q_s^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)-R^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)+\right.\right.\hfill \\ \hfill & \left.{\left.\gamma {\mathbb{E}}_{(s_{\mathrm{t}+1}^i,a_{\mathrm{t}+1}^i)}\left[Q_{\overline{s}}^{i,\min }\left(s_{\mathrm{t}+1}^i,a_{\mathrm{t}+1}^i\right)\right]\right)}^2\right],\hfill \end{array}$$

(35)

where $Q_{\overline{s}}^{i,\min }=min\left\{Q_{\overline{s}}^{i,j}(s_{t+1}^i,a_{t+1}^i),j=1,2\right\}$ is the minimal target soft Q-value from Double-Critic. Thus, the update differential expression of the critic could be written as follows:

$$\begin{array}{cc}\hfill {\nabla }_{{\theta }_c^i}{J}_{Q_s}^i\left({\theta }_c^i\right)=& \mathbb{E}\left[{\nabla }_{{\theta }_c^i}{Q}_s^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)\left(Q_s^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)-\left(R^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)+\right.\right.\right.\hfill \\ \hfill & \left.\left.\left.\gamma \left(Q_{\overline{s}}^{i,\min }\left(s_{\mathrm{t}+1}^i,a_{\mathrm{t}+1}^i\right)-{\alpha }^{i}ln\left({\pi }^i\left(a_{t+1}^i|s_{t+1}^i\right)\right)\right)\right)\right)\right],\hfill \end{array}$$

(36)

where ∇ means the differential operator. In addition, a parameterized Gaussian distribution sampling function $G_{{\theta }_a^i}\left(s_t^i\right)$ is utilized to model the action outputs. The actor loss function and its derivative are indicated as the following form to realize soft policy iteration:

$$J_{\pi }^i\left({\theta }_a^i\right)={\mathbb{E}}_{(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i)}\left[{\alpha }^{i}ln\left({\pi }_{{\theta }_a^i}^i\left(G_{{\theta }_a^i}\left(s_t^i\right)|s_t^i\right)\right)-Q_s^{i,\min }\left(s_{\mathrm{t}}^i,G_{{\theta }_a^i}\left(s_t^i\right)\right)\right],$$

(37)

$$\begin{array}{cc}\hfill {\nabla }_{{\theta }_a^i}{J}_{\pi }^i\left({\theta }_a^i\right)=& \mathbb{E}\left[{\nabla }_{{\theta }_a^i}{\alpha }^{t}ln\left({\pi }_{{\theta }_a^i}^i\left(a_t^i|s_t^i\right)\right)+\right.\hfill \\ \hfill & \left.{\nabla }_{{\theta }_a^i}{G}_{{\theta }_a^i}\left(s_t^i\right)\left({\nabla }_{a^i}{\alpha }^{t}ln\left({\pi }_{{\theta }_a^i}^i\left(a_t^i|s_t^i\right)\right)-Q_{\mathrm{s}}^{i,min}\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)\right)\right],\hfill \end{array}$$

(38)

With the intention to reduce hyper-parameters of the algorithm, the ${\alpha }^i$ should adaptively update based on the following equation to increase exploration at the primary stage of training and restrain it when the entropy $\mathcal{H}$ approaches the target value $\widehat{\mathcal{H}}$:

$$J\left({\alpha }^i\right)=-{\alpha }^i{\mathbb{E}}_{(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i)}\left[ln\left({\pi }_{{\theta }_a^i}^i\left(a_{\mathrm{t}}^i|s_{\mathrm{t}}^i\right)\right)+\widehat{\mathcal{H}}\right].$$

(39)

4.3.2. Actor and Critic Structures

Unlike conventional approaches of using parameterized approximators in RL methods for TTE, our proposed method integrates recurrent units into the actor structure as illustrated in Figure 4a. Its inputs include the sequential relative kinematic vector $L\times {\kappa }^i$ and obstacle observation $O_{obs}^i$. Regarding the network architecture, the H1 layer of the actor structure, as shown in Figure 4a, is implemented as GRU and comprises 256 hidden units with $tanh$ activation function for new gates. This design is similar to that of Long Short-Term Memory (LSTM) but with fewer parameters, and has demonstrated comparable performance. MLP hidden layers H2 and H3 extract features from the upper layer outputs. Finally, the outputs of O1 and O2 represent the means and log standard deviations of Gaussian distribution, respectively, preparing for action sampling. On the other hand, the critic known as soft Q-value approximator should estimate the Q-value based on the current state. Therefore, a full MLP construction shown in Figure 4b is used to achieve this goal. Its input consists of the last kinematic step ${\kappa }^i$ and observation $O_{obs}^i$.

4.3.3. Details

The training process of TMSAC is shown in Algorithm 1. At the preparatory phase, all agents are born in the simulation environment, including the target with a stochastic initial speed and position. The agent i would be selected to train by the following switching rule in every episode:

$$f_{\mathrm{sw}}\left(ep\right)=\mathrm{div}\left(\mathrm{mod}\left(ep,M_{\mathrm{r}}\right),M_{\mathrm{r}}/N\right),$$

(40)

where $ep$ is the current episode. And $M_{\mathrm{r}}$ represents the switching period. The AC parameters of each agent are alternately optimized for the local Nash Equilibrium (NE) while the parameters of other agents are frozen, following the aforementioned rule. In the beginning of one episode, state set $\mathcal{S}$ is stochastically initialized with the condition that no collision occurred when agents locate their positions firstly.

After the initialization, agents execute actions $A_{\mathrm{t}}$ sampled from their Actors, respectively, and interact with the environment to obtain new state ${\mathcal{S}}_{\mathrm{t}+1}$ and rewards ${\mathcal{R}}_{\mathrm{t}}$. Meanwhile, experience buffer $\mathcal{D}$ is used to store the step information set $\left\{{\mathcal{S}}_{\mathrm{t}},A_{\mathrm{t}},R_{\mathrm{t}},{\mathcal{S}}_{\mathrm{t}+1}\right\}$ after the above interaction is accomplished. With the interval of $N_{\mathrm{g}}$ steps, agent i chose by Equation (41) update its AC parameters based on loss functions mentioned in (1). Note that the expectation $\mathbb{E}$ in loss functions is approximated by averaging a batch of data. The data flowchart in TMSAC is presented in Figure 5.

Algorithm 1 Time-Series Multi-Agent Soft Actor–Critic.

Input: Randomly initial ${\theta }_{\mathrm{c}}^{1,1}\cdots {\theta }_{\mathrm{c}}^{N,1}$, ${\theta }_{\mathrm{c}}^{1,2}\cdots {\theta }_{\mathrm{c}}^{N,2}$, ${\theta }_{\mathrm{a}}^1\cdots {\theta }_{\mathrm{a}}^N$ and ${\overline{\theta }}_{\mathrm{c}}^i\leftarrow {\theta }_{\mathrm{c}}^i$;    for $ep=1$ to M do:       Initialize the environment and agents state $\mathcal{S}$;       Initialize the experience buffer $\mathcal{D}$;       Select agent $i=f_{\mathrm{sw}}\left(ep\right)$ for training;       for $step=1$ to T do:           Sample actions $A_{\mathrm{t}}=\{a_{\mathrm{t}}^1,\cdots ,a_{\mathrm{t}}^N\}$ from each $G_{{\theta }_a^i}\left(s_t^i\right)$;           Execute $A_{\mathrm{t}}$ and update state set ${\mathcal{S}}_{\mathrm{t}+1}$ from environment;           Obtain rewards ${\mathcal{R}}_{\mathrm{t}}=\{R^1,\cdots ,R^N\}$;           Pull $\left\{{\mathcal{S}}_{\mathrm{t}},A_{\mathrm{t}},R_{\mathrm{t}},{\mathcal{S}}_{\mathrm{t}+1}\right\}$ into the experience buffer $\mathcal{D}$;           if $step\%N_{\mathrm{g}}=0$:               Gather a batch of set from the experience buffer $\mathcal{D}$;               Update ${\theta }_{\mathrm{c}}^{i,j}$, ${\theta }_{\mathrm{a}}^i$ and ${\alpha }^i$ of agent i through Equations (35), (37) and (39);               Softly update ${\overline{\theta }}_{\mathrm{c}}^i=\tau {\theta }_{\mathrm{c}}^i+(1-\tau ){\overline{\theta }}_{\mathrm{c}}^i$ of agent i;           end if       end for    end for    Output: ${\theta }_{\mathrm{a}}^1\cdots {\theta }_{\mathrm{a}}^N$;

4.3.4. Convergence Discussion

Although Definition 2 reveals that using MARL methods to solve an MDP is equivalent to finding the local NE of a proposed SPG scenario, it is necessary to prove the convergence of the TMSAC algorithm, which is a stochastic approximation tool. Given by policy iteration in Equation (34), the policy updating scheme can be rewritten as Equation (41):

$${\pi }_i^{\mathrm{new}}=arg\underset{{\pi }_i^{\prime }\in {\mathbf{\Pi }}_i}{min}{D}_{\mathrm{KL}}\left({\pi }_i^{\prime }\left(G_{{\theta }_a^i}\left(s_t^i\right)|s_t^i\right)\parallel exp\left({\displaystyle \frac{1}{\alpha }}\left(Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)-V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)\right)\right)\right).$$

(41)

where $D_{\mathrm{KL}}$ is the Kullback–Leibler divergence, and $Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)$ and $V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)$ are the soft state-action value and state value of agent i’s policy ${\pi }_i^{\mathrm{old}}$. Let ${\pi }_i^{\mathrm{old}}\in {\Pi }_i$, and ${\pi }_i^{\mathrm{new}}$ represents the optimizer of Equation (41). Only if $Q_s^{{\pi }_i^{\mathrm{new}}}\left(s_t^i,a_t^i\right)\ge Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)$ for all agents can the proposed algorithm potentially guarantee convergence towards the MPE.

Considering that $J_{{\pi }_i^{\mathrm{old}}}\left({\pi }_i^{\mathrm{old}}\right)\le J_{{\pi }_i^{\mathrm{old}}}\left({\pi }_i^{\mathrm{new}}\right)$, then the following inequality can be established:

$$\begin{array}{c}\hfill {\mathbb{E}}_{\left(s_t^i,a_t^i\sim {\pi }_i^{\mathrm{new}}\right)}\left[{\alpha }^{i}ln\left({\pi }_i^{\mathrm{new}}\left(G_{{\theta }_a^i}\left(s_t^i\right)|s_t^i\right)\right)-Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)+V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)\right]\\ \hfill \le {\mathbb{E}}_{\left(s_t^i,a_t^i\sim {\pi }_i^{\mathrm{old}}\right)}\left[{\alpha }^{i}ln\left({\pi }_i^{\mathrm{old}}\left(G_{{\theta }_a^i}\left(s_t^i\right)|s_t^i\right)\right)-Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)+V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)\right].\end{array}$$

(42)

As $V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)$ is only affected by $s_t^i$, when subtracting $V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)$ from both sides of the Inequality (42), it reduces to

$${\mathbb{E}}_{\left(s_t^i,a_t^i\sim {\pi }_i^{\mathrm{new}}\right)}\left[Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)-{\alpha }^{i}ln\left({\pi }_i^{\mathrm{new}}\left(G_{{\theta }_a^i}\left(s_t^i\right)|s_t^i\right)\right)\right]\ge {\mathbb{E}}_{\left(s_t^i\right)}\left[V^{{\pi }_i^{\mathrm{old}}}\left(s_t^i\right)\right].$$

(43)

Given by the soft Bellman operator mentioned in Section 4.3, the deductive inequality is shown as

$$\begin{array}{cc}\hfill Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_t^i,a_t^i\right)& =R^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)+\gamma V^{{\pi }_i^{\mathrm{old}}}\left(s_{t+1}^i\right)\hfill \\ \hfill & \le R^i\left(s_{\mathrm{t}}^i,a_{\mathrm{t}}^i\right)+\gamma {\mathbb{E}}_{s_{t+1}^i}\left[Q_s^{{\pi }_i^{\mathrm{old}}}\left(s_{t+1}^i,a_{t+1}^i\right)-{\alpha }^{i}ln{\pi }_i^{\mathrm{new}}\left(G_{{\theta }_a^i}\left(s_{t+1}^i\right)|s_{t+1}^i\right)\right]\hfill \\ \hfill & \le Q_s^{{\pi }_i^{\mathrm{new}}}\left(s_t^i,a_t^i\right),\hfill \end{array}$$

(44)

where the convergence towards the optimal policy of each agent is proven.

5. Simulation and Discussion

In this section, we conduct experiments using the proposed algorithm under various environmental and conditional settings. More importantly, we provide comparative analyses between TMSAC and other continuous MARL methods utilized in the SPG. The implementation is developed using Python and the PyTorch deep learning framework for constructing learning models, while the Matplotlib library is used for visualizing the results. The main computing hardware in the training platform is NVIDIA RTX 3090 GPU and AMD Ryzen 5950X CPU. Additionally, a visual simulation that takes into account the dynamic characteristics of the agents is provided, demonstrating the real-time deployability of the algorithm. Since the actor employs a relatively simple neural network architecture, it has low computational hardware requirements during the inference phase after training. This makes it suitable for deployment on embedded platforms.

5.1. Algorithm Performance

Considering the small-size quadrotor as trackers and target in all described simulation scenarios, the environment and agent capability parameters are indicated in

Table 1. Environment and agent parameters.

Names	Values
Plane size	1000 m × 1000 m
Tracker max speed $v_{\max }$	10 m/s
Tracker max angular speed ${\dot{\psi }}_{\max }$	15 deg/s
Target max speed $v_{\max }^t$	10 m/s
Target max angular speed ${\dot{\psi }}_{\max }^t$	10 deg/s
Trackers number $N_{\mathrm{UAV}}$	4 or 5
Target number	1
Agent safe interval $l_{\mathrm{int}}$	20 m
Obstacle detection distance $l_{\mathrm{obs}}$	8 m
Obstacle detection points $N_{\mathrm{obs}}$	12
Obstacles shape and size	random

. The algorithm hyper-parameters for agents are shown in

Table 2. Training hyper-parameters.

Names	Values
Total episodes M	$2\times {10}^4$
Total steps in each episode T	150
Switching period $M_{\mathrm{r}}$	20 × N
Learning interval $N_{\mathrm{g}}$	5
Batch size	1024
Target entropy $\widehat{\mathcal{H}}$	−2
Discount factor $\gamma$	0.988
Buffer size	$2^{17}$
Sequential length	10
Optimizer	Adam

. To verify the effectiveness of the algorithm under different environmental conditions, we first conduct simulations in an obstacle-free environment. The MUAVS consists of four to five UAVs. Based on the performance parameters of small quadrotors, the maximum flight speed of each agent is set to 10 m/s, and the maximum angular velocity is set to 15 deg/s. We deploy the MUAVS in a 1000 m × 1000 m two-dimensional simulation area without any obstacle disturbing tracking process, to verify the TTE performance preliminarily. It is important to emphasize that the initial positions and orientations of both the UAV swarm and the target are randomly generated, and the agent’s safe distance $l_{\mathrm{int}}$ is set to 20 m to mitigate the risk of collisions.

As the verification results show in Figure 6, four typical tracking and encirclement evaluation scenarios are presented within four or five trackers after the training. On the trackers side, their main objective is to proactively converge towards the target and encircle it in a coordinated manner. The punishment incurred from the safe interval reward has incentivized the UAVs in Figure 6c to learn to encircle the target instead of getting too close to it, resulting in higher rewards. This demonstrates the practicality and intelligence of the algorithm. On the target side, it exhibits an evasion tendency while trackers are getting closer in all presented scenarios that conform to a rational response in a real pursuit–evasion situation. Under four different initialization conditions, all TTE tasks are completed within 100 steps (seconds). UAVs selected flight speeds greater than 10 m/s in the early stages of tracking to quickly approach the target, and matched speed with the target after approaching.

More importantly, the TTE performance of TMSAC in the unknown environment is revealed in Figure 6 which employs four types of obstacles distribution in simulation while four UAVs are involved in the tracker group. In Figure 6e,f, there are several isolated geometrical objects deployed in the area, respectively, meanwhile trackers and evader are born in the rest of the open area. On the other hand, narrow channels with different humps are set in the environment shown in Figure 6g,h. The results demonstrate that the agents effectively avoid collisions during both the pursuit and evasion phases.

5.2. Comparison with Other Methods

Intending to demonstrate the effectiveness and efficiency of the proposed algorithm, we also implement two additional baseline methods, MADDPG and Multi-Agent Soft Actor–Critic (MASAC), with nearly identical hyper-parameters in the same environment. In particular, the MASAC has the same algorithm framework as MADDPG but replaces the loss function using the Max Entropy Learning theory. Both of the two methods, MADDPG and MASAC, adopt fully connected MLPs as their actors and critics, which have similar structures to TMSAC but do not include any GRU layers. Different numbers of UAVs are deployed in the simulation as trackers, and training rewards are shown in Figure 7, where UAV rewards are the average values of all swarm members. According to the rewards figures, there are several learning characteristics of algorithms that can be summarized as follows:

• Higher rewards value. In all cases of comparison, TMSAC obtains 5–20% more rewards than other methods with the same environment and hyper-parameters settings, both for trackers and target, simultaneously. That means strategies approximated by the TMSAC algorithm are closer to MPE as mentioned in Section 3.
• Adaptivity for more agents scenarios. For fully distributed algorithms, increasing the number of agents makes it more challenging to search for NE strategies. From Figure 7a–e, the advantage of TMSAC is progressively obvious on the trackers side. The recurrent network structure gives a measure of prediction ability for agent to avoid potential collision, and select effective actions for the predictive tracking. The mounting agents group brings more kinematic information into state vectors. The results demonstrate that as the number of trackers in the environment increases, TMSAC exhibits enhanced robustness in identifying near-Nash Equilibrium policies.

After evaluating the training efficiency, we analyze the success rate and the number of steps required for the successful tracking process under different conditions. To evaluate the algorithms, we define two performance metrics:

• Full Success Rate (FSR): This metric measures the percentage of scenarios in which all trackers in the MUAVS successfully capture their targets within a limited number of time steps.
• Half Success Rate (HSR): This metric quantifies the percentage of scenarios, in which at least half of the trackers achieve their goals.

The criteria for tracking success are specified as $d_{\mathrm{t}}^i\in l_{\mathrm{int}}\pm 5\mathrm{m}$. And the time steps usage represents the percentage of tracking success steps over total steps T. As shown in Figure 8a, with the random obstacles, our algorithm achieves higher scores averaged from 100 test episodes both in the FSR and TSR metrics. Meanwhile, as can be seen in Figure 8b, the slight advantage of TMSAC in time steps usage is verified.

Furthermore, we evaluate the performance of different algorithms as the number of UAVs is increased. As shown in Figure 9, we gradually increase the number of deployed UAVs in an environment with random obstacles. The TMSAC algorithm demonstrates superior scalability compared to the other two algorithms, maintaining a stable task success rate as the number of game participants is increased. Specifically, when the number of UAVs is less than 10, the success rate of TTE tasks executed by navigators trained with TMSAC and MASAC algorithms is not significantly affected. When the number of UAVs exceeds 10, the environmental complexity begins to significantly impact the task success rates of all three algorithms. When the number of UAVs reaches 20, the TMSAC task success rate drops by 8% compared to the conditions of the 5 UAVs, MASAC drops by 10%, and MADDPG drops by nearly 20%. It should be noted that all three algorithms use the CEDT architecture, and increasing the number of agents does not affect the algorithms’ convergence. However, under the same training conditions, different algorithms show varying levels of performance.

5.3. Visual Simulation

Based on the above works, we build a comprehensive dynamic simulation platform derived from gym-pybullet-drones [47], which is a multifunctional end-to-end training environment for UAV intelligent control methods. It uses the PyBullet physical engine as a dynamics solver. The proficiently trained actors function as navigators in the platform, complemented by the use of Proportional–Integral–Derivative (PID) attitude and velocity controllers to establish a closed-loop negative feedback control system. For the real-time dynamic simulation,

Table 3. Physical parameters.

Names	Values
Gravity acceleration g	9.81
UAV mass m	0.25 kg
UAV arm length l and b	0.2 m
X-axis inertias $J_{\mathrm{x}}$	$1.4\times {10}^{-4}$ $\mathrm{kg}·{\mathrm{m}}^2$
Y-axis inertias $J_{\mathrm{y}}$	$1.4\times {10}^{-4}$ $\mathrm{kg}·{\mathrm{m}}^2$
Z-axis inertias $J_{\mathrm{z}}$	$2.2\times {10}^{-4}$ $\mathrm{kg}·{\mathrm{m}}^2$
Max propeller force	7.23 N

describes the physical parameters of agents both for the MUAVS and target.

The system block diagram is shown in Figure 10, where the SPG-based navigator receives embedded kinematics and environment observation, and outputs the expectations of plane velocity $v_{\mathrm{ref}}^i$ and yaw velocity ${\dot{\psi }}_{\mathrm{ref}}^i$. The classic PID velocity and attitude controller calculate propeller forces that simplify the expression without motor rotation speed. Then, the UAV dynamic model introduced in Section 2 updates three axes’ velocities and angular velocities for the SPG relative kinematics. Given by this framework, the visual simulation displayed in Figure 11 is executed at an iteration frequency of (a) 120 Hz and a control frequency of (b) 60 Hz to evaluate the performance of the actors under real-time conditions.

In one iteration of the virtual experiment, the dynamic responses of all agents’ positions and attitudes in all three axes are recorded and presented in Figure 12. The initial conditions are exactly the same as those of the above simulations. It can be seen that both the X and Y positions converge to smaller biases, with velocities and angular velocities controlled within reasonable ranges. Note that the yaw angle is defined in $[-\pi ,\pi ]$; therefore, there exist some situations of numerical mutation in the yaw angle records.

6. Conclusions

Facing the unknown environment and target evasion problems, we proposed the TMSAC, a multi-agent reinforcement learning algorithm to approximate agent strategies that converge to the MPE in the SPG abstracted from TTE scenarios. By integrating a reward function design compliant with the SPG paradigm and a state-action space design meeting practical deployment requirements, our proposed algorithm demonstrates exceptional performance in MUAVS executing TTE tasks, while also showing potential for scalability and large-scale applications. Furthermore, a series of simulations was conducted to validate the characteristics of the proposed algorithm. These simulations involved performance tests and virtual dynamic experiments. Simulation results demonstrate that even under conditions of unknown target dynamic escape and environmental disturbances, MUAVS can still effectively track the target and achieve circumnavigation flight. Comparison of the training rewards demonstrated that the proposed TMSAC algorithm outperformed the two baseline algorithms in terms of both training rewards and tracking success rate. Additionally, the dynamic experiments provided evidence of the real-time deployability of the TMSAC method.

In terms of future work, the key issue is to enhance approaches for solving the Nvs.M TTE problem by incorporating task allocation and target-choosing mechanisms. This will broaden the scope of approaches and enable them to operate in more complex environments. By incorporating hierarchical learning techniques, the model and algorithm proposed in this paper show potential for further expansion and research.