A Novel Method for a Pursuit–Evasion Game Based on Fuzzy Q-Learning and Model-Predictive Control

Abstract

This paper explores a pursuit–evasion game (PEG) based on quadrotors by combining fuzzy Q-learning (FQL) and model-predictive control (MPC) algorithms. Initially, the FQL algorithm is employed to perceive, make decisions, and predict the trajectory of the evader. Based on the position and velocity information of both players in the game, the pursuer quadrotor determines its action strategy using the FQL algorithm. Subsequently, a state feedback controller is designed using the MPC algorithm, with reference inputs derived from the FQL algorithm. Within each MPC cycle, the FQL algorithm dynamically provides reference inputs to the MPC, thereby enhancing its robust control and dynamic optimization for the quadrotor. Finally, simulation results verify the effectiveness of the proposed algorithm.

1. Introduction

A pursuit–evasion game (PEG) involves two opposing groups: the pursuers and the evaders. It mainly focuses on how pursuers cooperate to capture evaders and what strategies evaders should adopt to avoid capture or prolong the capture time [1]. In recent years, PEGs have been widely applied in military and civil fields, such as missile defense and interception [2], air combat [3], search and rescue [4], and transportation management [5]. Various forms of PEGs have emerged, including single-to-single [6], single-to-multiple [7], multiple-to-single [8], multiple-to-multiple [9], and other scenarios. At the same time, a variety of algorithms have been presented and discussed, such as the geometric method [10,11], the differential game method [12], the classical control method [13], and the reinforcement learning (RL) method [14].

Compared with other classical control algorithms, the model-predictive control (MPC) algorithm has garnered significant research attention in PEG problems due to its flexible optimization capabilities, constraint-handling ability, adaptability, and robustness. For a PEG of fixed-wing unmanned aerial vehicles (UAVs), the authors proposed a nonlinear MPC algorithm that incorporated the safety and functional constraints of the UAV into the cost function, achieving autonomous obstacle avoidance and pursuit maneuvers [15]. Furthermore, the MPC algorithm was applied to the PEG of unicycle models in an obstacle environment, and the robustness and security of the algorithm were demonstrated [16]. Considering obstacles, the author proposed a novel PEG algorithm that can alternately convert strategies between game theory and MPC algorithms, effectively reducing computational complexity [17]. Considering the constraint of limited information, a PEG algorithm based on MPC was proposed, which enables players to use limited information to predict the opponent’s movement strategy in an obstacle environment [18]. To tackle the two-target, two-attacker PEG problem with state and control constraints, a game strategy based on nonlinear MPC and Apollonius circle (AC) was developed [19]. Additionally, based on the unmanned surface vessel platform, a robust adversarial PEG algorithm based on MPC was proposed by combining the particle swarm optimization algorithm with the cost function [20].

Traditional PEG algorithms, such as MPC, typically rely on precise mathematical models for designing control laws. However, practical applications often involve random factors such as sensor errors and environmental disturbances, which make it challenging to acquire accurate models and bring significant challenges to the control of agents. In this context, RL algorithms have become an important method to solve PEG problems. A distributed deep Q network (DQN) algorithm was proposed to enhance the collaborative pursuit or evasion strategy between agents through a decentralized learning framework [21]. The authors proposed a distributed cooperative PEG algorithm based on RL and designed communication networks to save communication and computational resources [22]. Kartal et al. [23] applied an integral reinforcement learning algorithm to achieve optimal action selection for online and real-time PEG scenarios. In [24], they introduced a novel two-stage pursuit strategy employing the deep deterministic policy gradient (DDPG) algorithm to resolve PEG challenges using real-time feedback information, particularly under conditions of incomplete information. Zhang et al. [25] have integrated a prediction network into an RL framework, enhancing the algorithm’s performance and agents’ predictive capabilities in the PEG scenario. Selvakumar et al. [26] utilized min–max Q-learning and matrix game theory to transform a non-zero-sum PEG into multiple two-player static zero-sum PEGs and then determined the optimal actions for agents in these complex scenarios.

Traditional RL algorithms, deeply rooted in rigorous mathematical foundations, provide a framework for the precise modeling and analysis of complex problems. However, they often encounter challenges related to interpretability and adaptability. Fuzzy reinforcement learning has emerged as a promising solution to effectively handle imprecise and ambiguous information, especially in contexts where uncertain environmental perception and decision-making processes are prevalent. The inherent flexibility of fuzzy logic enables it to adeptly handle intricate relationships, proving particularly advantageous when dealing with complexities like non-linearity and non-convexity that arise in diverse problem domains. Notably, the rules and inference processes produced by fuzzy reinforcement learning tend to be intuitive and easily comprehensible. Combining the strengths of both fuzzy reinforcement learning and traditional approaches yields promising results. Consequently, a large number of RL algorithms incorporating fuzzy logic principles have emerged. Based on the quadrotor platform, a novel control algorithm was proposed that combines fuzzy logic control with RL techniques and can adapt to varied initial positions and noise conditions [27]. The author utilized the fuzzy Q-learning (FQL) algorithm to acquire control policies for agents and also used the formation control mechanism to capture a superior evader, effectively addressing the multiple-to-single PEG problem [28]. Yue et al. [29] combined the fuzzy inference system and RL to address the PEG problem in the continuous action domain, and the proposed algorithm can be applied to high-dimensional situations. A PEG technique combining fuzzy systems with Q-learning was proposed for UAVs, which solved the challenges of continuous state and action spaces in Q-learning [30].

In summary, the MPC algorithm relies on precise mathematical models and target reference inputs. The traditional RL algorithm suffers from poor interpretability and adaptability. Fuzzy RL algorithm can effectively deal with imprecise and fuzzy information. Therefore, in this paper, we integrate the FQL algorithm with the MPC algorithm to address the PEG problem based on quadrotors and provide a robust framework for its solution. We utilize the perception and decision-making capabilities of FQL to formulate and evaluate strategies for the quadrotor. Meanwhile, we leverage the robust constraint handling and dynamic optimization capabilities of the MPC algorithm to achieve robust control for the quadrotor. First, we employ the FQL algorithm to perform perception and decision-making tasks, while evaluating the target’s position and velocity information. Fuzzy logic is employed to process ambiguous information and generate appropriate action strategies, enabling the quadrotor to make suitable decisions in diverse PEG scenarios for both tracking and evasion. Then, we utilize the action strategy generated by FQL as a reference state for MPC to design a state feedback controller. During each control cycle, we minimize an objective function that considers control input constraints and obstacle avoidance constraints, ensuring that the quadrotor does not violate physical constraints and safety requirements in the PEG scenario. The main contributions of this paper are as follows:

(1)

For the first time, the Apollonius circle is extended from 2D space to 3D space, and an analytical equation is provided. A novel learning algorithm based on the FQL framework is proposed, and a reward function based on the idea of artificial potential field is designed for the learning algorithm, which improves the convergence speed and learning performance of the algorithm.

(2)

Distinct from existing methods, with the strategy derived from the FQL algorithm as the reference signal, we model the trajectory tracking problem as an MPC optimization problem for decoupled quadrotor. We designed a state feedback controller for the quadrotor that takes into account control input constraints and obstacle avoidance constraints and analyzed the feasibility and stability of the solution.

(3)

We have verified the learning performance of FQL in 3D scenarios. The trajectory tracking capability of the designed state feedback algorithm has been validated through 3D trajectory tracking results. Finally, the combination of FQL and MPC algorithms is employed to control the quadrotor to achieve a PEG on the Gazebo platform.

The rest of this paper is organized as follows. Section 2 introduces the perception and decision based on the FQL algorithm. In Section 3, we present the quadrotor control based on the MPC algorithm. Section 4 provides detailed simulation results to verify the effectiveness of the algorithm. Finally, Section 5 concludes the paper.

2. Perception and Decision Based on Fuzzy Q-Learning

2.1. The Model of Pursuit–Evasion Game

The pursuer generates its trajectory prediction by assessing factors such as the evader’s position and velocity. We employ the FQL algorithm to learn the evader’s movement strategy and make decisions. In this section, we only require knowledge of the evader’s control strategy without involving specific dynamic controls introduced in the following sections. Therefore, to facilitate the learning of control strategies, we omit the complex control model of the quadrotor and focus on a simpler model. As shown in Figure 1, the agent’s kinematic model is described as (1)

$$\begin{array}{cc}\hfill x_{t+1}& =x_t+vsin\theta cos\alpha \hfill \\ \hfill y_{t+1}& =y_t+vsin\theta sin\alpha \hfill \\ \hfill z_{t+1}& =z_t+vcos\theta ,\hfill \end{array}$$

(1)

where $(x_t,y_t,z_t)$ and v represent the position and the velocity of agent, respectively. $\alpha$ is the angle between the projection of velocity v onto the x-y plane and the x-axis, and $\theta$ is the angle between velocity v and the z-axis. $U={[\alpha ,\theta ]}^{\top }$ is the steering angle of the agent and also the output of the controller. To ensure that the agent’s motion adheres to real-world constraints, the increments $\Delta \alpha$ and $\Delta \theta$ of the agent’s steering angle are limited in the interval of $\left[-{\displaystyle \frac{\pi }{4}},{\displaystyle \frac{\pi }{4}}\right]$. The ability of agents to engage in capture or evasion depends on the unique input information of the algorithm, namely, the positions and velocities of agents. Consequently, we introduce the Apollonius circle to depict the agent’s dominant area and determine their movement strategies based on their geometric relationships.

It is assumed that the pursuer’s velocity $v_P$ is faster than the evader’s velocity $v_E$, that is $\left|v_P\right|>\left|v_E\right|$. As shown in Figure 2, we extended the Apollonius circle from 2D space to its generalized form in 3D space, providing the analytical equation of the generalized Apollonius circle. A, B, C and D are points located on the generalized Apollonius circle. The lines $PA$ and $PB$ are tangent to sphere $O_{AC}$, that is, $PA\perp A{O}_{AC}$ and $PB\perp B{O}_{AC}$, separating the capture and the escape regions of agents. Given the positions of the pursuer and the evader denoted as $(x_t^P,y_t^P,z_t^P)$ and $(x_t^E,y_t^E,z_t^E)$, respectively, and the speed ratio a of the agents, we have (2)

$$a={\displaystyle \frac{\left|v_E\right|}{\left|v_P\right|}}={\displaystyle \frac{\left|EA\right|}{\left|PA\right|}}={\displaystyle \frac{\left|EB\right|}{\left|PB\right|}}={\displaystyle \frac{\left|EC\right|}{\left|PC\right|}}<1,$$

(2)

where $\left|EA\right|$ denotes the Euclidean distance between the position of the evader E and the point A. The center and radius of the generalized Apollonius circle can be calculated by (3)

$$\begin{array}{cc}\hfill O_{AC}& =\left({\displaystyle \frac{x_t^E-a^2{x}_t^P}{1-a^2}},{\displaystyle \frac{y_t^E-a^2{y}_t^P}{1-a^2}},{\displaystyle \frac{z_t^E-a^2{z}_t^P}{1-a^2}}\right)\hfill \\ \hfill R_{AC}& =\left|{\displaystyle \frac{a}{1-a^2}}\right|\sqrt{{\left(x_t^P-x_t^E\right)}^2+{\left(y_t^P-y_t^E\right)}^2+{\left(z_t^P-z_t^E\right)}^2}.\hfill \end{array}$$

(3)

Within the space enclosed by the yellow sphere in Figure 2, the evader consistently arrives earlier than the pursuer, establishing it as the evader’s dominant region. Conversely, the external area represents the pursuer’s dominant region.

To ensure that the pursuer always moves towards capturing the evader and that the evader constantly avoids the pursuer’s pursuit, we define the agents’ motion strategies. For a random point C on the generalized Apollonius circle, we have (4)

$${\displaystyle \frac{sin{\alpha }_1}{sin{\alpha }_2}}={\displaystyle \frac{\left|v_P\right|}{\left|v_E\right|}},$$

(4)

where ${\alpha }_1$ is the angle difference between the vector $\overrightarrow{EP}$ and the vector $\overrightarrow{EC}$, and ${\alpha }_2$ is the angle difference between the vector $\overrightarrow{PE}$ and the vector $\overrightarrow{PC}$. We obtain (5)

$${\alpha }_2={sin}^{-1}\left({\displaystyle \frac{\left|v_E\right|}{\left|v_P\right|}}sin{\alpha }_1\right).$$

(5)

According to (5), when $sin{\alpha }_1=1$, the boundary value of the optimal solution for ${\alpha }_2$ is

$${\alpha }_2^{*}={sin}^{-1}\left({\displaystyle \frac{\left|v_E\right|}{\left|v_P\right|}}\right).$$

(6)

Therefore, the optimal motion space for both agents lies within the conical region. In a PEG scenario with $\left|v_P\right|>\left|v_E\right|$, the pursuer moves within the envelope defined by the cone, satisfying the condition ${\alpha }_2\le {\alpha }_2^{*}={sin}^{-1}\left({\displaystyle \frac{\left|v_E\right|}{\left|v_P\right|}}\right)$ to efficiently capture the evader. As for the evader, it maneuvers within the remaining cone envelope, satisfying ${\alpha }_1\ge {\displaystyle \frac{\pi }{2}}$. The terminal conditions are defined as follows: the PEG ends at time $t_f$ if the distance between the pursuer P and the evader E satisfies $d_{PE}\left(t_f\right)\le d_s$, where $d_s$ represents the capture distance. Additionally, if the game duration exceeds the designated threshold, the game is deemed over.

2.2. Fuzzy Q-Learning Algorithm

In this paper, the PEG environment is not a typical discrete state–action space but rather a continuous environment. Thus, we need an approximator to map the state into an action. For classical RL, the goal is to obtain the optimal expected return to evaluate the performance of the control strategy. The expected value function is

$$Q(s,a)=\mathbb{E}\left\{\sum _{k=0}^{\infty }{\gamma }^t{r}_{t+k+1}|s_t=s,a_t=a\right\}.$$

(7)

In the Q-learning algorithm, it is necessary to consider the challenges arising from dealing with continuous state or action spaces. Therefore, the FQL algorithm, which generates global continuous actions for agents based on a predefined discrete action set, can effectively address the aforementioned issues [14]. We summarize the update process of the FQL algorithm as shown in Figure 3. We assume that the agent has n inputs $\overline{x}=[x_1,\cdots ,x_n]$, m actions $A=\{a_1,\cdots ,a_m\}$, and L fuzzy rules. For each fuzzy rule $l(l=1,\cdots ,L)$, it can be described by an IF-THEN statement

$$R_l:IF{x}_{1}is{F}_1^l,\cdots ,AND{x}_{n}is{F}_n^{l}THEN{f}^l=a^l,$$

(8)

where $x_i$ is the ith input of the $\overline{x}$ and $f^l$ is the rule’s output related to the corresponding fuzzy variable. $F_i^l$ is the fuzzy set related to the input $x_i$, and $a^l$ is the action selected from the discrete action set A for rule l. During the training process, it is necessary to design a reasonable action selection mechanism to overcome the exploration–exploitation dilemma.

In the Softmax exploration strategy, the action with higher Q-value has a higher probability of being selected. For rule l, the probability of an action $a^l$ being selected is (9)

$$Pr\left(a^l\right)={\displaystyle \frac{exp\left(\frac{Q(l,a^l)}{\mathcal{T}}\right)}{{\displaystyle \sum _{k=1}^{\left|A\right|}}exp\left(\frac{Q(l,a^k)}{\mathcal{T}}\right)}},$$

(9)

where $Q(l,a^l)$ is the Q-value of action $a^l$ given rule l and $\left|A\right|$ is the size of the action space. $\mathcal{T}$ is the softmax temperature. To reduce the algorithmic complexity and adapt to different task requirements, another action-selection mechanism is the $\epsilon -\mathrm{greedy}$ method, which is described as (10)

$$a^l=\left\{\begin{array}{cc}\forall a\in A& \epsilon \\ \underset{a\in A}{argmax}Q(l,a)& 1-\epsilon \end{array}\right.,$$

(10)

where $\epsilon$ is the probability of randomly selecting action $a\in A$ given rule l. The probability of selecting action a with the maximum Q-value is $1-\epsilon$. After the action is selected for each rule, the global continuous action at time t is

$$a_t\left({\overline{x}}_t\right)={\displaystyle \frac{{\displaystyle \sum _{l=1}^L}\left(\left({\displaystyle \prod _{i=1}^n}{\mu }^{F_i^l}\left(x_i\right)\right)·a^l\right)}{{\displaystyle \sum _{l=1}^L}\left({\displaystyle \prod _{i=1}^n}{\mu }^{F_i^l}\left(x_i\right)\right)}}={\displaystyle \sum _{l=1}^L}{\Phi }_t^l{a}_t^l,$$

(11)

where ${\Phi }_t^l$ is the firing strength for rule l at time t. The firing strength serves to quantify the degree to which the rule is activated under given input conditions, determining which rules will be prioritized and utilized to generate the final output. In a fuzzy logic system, the process begins by fuzzifying the system’s input variables through membership functions, yielding corresponding membership values for each input variable. Subsequently, these membership values are matched against the conditions of the rules within the fuzzy rule base. The firing strength of a rule is then calculated based on whether the conditions of the rule are met. The firing strength is defined as follows

$${\Phi }_t^l={\displaystyle \frac{{\displaystyle \prod _{i=1}^n}{\mu }^{F_i^l}\left(x_i\right)}{{\displaystyle \sum _{l=1}^L}\left({\displaystyle \prod _{i=1}^n}{\mu }^{F_i^l}\left(x_i\right)\right)}},$$

(12)

where ${\mu }^{F_i^l}$ is the membership degree of the fuzzy set $F_i^l$, and it can be calculated by the Gaussian membership function or the triangular membership function. The global Q-function is

$$Q_t\left({\overline{x}}_t\right)=\sum _{l=1}^L{\Phi }_t^l{Q}_t(l,a_t^l),$$

(13)

where $Q_t(l,a_t^l)$ is the Q-value after taking action $a_t^l$ for rule l at time step t. The global Q-function with the maximum Q-value is

$$Q_t^{*}\left({\overline{x}}_t\right)=\sum _{l=1}^L{\Phi }_t^l\underset{a\in A}{max}{Q}_t(l,a).$$

(14)

The temporal difference (TD) error is defined as

$${\tilde{\epsilon }}_{t+1}=r_{t+1}+\gamma Q_t^{*}\left({\overline{x}}_{t+1}\right)-Q_t\left({\overline{x}}_t\right),$$

(15)

where $\gamma$ is the discount factor and $r_{t+1}$ is the reward at time step t. The update law for the Q-function is

$$\begin{array}{cc}\hfill Q_{t+1}(l,a_t^l)& =Q_t(l,a_t^l)+{\alpha }_q{\tilde{\epsilon }}_{t+1}{\Phi }_t^l\hfill \\ \hfill & =Q_t(l,a_t^l)+{\alpha }_q\left[r_{t+1}+\gamma Q_t^{*}(l,a_{t+1}^l)-Q_t(l,a_t^l)\right]{\Phi }_t^l,\hfill \end{array}$$

(16)

where ${\alpha }_q$ is the learning rate.

In this paper, we define four inputs for the FQL used by each agent, and the inputs for the pursuer are (17)

$${\overline{x}}_P=[d_{PE},{\delta }_P,d_{PO},{\delta }_{PO}],$$

(17)

where $d_{PE}$ represents the distance between the pursuer and the evader. The term ${\delta }_P$ denotes the angle difference between the heading of the pursuer ${\overrightarrow{v}}_P\left(t\right)$ and a straight line from the pursuer to the evader $\overrightarrow{PE}\left(t\right)$. Additionally, $d_{PO}$ stands for the distance between the pursuer and the obstacle, and ${\delta }_{PO}$ indicates the angle difference between the heading of the pursuer ${\overrightarrow{v}}_P\left(t\right)$ and a straight line from the pursuer to the obstacle $\overrightarrow{PO}\left(t\right)$. For the evader, the inputs of the FQL system are defined as (18)

$${\overline{x}}_E=[d_{PE},{\delta }_E,d_{EO},{\delta }_{EO}],$$

(18)

where ${\delta }_E$ denotes the angle difference between the heading of the evader ${\overrightarrow{v}}_E\left(t\right)$ and a straight line from the evader to the pursuer $\overrightarrow{EP}\left(t\right)$. $d_{EO}$ is the distance between the evader and the obstacle. The term ${\delta }_{EO}$ indicates the angle difference between the heading of the evader ${\overrightarrow{v}}_E\left(t\right)$ and a straight line from the evader to the obstacle $\overrightarrow{EO}\left(t\right)$.

As we all know, the reward function plays a crucial role in RL, as it significantly influences both the convergence speed and overall performance of the RL algorithm. An appropriate reward function can greatly assist the agent in acquiring an accurate strategy. For instance, by allowing the reward to increase as the agent approaches the optimal solution, the algorithm can converge quickly. To achieve better learning performance, we need to use prior knowledge to design a reasonable reward function for the agent.

In this paper, we design a reward function based on the idea of the artificial potential field. The obstacle exerts repulsive force on the pursuer, while the evader exerts attraction on the pursuer. We model the repulsion exerted by obstacles on the pursuer using an exponential function. The reward function for repulsion can be designed as

$$r_{PO}=1-exp(-{\alpha }_r\ast \Delta d_{PO}),$$

(19)

where $\Delta d_{PO}=d_{PO}(t+1)-d_{PO}\left(t\right)$, and ${\alpha }_r$ is the repulsion coefficient controlling the strength of the repulsive force. Smaller distance results in higher repulsion and consequently smaller reward. The evader has a similar form of reward function when approaching obstacles. The reward function for attraction can be formulated as

$$r_{PE}=exp\left(-{\beta }_a\ast \Delta d_{PE}\right)-1,$$

(20)

where $\Delta d_{PE}=d_{PE}(t+1)-d_{PE}\left(t\right)$, and ${\beta }_a$ is the attraction coefficient governing the strength of the attractive force. Smaller distance results in stronger attraction and larger rewards. For the evader, the reward function (20) is inverted. Figure 4 depicts the curves of reward functions (19) and (20) under different coefficients. It can be observed that when the independent variable, i.e., the distance increment, is zero, the value of the reward function is also zero, which aligns with the requirement for realistic reward values. We incorporate an additional term to the reward function that considers whether the pursuer successfully captures the evader

$$r_s={\gamma }_s\ast g_s,$$

(21)

where ${\gamma }_s$ is the coefficient for the success-based reward, and $g_s$ is an indicator function that $g_s=1$ when the pursuer captures the evader and $g_s=0$ otherwise. This component incentivizes the pursuer to complete its objective successfully. We formulate a comprehensive reward function using a weighted balance between these forces.

$$r_{total}=w_r\ast r_{PO}+w_a\ast r_{PE}+r_s,$$

(22)

where $w_r$ and $w_a$ represent the weights for balancing the influence of repulsion and attraction. We can adjust these weights based on experience, domain knowledge, and the specific requirements of the practical problem. By adjusting these weights, we can control the overall behavior of the reward function, enabling the pursuer to pursue the evader flexibly in a complex environment.

Throughout the learning process, each agent selects an action based on the current state. After taking the action and receiving the corresponding reward, the Q-function updates by (16). The learning process then continues until the final state is reached, that is, the evader is captured by the pursuer or the maximum learning time is reached. Using the FQL algorithm, we obtain the motion information of the evader, including position and velocity information. We define the trajectory of the evader as a parameterized regular curve in the output space of the position system

$${\mathcal{P}}_E=\left\{{\widehat{p}}_t\in {\mathbb{R}}^3|{\widehat{p}}_t={[x_r\left(t\right),y_r\left(t\right),z_r\left(t\right)]}^{\top }\right\},$$

(23)

where ${\widehat{p}}_t$ is a time-dependent trajectory parameter. Our objective is to design a robust controller for the pursuer quadrotor with control constraints, enabling the stable and precise tracking of the trajectory ${\mathcal{P}}_E$ with evolving parameters ${\widehat{p}}_t$ of the evader quadrotor. We utilize FQL and MPC algorithms to collaboratively control the quadrotors in the PEG scenario. The control framework of the quadrotor based on FQL and MPC is shown in Figure 5. We employed the FQL algorithm to train and acquire the agent’s control strategy. Then, the obtained motion trajectory of the evader from FQL served as a reference input for the MPC algorithm. In each control cycle, the MPC algorithm combines the current state of the quadrotor, including position and velocity, as well as the reference trajectory provided by FQL, to generate control signals for the quadrotor.

3. Quadrotor Control Based on MPC

3.1. The Model of Quadrotor and Control Objective

As shown in Figure 6, in the inertial coordinate system $\mathcal{I}=\left\{I_x,I_y,I_z\right\}$, the position of the quadrotor’s center of mass is $\ell ={[x,y,z]}^{\top }$. The linear velocity and the angular velocity are expressed as $\dot{\eta }={[\dot{x},\dot{y},\dot{z},\dot{\psi }]}^{\top }$. In the body’s coordinate system $\mathcal{B}=\left\{B_x,B_y,B_z\right\}$, the velocity and angular velocity are expressed as $v_B={[u,v,w,r]}^{\top }$. In the initial state, we assume that the body’s coordinate system coincides with the inertial coordinate system. As the quadrotor’s movement state changes, the body’s coordinate system and the inertial coordinate system will no longer coincide. To align the body’s coordinate system with the inertial coordinate system, it is necessary to rotate the inertial coordinate system around the coordinate axis $I_x,I_y$, and $I_z$, respectively, by the roll angle $\varphi$, pitch angle $\theta$, and yaw angle $\psi$. In this paper, we do not consider the complex dynamic model based on $Euler--Lagrange$ equations, but only focus on the translational motion of the quadrotor. We have

$$\dot{\eta }=\mathbf{R}\left(\psi \right)v_B,$$

(24)

where $\mathbf{R}\left(\psi \right)=\left[\begin{array}{cccc}cos\left(\psi \right)& -sin\left(\psi \right)& 0& 0\\ sin\left(\psi \right)& cos\left(\psi \right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$. We design a four-degree-of-freedom quadrotor model that is decoupled in terms of horizontal and altitude control, as follows:

$$\dot{\mathbf{x}}=\left[\begin{array}{c}\dot{\eta }\\ {\dot{v}}_B\end{array}\right]=\left[\begin{array}{c}\mathbf{R}\left(\psi \right)v_B\\ \mathbf{F}{v}_B+\mathbf{G}\mathbf{u}\end{array}\right]=\mathbf{g}(\mathbf{x},\mathbf{u}),$$

(25)

where $\mathbf{F}=diag\left({\displaystyle \frac{-1}{{\tau }_x}},{\displaystyle \frac{-1}{{\tau }_y}},{\displaystyle \frac{-1}{{\tau }_z}},{\displaystyle \frac{-1}{{\tau }_{\psi }}}\right)$ and $\mathbf{G}=diag\left({\displaystyle \frac{k_x}{{\tau }_x}},{\displaystyle \frac{k_y}{{\tau }_y}},{\displaystyle \frac{k_z}{{\tau }_z}},{\displaystyle \frac{k_{\psi }}{{\tau }_{\psi }}}\right)$ are system identification parameters used to reduce the impact of system complexity and the lagged response of the quadrotor. ${\tau }_x,{\tau }_y,{\tau }_z,{\tau }_{\psi }$ are time constants, and $k_x,k_y,k_z,k_{\psi }$ are gains. $\mathbf{u}$ is the control signal. To facilitate the design of the control law, we present the following assumptions.

Assumption 1. 

The quadrotor’s body is considered a rigid structure, and other external force disturbances are neglected. All control parameters of the quadrotor can be accurately measured.

Assumption 2. 

The reference trajectory ${\widehat{p}}_t$ satisfies the following conditions: $\left|x_r\left(t\right)\right|\le \overline{x}$, $\left|y_r\left(t\right)\right|\le \overline{y}$, $\left|z_r\left(t\right)\right|\le \overline{z}$, $\left|{\dot{x}}_r\left(t\right)\right|\le {\overline{x}}_1$, $\left|{\dot{y}}_r\left(t\right)\right|\le {\overline{y}}_1$, $\left|{\dot{z}}_r\left(t\right)\right|\le {\overline{z}}_1$, $\left|{\ddot{x}}_r\left(t\right)\right|\le {\overline{x}}_2$, $\left|{\ddot{y}}_r\left(t\right)\right|\le {\overline{y}}_2$, $\left|{\ddot{z}}_r\left(t\right)\right|\le {\overline{z}}_2$. The reference system state increment is ${\mathbf{x}}_r\left(t\right)=\left[{\eta }_r\left(t\right);v_{Br}\left(t\right)\right]$; we have (26)

$$\begin{array}{cc}\hfill {\eta }_r\left(t\right)& ={\left[x_r\left(t\right),y_r\left(t\right),z_r\left(t\right),{\psi }_r\left(t\right)\right]}^{\top }\hfill \\ \hfill v_{Br}\left(t\right)& ={\left[u_r\left(t\right),v_r\left(t\right),w_r\left(t\right),r_r\left(t\right)\right]}^{\top }\hfill \end{array}$$

(26)

Lemma 1. 

The reference state ${\eta }_r\left(t\right)$, as well as its first-order derivative ${\dot{\eta }}_r\left(t\right)$ and the second-order derivative ${\ddot{\eta }}_r\left(t\right)$, are all bounded.

Proof. 

According to (24), we have $\dot{\psi }=r_r$, and it satisfies the following relation:

$$\begin{array}{cc}\hfill \left|{\dot{\psi }}_r\left(t\right)\right|& \le {\displaystyle \frac{{\overline{x}}_1{\overline{y}}_2+{\overline{x}}_2{\overline{y}}_1}{{\overline{x}}_1^2+{\overline{y}}_1^2}}\hfill \\ \hfill \left|{\ddot{\psi }}_r\left(t\right)\right|& \le {\displaystyle \frac{{\overline{x}}_1{\overline{y}}_3+{\overline{x}}_3{\overline{y}}_1}{{\overline{x}}_1^2+{\overline{y}}_1^2}}+{\displaystyle \frac{2({\overline{x}}_1{\overline{x}}_2+{\overline{y}}_1{\overline{y}}_2)({\overline{x}}_1{\overline{y}}_2+{\overline{x}}_2{\overline{y}}_1)}{{({\overline{x}}_1^2+{\overline{y}}_1^2)}^2}}.\hfill \end{array}$$

(27)

If Assumption 2 holds, we have

$$\begin{array}{cc}\hfill \overline{\eta }& =max\left\{\overline{x},\overline{y},\overline{z},{\psi }_r\right\}\hfill \\ \hfill {\overline{\eta }}_1& =max\left\{{\overline{x}}_1,{\overline{y}}_1,{\overline{z}}_1,{\displaystyle \frac{{\overline{x}}_1{\overline{y}}_2+{\overline{x}}_2{\overline{y}}_1}{{\overline{x}}_1^2+{\overline{y}}_1^2}}\right\}\hfill \\ \hfill {\overline{\eta }}_2& =max\left\{{\overline{x}}_2,{\overline{y}}_2,{\overline{z}}_2,{\displaystyle \frac{{\overline{x}}_1{\overline{y}}_3+{\overline{x}}_3{\overline{y}}_1}{{\overline{x}}_1^2+{\overline{y}}_1^2}}+{\displaystyle \frac{2({\overline{x}}_1{\overline{x}}_2+{\overline{y}}_1{\overline{y}}_2)({\overline{x}}_1{\overline{y}}_2+{\overline{x}}_2{\overline{y}}_1)}{{({\overline{x}}_1^2+{\overline{y}}_1^2)}^2}}\right\}.\hfill \end{array}$$

(28)

By combining (26) and (28), we can deduce that ${\eta }_r,{\dot{\eta }}_r,{\ddot{\eta }}_r$ are bounded and satisfy

$${∥{\eta }_r∥}_{\infty }\le \overline{\eta },{∥{\dot{\eta }}_r∥}_{\infty }\le {\overline{\eta }}_1,{∥{\ddot{\eta }}_r∥}_{\infty }\le {\overline{\eta }}_2.$$

(29)

The proof is completed.    □

If the reference trajectory of the evader and the model of quadrotor are given by (23) and (25), respectively, the control problem for the pursuer quadrotor can be described as Problem 1.

Problem 1. 

Given a reference trajectory depicted by (23), the system described by (25), the control objective is to compute appropriate control inputs $∥\mathbf{u}\left(t\right)∥$ that satisfy

$$\underset{t\to t+T}{lim}∥\ell \left(t\right)-{\widehat{p}}_t∥=0,$$

(30)

where T is the control cycle and $\ell \left(t\right)$ is the position of the quadrotor. The control signal satisfies constraint ${∥\mathbf{u}\left(t\right)∥}_{\infty }\le u_{max}$, and $u_{max}$ is the maximum amplitude of the control signal.

3.2. Design of the State Feedback Controller Based on MPC

We model the tracking of the quadrotor to the reference trajectory in Problem 1 as an MPC optimization problem:

$$\begin{array}{}\mathrm{(31a)}& \underset{\widehat{\mathbf{u}}(·)}{min}\overline{J}& =& {\int }_{t_k}^{t_k+T}{∥\tilde{\mathbf{x}}(\tau ;t_k)∥}_Q^2+J\left(\widehat{\mathbf{x}}(\tau ;t_k)\right)+{∥\widehat{\mathbf{u}}(\tau ;t_k)∥}_R^{2}d\tau \hfill \mathrm{(31b)}& & s.t.& \dot{\widehat{\mathbf{x}}}(\tau ;t_k)=\mathbf{g}(\widehat{\mathbf{x}}(\tau ;t_k),\widehat{\mathbf{u}}(\tau ;t_k))\mathrm{(31c)}& & & \widehat{\mathbf{x}}(t_k;t_k)=\mathbf{x}\left(t_k\right)\mathrm{(31d)}& & & {∥\widehat{\mathbf{u}}(\tau ;t_k)∥}_{\infty }\le u_{max},\tau \in [t_k,t_k+T]\mathrm{(31e)}& & & {\displaystyle \frac{\partial V}{\partial \mathbf{x}}}\mathbf{g}(\widehat{\mathbf{x}}(t_k;t_k),\widehat{\mathbf{u}}(t_k;\widehat{\mathbf{x}}(t_k;t_k)))\le {\displaystyle \frac{\partial V}{\partial \mathbf{x}}}\mathbf{g}(\widehat{\mathbf{x}}(t_k;t_k),\mathbf{h}(t_k;\widehat{\mathbf{x}}(t_k;t_k))),\end{array}$$

where $\tilde{x}\left(\tau \right)=\widehat{x}\left(\tau \right)-x_r\left(\tau \right)$ is the trajectory tracking error. Q and R are weight matrices. $J\left(x\right(t\left)\right)$ is a cost function for obstacle avoidance:

$$J\left(x\left(t\right)\right)={\displaystyle \frac{\lambda }{1+e^{\Delta d_o\left(t\right)}}}-{\displaystyle \frac{\lambda }{2}},$$

(32)

where $\lambda$ is an adjustable weight parameter representing the penalty for collisions, $\Delta d_o\left(t\right)=d_o^2-{∥\ell \left(t\right)-p_o∥}_2^2$, where $d_o$ is the safety distance for obstacle avoidance and $p_o$ is the location of the obstacle. The function $\mathbf{h}\left(\mathbf{x}\right)$ is a tracking control law based on Lyapunov, and V is the corresponding Lyapunov function. (31e) represents the stability constraint term, which dynamically adjusts the prediction horizon by comparing the program calculation time and the sampling period. The introduction of the stability constraint term can ensure the stability of the control system while achieving a dynamic adjustment between control performance and computational efficiency. We define the sampling period as $\delta =T/N$, where N represents the prediction horizon. The update process of the trajectory tracking and collision-avoidance controller is shown in Figure 7. First, we construct the objective function $\overline{J}$ and obtain the state measurement $\mathbf{x}\left(t_k\right)$ of the quadrotor. Second, we set ${\widehat{\mathbf{x}}}_0=\mathbf{x}\left(t_k\right)$ as the initial prediction value, solve the MPC optimization problem (31), and obtain the optimal control sequence $U^{*}={\left[{\widehat{\mathbf{u}}}_0^{*},{\widehat{\mathbf{u}}}_1^{*},\cdots ,{\widehat{\mathbf{u}}}_{N-1}^{*}\right]}^{\top }$. During the sampling period, the control signal $\mathbf{u}\left(t\right)={\widehat{\mathbf{u}}}_0^{*}$ is used to control the quadrotor. Third, the increment of the prediction horizon is calculated as $\Delta N=\left|t_{cpu}-T\right|$, where $t_{cpu}$ represents the program processing time. If $t_{cpu}>\delta$, let $N=N-\Delta N$; if $t_{cpu}<\delta$, let $N=N+\Delta N$; otherwise N remains unchanged. Finally, update $k\leftarrow k+1$ and $t_{k+1}\leftarrow t_k+\delta$.

Below, we present the design and stability analysis of the state feedback control law $\mathbf{h}\left(\mathbf{x}\right)$. We design a Lyapunov function as follows:

$$V_{\ell }={\displaystyle \frac{1}{2}}{\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell },$$

(33)

where ${\mathbf{e}}_{\ell }={\eta }_r-\eta$ is the position error. The first derivative of $V_{\ell }$ is

$${\dot{V}}_{\ell }={\mathbf{e}}_{\ell }^{\top }{\dot{\mathbf{e}}}_{\ell }={\mathbf{e}}_{\ell }^{\top }\left[{\dot{\eta }}_r-\mathbf{R}\left(\psi \right)v_B\right]+\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell }-\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell },$$

(34)

where $\omega >0$ is the parameter for the control law $\mathbf{h}\left(\mathbf{x}\right)$. We define the velocity error as

$${\mathbf{e}}_v={\dot{\eta }}_r-\mathbf{R}\left(\psi \right)v_B+\omega {\mathbf{e}}_{\ell }={\dot{e}}_{\ell }+\omega e_{\ell }.$$

(35)

The first derivative of ${\mathbf{e}}_v$ is

$$\begin{array}{cc}\hfill {\dot{\mathbf{e}}}_v& ={\ddot{\eta }}_r-\dot{\mathbf{R}}\left(\psi \right)v_B-\mathbf{R}\left(\psi \right){\dot{v}}_B+\omega {\dot{\mathbf{e}}}_{\ell }\hfill \\ \hfill & ={\ddot{\eta }}_r-\dot{\mathbf{R}}\left(\psi \right)v_B-\mathbf{R}\left(\psi \right){\dot{v}}_B+\omega ({\mathbf{e}}_v-\omega {\mathbf{e}}_{\ell })\hfill \\ \hfill & ={\ddot{\eta }}_r-\dot{\mathbf{R}}\left(\psi \right)v_B-\mathbf{R}\left(\psi \right)(\mathbf{F}{v}_B+\mathbf{Gu})+\omega ({\mathbf{e}}_v-\omega {\mathbf{e}}_{\ell }).\hfill \end{array}$$

(36)

By substituting (35) into (34), we obtain

$${\dot{V}}_{\ell }={\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_v-\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell }.$$

(37)

The Lyapunov function V is chosen as (38)

$$V=V_{\ell }+{\displaystyle \frac{1}{2}}{\mathbf{e}}_v^{\top }{\mathbf{e}}_v.$$

(38)

The first derivative of V is

$$\begin{array}{cc}\hfill \dot{V}& ={\dot{V}}_{\ell }+{\mathbf{e}}_v^{\top }{\dot{\mathbf{e}}}_v\hfill \\ \hfill & =-\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell }+{\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_v+{\mathbf{e}}_v^{\top }{\dot{\mathbf{e}}}_v\hfill \\ \hfill & =-\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell }-\rho {\mathbf{e}}_v^{\top }{\mathbf{e}}_v+\rho {\mathbf{e}}_v^{\top }{\mathbf{e}}_v+{\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_v+{\mathbf{e}}_v^{\top }{\dot{\mathbf{e}}}_v,\hfill \end{array}$$

(39)

where $\rho$ is another parameter of the control law $\mathbf{h}\left(\mathbf{x}\right)$. To ensure that under the control law $\mathbf{h}\left(\mathbf{x}\right)$, the Lyapunov function V satisfies (40),

$$\dot{V}(\mathbf{x},\mathbf{h}\left(\mathbf{x}\right))=-\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell }-\rho {\mathbf{e}}_v^{\top }{\mathbf{e}}_v\le 0.$$

(40)

We deduce the control law $\mathbf{h}\left(\mathbf{x}\right)$ as

$$\mathbf{h}\left(\mathbf{x}\right)=-{\mathbf{G}}^{-1}\mathbf{F}{v}_B-{\mathbf{G}}^{-1}{\mathbf{R}}^{\top }\left(\psi \right)\mu ,$$

(41)

where $\mu$ is defined as follows

$$\mu ={\ddot{\eta }}_r-\dot{\mathbf{R}}\left(\psi \right)v_B+(\omega +\rho ){\mathbf{e}}_v+(1-{\omega }^2){\mathbf{e}}_{\ell }.$$

(42)

Under the action of the control signal $\mathbf{u}=\mathbf{h}\left(\mathbf{x}\right)$, the derivative of the Lyapunov function is

$$\dot{V}(\mathbf{x},\mathbf{u})=-\omega {\mathbf{e}}_{\ell }^{\top }{\mathbf{e}}_{\ell }-\rho {\mathbf{e}}_v^{\top }{\mathbf{e}}_v+{\mathbf{e}}_v^{\top }[\mu -\mathbf{R}\left(\psi \right)(\mathbf{F}{v}_B+\mathbf{Gu})].$$

(43)

Therefore, the stability constraint can be expressed as

$${\widehat{\mathbf{e}}}_v^{\top }(t_k;t_k)\left[\widehat{\mu }(t_k;t_k)-\mathbf{R}\left(\widehat{\psi }(t_k;t_k)\right)\left(\mathbf{F}{\widehat{v}}_B(t_k;t_k)+\mathbf{G}\widehat{\mathbf{u}}(t_k;\widehat{\mathbf{x}}(t_k;t_k))\right)\right]\le 0.$$

(44)

3.3. Feasibility Analysis of the State Feedback Control Law

Under the condition of satisfying the input constraint ${∥\mathbf{h}\left(\mathbf{x}\right)∥}_{\infty }\le u_{max}$, we prove the existence of a feasible solution to the MPC optimization problem using Lemma 2 and Theorem 1.

Lemma 2. 

According to Lemma 1, under the control law $\mathbf{u}=\mathbf{h}\left(\mathbf{x}\right)$, the velocity $v_B$ of the quadrotor satisfies

$${∥v_B∥}_{\infty }\le \sqrt{2}\left({\overline{\eta }}_1+{∥{\mathbf{e}}_v\left(t_0\right)∥}_2+\omega {∥{\mathbf{e}}_{\ell }\left(t_0\right)∥}_2\right),$$

(45)

where ${\mathbf{e}}_v\left(t_0\right)$ and ${\mathbf{e}}_{\ell }\left(t_0\right)$ represent the position error and velocity error at time $t_0$, respectively.

Proof. 

According to (24), we have

$${∥v_B∥}_{\infty }={∥{\mathbf{R}}^{\top }\left(\psi \right)\dot{\eta }∥}_{\infty }\le {∥{\mathbf{R}}^{\top }\left(\psi \right)∥}_{\infty }{∥\dot{\eta }∥}_{\infty }\le \sqrt{2}{∥\dot{\eta }∥}_{\infty }.$$

(46)

According to Lemma 1, as well as ${\mathbf{e}}_{\ell }={\eta }_r-\eta$ and ${\dot{\mathbf{e}}}_{\ell }={\mathbf{e}}_v-\omega {\mathbf{e}}_{\ell }$, we have

$${∥\dot{\eta }∥}_{\infty }=∥{\dot{\eta }}_r-{\dot{\mathbf{e}}}_{\ell }∥\le {\overline{\eta }}_1+{∥{\dot{\mathbf{e}}}_{\ell }∥}_{\infty }\le {\overline{\eta }}_1+{∥{\mathbf{e}}_v∥}_{\infty }+\omega {∥{\mathbf{e}}_{\ell }∥}_{\infty }.$$

(47)

Due to $\dot{V}<0$, there exist ${∥{\mathbf{e}}_{\ell }\left(t_k\right)∥}_{\infty }\le {∥{\mathbf{e}}_{\ell }\left(t_k\right)∥}_2\le {∥{\mathbf{e}}_{\ell }\left(t_0\right)∥}_2$ and ${∥{\mathbf{e}}_v\left(t_k\right)∥}_{\infty }\le {∥{\mathbf{e}}_v\left(t_k\right)∥}_2\le {∥{\mathbf{e}}_v\left(t_0\right)∥}_2$. Therefore, according to ${∥\dot{\eta }∥}_{\infty }\le {\overline{\eta }}_1+{∥{\mathbf{e}}_v\left(t_0\right)∥}_2+\omega {∥{\mathbf{e}}_{\ell }\left(t_0\right)∥}_2$, we can conclude that (45) holds. The proof is completed.    □

Theorem 1. 

Let $\overline{g}={∥{\mathbf{G}}^{-1}∥}_{\infty }$ and $\overline{f}={∥\mathbf{F}∥}_{\infty }$, the control parameters ω and ρ for control law $\mathbf{h}\left(\mathbf{x}\right)$ are positive constants. If (48) is satisfied,

$$\sqrt{2}·\overline{g}\left(\overline{f}·l+{\overline{\eta }}_2+2\sqrt{2}·l^2+m\right)\le u_{max},$$

(48)

where $l={\overline{\eta }}_1+{∥{\mathbf{e}}_v\left(t_0\right)∥}_2+\omega {∥{\mathbf{e}}_{\ell }\left(t_0\right)∥}_2$ and $m=(\omega +\rho ){∥{\mathbf{e}}_v\left(t_0\right)∥}_2+(1-{\omega }^2){∥{\mathbf{e}}_{\ell }\left(t_0\right)∥}_2$. Then, the optimization problem (31) has a feasible solution.

Proof. 

According to (41), we have

$${∥\mathbf{h}\left(\mathbf{x}\right)∥}_{\infty }\le {∥{\mathbf{G}}^{-1}∥}_{\infty }{∥\mathbf{F}∥}_{\infty }{∥v_B∥}_{\infty }+{∥{\mathbf{G}}^{-1}∥}_{\infty }{∥{\mathbf{R}}^{\top }\left(\psi \right)∥}_{\infty }{∥\mu ∥}_{\infty }.$$

(49)

According to the definition of $\mathbf{R}\left(\psi \right)$, we have

$${∥{\mathbf{R}}^{\top }\left(\psi \right)∥}_{\infty }=max\left\{\left|sin\left(\psi \right)\right|+\left|cos\left(\psi \right)\right|,1\right\}\le \sqrt{2}.$$

(50)

According to (42), we have

$${∥\mu ∥}_{\infty }\le {\overline{\eta }}_2+{∥-\dot{\mathbf{R}}\left(\psi \right)∥}_{\infty }{∥v_B∥}_{\infty }+(\omega +\rho ){∥{\mathbf{e}}_v∥}_{\infty }+(1-{\omega }^2){∥{\mathbf{e}}_{\ell }∥}_{\infty },$$

(51)

where $-\dot{\mathbf{R}}\left(\psi \right)=\left[\begin{array}{cccc}sin\left(\psi \right)r& cos\left(\psi \right)r& 0& 0\\ -cos\left(\psi \right)r& sin\left(\psi \right)r& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$, and ${∥-\dot{\mathbf{R}}\left(\psi \right)∥}_{\infty }\le \sqrt{2}{∥v_B∥}_{\infty }$. According to Lemma 2, we can obtain

$${∥\mu ∥}_{\infty }\le {\overline{\eta }}_2+2\sqrt{2}·l^2+m.$$

(52)

By substituting (45), (50), and (52) into (49), we obtain

$${∥\mathbf{h}\left(\mathbf{x}\right)∥}_{\infty }\le \sqrt{2}·\overline{g}\left(\overline{f}·l+{\overline{\eta }}_2+2\sqrt{2}·l^2+m\right).$$

(53)

Therefore, if (48) is satisfied, we have ${∥\mathbf{h}\left(\mathbf{x}\right)∥}_{\infty }\le u_{max}$. Consequently, a feasible solution exists for the MPC optimization problem (31). The proof is completed.    □

Through the above analysis, we know that there exists a feasible solution to the optimization problem (31). Next, we analyze the stability of the closed-loop system.

3.4. Stability Analysis of the State Feedback Control Law

We analyze the stability of the closed-loop system using the following Theorem 2.

Theorem 2. 

The state feedback control law $\mathbf{h}\left(\mathbf{x}\right)$ enables the system tracking trajectory to converge asymptotically to the reference trajectory.

Proof. 

For the Lyapunov function V in (38), there exists a positive definite function ${\mathcal{F}}_{i=1,2,3}(·)$ that satisfies the following inequality [31]:

$$\begin{array}{cc}\hfill & {\mathcal{F}}_1\left(∥\mathbf{x}∥\right)\le V\left(\mathbf{x}\right)\le {\mathcal{F}}_2\left(∥\mathbf{x}∥\right)\hfill \\ \hfill & {\displaystyle \frac{\partial V\left(\mathbf{x}\right)}{\partial \mathbf{x}}}\mathbf{g}(\mathbf{x},\mathbf{h}\left(\mathbf{x}\right))\le -{\mathcal{F}}_3\left(∥\mathbf{x}∥\right).\hfill \end{array}$$

(54)

The Lyapunov function of the closed-loop system satisfies

$${\displaystyle \frac{\partial V\left(\mathbf{x}\right)}{\partial \mathbf{x}}}\mathbf{g}(\mathbf{x},\mathbf{u})\le {\displaystyle \frac{\partial V\left(\mathbf{x}\right)}{\partial \mathbf{x}}}\mathbf{g}(\mathbf{x},\mathbf{h}\left(\mathbf{x}\right))\le -{\mathcal{F}}_3\left(∥\mathbf{x}∥\right).$$

(55)

According to the Lyapunov theory, under the control algorithm proposed in this paper, the closed-loop system is asymptotically stable at the equilibrium point $\tilde{\mathbf{x}}=\mathbf{0}$. The proof is completed.    □

Finally, we provide the overall algorithmic flowchart of the paper, as shown in Algorithm 1.

Algorithm 1 PEG Algorithm Based on FQL and MPC

1: Initialization: inputs ${\overline{x}}_0$, rule $l(l=1,2,\cdots ,L)$, action $A=\{a_1,\cdots ,a_m\}$, learning rate ${\alpha }_q$, discount factor $\gamma$, reward function r, and Q-function $Q(·)=0$   2: for each time step t do   3:    Obtain inputs ${\overline{x}}_t=[x_1,\cdots ,x_n]$ at time t   4:    Choose an action $a^l$ for rule l based on (9) or (10)   5:    Calculate the global continuous action $U_t\left({\overline{x}}_t\right)$ by (11)   6:    Calculate the global Q-function $Q_t\left({\overline{x}}_t\right)$ by (13)   7:    Take the global action $U_t\left({\overline{x}}_t\right)$   8:    Obtain the reward $r_{t+1}$ and the inputs ${\overline{x}}_{t+1}$ at time $t+1$   9:    Calculate the maximum Q-value $Q_t^{*}\left({\overline{x}}_t\right)$ by (14) 10:    Calculate the TD error ${\tilde{\epsilon }}_{t+1}$ by (15) 11:    Update the Q-function $Q_{t+1}(l,a_t^l)$ by (16) for rule l 12: end for 13: Get the motion information of the evader ${\mathcal{P}}_E$ with the trajectory parameter ${\widehat{p}}_t$ by (23) 14: for control instant $t_k$ do 15:    Construct the objective function $\overline{J}$ and define $\delta =T/N$ 16:    Obtain the state measurement $\mathbf{x}\left(t_k\right)$ 17:    Let ${\widehat{\mathbf{x}}}_0=\mathbf{x}\left(t_k\right)$, and solve the MPC optimization problem (31) 18:    Obtain the optimal control sequence $U^{*}={\left[{\widehat{\mathbf{u}}}_0^{*},{\widehat{\mathbf{u}}}_1^{*},\cdots ,{\widehat{\mathbf{u}}}_{N-1}^{*}\right]}^{\top }$ 19:    Execute the control signal $\mathbf{u}\left(t\right)={\widehat{\mathbf{u}}}_0^{*}$ and record program processing time $t_{cpu}$ 20:    Calculate the increment of the prediction horizon $\Delta N=\left|t_{cpu}-T\right|$ 21:    if $t_{cpu}>\delta$ then 22:      $N=N-\Delta N$ 23:    else if $t_{cpu}<\delta$ then 24:      $N=N+\Delta N$ 25:    else 26:      $N=N$ 27:    end if 28:    Update $k\leftarrow k+1$ and $t_{k+1}\leftarrow t_k+\delta$. 29: end for

4. Simulation Results and Analysis

First, we evaluated the learning performance of the FQL algorithm in different PEG scenarios. Then, we validated the effectiveness of the MPC algorithm in controlling the quadrotor. Finally, the performance of the FQL and MPC algorithms to jointly control the quadrotor PEG was verified in the Gazebo platform.

4.1. Simulation of Quadrotor PEG Based on FQL

We define a 3D environment for the PEG, where agents move in a space with a size of $35\mathrm{m}\times 35\mathrm{m}\times 15\mathrm{m}$ at maximum speeds of $v_P=1.1\mathrm{m}/\mathrm{s}$ and $v_E=1\mathrm{m}/\mathrm{s}$, respectively. According to (22), the parameters of the reward functions are set as ${\alpha }_r=10,{\beta }_a=5,{\gamma }_s=20,w_r=5,w_a=10$. The number of episodes is set to 200, and the number of plays in each game is 500. The sample time is $0.1\mathrm{s}$. The discount factor is $\gamma =0.95$, and the learning rate is ${\alpha }_q=0.001$. The game terminates when the pursuer captures the evader, that is, $d_{PE}\le d_s=0.5\mathrm{m}$, or when the time exceeds $100\mathrm{s}$. To reduce the computational load, we select the triangular membership function. Specifically, the FQL has five triangular membership functions for each input. The distance inputs are in the interval of $[0,35]$, and the angle inputs are in the interval of $[-\pi ,\pi ]$. Consequently, we have $5\times 5\times 5\times 5=625$ rules. In the 3D environment, we select the controller’s output from a set of tuples $\left\{(\alpha ,\theta )|\Delta \alpha ,\Delta \theta \in [-{\displaystyle \frac{\pi }{4}},{\displaystyle \frac{\pi }{4}}]\right\}$.

We set up four different scenarios, in which the pursuer starts from the initial position $[7.5,7.5,0]$, and the evader needs to start from the position $[7.5,27.5,0]$, avoid obstacles, and then reach the target point $[30,5,10]$. The objective of the pursuer is to capture the evader before it reaches the target point. The obstacles are randomly placed and represented by cuboids of dimensions $1.5\mathrm{m}\times 1.5\mathrm{m}\times 15\mathrm{m}$. The trajectories of pursuers and evaders are illustrated in Figure 8, where Figure 8a–d represent the 3D trajectories of agents, while Figure 8e–h represent the projections of the 3D trajectories onto the x-y plane. We can observe that in the 3D environment, the pursuer successfully captured the evader as well, without any collisions occurring. As shown in Figure 9, we calculated the distances between the pursuer and the evader in four scenarios. After the PEG ends, the distances between the pursuer and the evader are $0.42\mathrm{m}$, $0.44\mathrm{m}$, $0.32\mathrm{m}$, and $0.35\mathrm{m}$, respectively, all of which are less than the set distance $d_s=0.5\mathrm{m}$.

To verify the superiority of the reward function designed in this paper, we compared the learning performance of the algorithm under different reward functions. In the same scenarios, we compared the performance of the reward function proposed in [32] with that designed in this paper. We repeatedly executed the FQL algorithm with two different reward functions, 20 times each, and then calculated the average path length and capture time of the pursuer, as shown in Figure 10. It can be observed that in all four scenarios, the reward function designed in this paper outperforms the one in [32] in terms of both path length and capture time. The primary reason is that the reward function designed in this paper is a nonlinear function that better aligns with the movement status of the agent. As shown in Figure 4, when the distance between the pursuer and the evader decreases, the reward value increases exponentially rather than linearly, which helps motivate the pursuer to capture the evader more quickly. The simulation results validate the effectiveness of the FQL algorithm and the reward function designed in this paper.

4.2. Simulation of Quadrotor Control Based on MPC

We define the parameters of the quadrotor and the MPC algorithm as follows. For the quadrotor parameters, the gains are set to $k_x=k_y=k_z=k_{\psi }=1$, and the time constants are set to ${\tau }_x=0.8$, ${\tau }_y=0.7$, ${\tau }_z=0.5$, ${\tau }_{\psi }=0.5$. For the MPC algorithm parameters, the sampling period is $\delta =0.1\mathrm{s}$, the weight matrices are $Q=diag({10}^3,{10}^3,{10}^3,{10}^3,{10}^2,{10}^2,{10}^2,{10}^2)$ and $R=diag(1,1,1,1)$, and the parameters for the control law are $\omega =0.5$ and $\rho =1$. For the cost function for obstacle avoidance, the adjustable weight parameter representing the penalty for collisions is $\lambda =500$, and the safety distance for obstacle avoidance is $d_o=0.5\mathrm{m}$. To verify the effectiveness of the altitude controller based on the MPC, we compared the proposed algorithm with the standard LQR algorithm [33,34,35]. We assumed that the quadrotor takes off from $(0,0,0)$ at the origin of the coordinate system and flies along the reference trajectory described in (56)

$${\mathcal{P}}_E=\left\{{\widehat{p}}_t\in {\mathbb{R}}^3\left|\widehat{p}=\left[{\displaystyle \frac{8\ast cos\left(\frac{\pi }{20}t\right)}{1+{sin}^2\left(\frac{\pi }{20}t\right)}},{\displaystyle \frac{8\ast sin\left(\frac{\pi }{20}t\right)cos\left(\frac{\pi }{20}t\right)}{1+{sin}^2\left(\frac{\pi }{20}t\right)}},2\ast (1-cos\left({\displaystyle \frac{\pi }{20}}t\right))\right]\right.\right\}.$$

(56)

The entire simulation process lasted for $40\mathrm{s}$. The result of 3D trajectory tracking of the quadrotor is shown in Figure 11, where (a) represents the 3D trajectory and (b) is the projection of the 3D trajectory onto the x-y plane. The circle represents the starting position, and the star represents the destination position. It can be seen that both algorithms can guide the quadrotor to track the reference trajectory. As shown in Figure 12, we analyzed the convergence of position and angle tracking. It is clear that the algorithm proposed in this paper exhibits less overshoot and a better response speed.

To further verify the stability of the proposed MPC algorithm, we introduced external disturbances into the simulation and compared it with the standard LQR algorithm. Based on the quadrotor model (25), the height control and horizontal control are decoupled. Therefore, we consider the anti-interference performance of the control algorithm in the horizontal plane. We assume that the quadrotor starts from the initial position $(2,1)$ and flies along the reference trajectory described in (57):

$$\begin{array}{cc}\hfill x_t& =2.5+1.5·cos({\displaystyle \frac{\pi }{20}}t-{\displaystyle \frac{\pi }{2}})\hfill \\ \hfill y_t& =2.5+1.5·sin({\displaystyle \frac{\pi }{20}}t-{\displaystyle \frac{\pi }{2}}).\hfill \end{array}$$

(57)

The entire simulation process lasted for $40\mathrm{s}$. External interferences with a period of $2\mathrm{s}$ and an amplitude of $0.5\mathrm{m}$ were introduced into the yaw channel at $t=12\mathrm{s}$ and $t=25\mathrm{s}$. As shown in Figure 13, the quadrotor moved in a counterclockwise direction as indicated by the arrow in the figure. We analyzed the performance of the proposed MPC and standard LQR algorithms in controlling the quadrotor’s trajectory tracking under external disturbances. The trajectory tracking curve demonstrates that although both algorithms exhibit excellent trajectory tracking performance and a notable capability to suppress external interference, our proposed MPC algorithm exhibits less overshoot and smaller tracking error. We analyzed the quadrotor’s trajectory tracking results along the coordinate axes, as shown in Figure 14. Specifically, Figure 14a,b present the trajectory tracking results, Figure 14c,d display the trajectory tracking errors, Figure 14e,f show the velocity tracking results, and Figure 14g,h represent the velocity tracking errors. Through Figure 14, it can be seen that our proposed algorithm exhibits superior robustness and resistance to external disturbances.

4.3. Simulation of Quadrotor PEG Based on FQL and MPC

We validated the performance of using FQL and MPC to control quadrotors in PEG on the Gazebo platform. The hardware configuration consisted of an Intel i7-9700 processor and a GeForce RTX 3090 graphics card. The learning algorithm operated on a central computer, which subsequently transmitted the strategies to quadrotors. The MPC algorithm ran in MATLAB, and the control of quadrotors was implemented in the Ubuntu 18.04 operating system based on the robot operating system (ROS). The scenario size was set to $35\mathrm{m}\times 35\mathrm{m}\times 15\mathrm{m}$, and the obstacles were represented by cuboids of dimensions $1.5\mathrm{m}\times 1.5\mathrm{m}\times 15\mathrm{m}$. The parameters of the learning algorithm, MPC algorithm, and the quadrotor were the same as the configurations described in Section 4.1 and Section 4.2. Figure 15 presents screenshots of key moments in four simulation scenarios. We selected one quadrotor as the pursuer, marked by the red circle, and another quadrotor as the evader, marked by the blue circle. The flight trajectories of quadrotors in four scenarios are shown in Figure 16, where Figure 16a–d represent the 3D trajectories of quadrotors, while Figure 16e–h represent the projections of the 3D trajectories onto the x-y plane. Notably, the quadrotor successfully completed the pursuit–evasion task while effectively avoiding obstacles. The effectiveness and feasibility of FQL and MPC joint control of quadrotors to carry out PEG were verified.

5. Conclusions

In this paper, we propose a novel control algorithm that combines the FQL and MPC methods for the joint control of quadrotors in a PEG. The proposed algorithm effectively leverages the perception and decision-making capabilities of FQL and the robust control capabilities of MPC. Initially, we used the generalized Apollonius circle in 3D space to determine the agent’s movement strategy. We trained the agent using FQL algorithm and designed a reward function based on the concept of artificial potential field to enhance learning performance. Upon obtaining the evader’s reference trajectory, we designed an MPC algorithm to control quadrotors in the PEG scenario. We developed a state feedback control algorithm for the decoupled quadrotor model, considering constraints on control input and obstacle avoidance. Then, we analyzed the feasibility and stability of the state feedback control law. The simulation results evaluate the control performance of the FQL algorithm and the MPC algorithm separately. Subsequently, utilizing the Gazebo platform, we controlled the quadrotors in the PEG using both the FQL algorithm and the MPC algorithm, and our proposed algorithms exhibited excellent performance. This paper provides valuable insights for combining RL algorithms and traditional control techniques and provides a promising reference for future research. In the future, we will conduct demonstration and validation on a quadrotor platform to compare the performance of the algorithm proposed in this paper with other algorithms. Furthermore, we will also explore the integration of actor–critic structured learning algorithms with MPC algorithms. Additionally, we will also investigate the integration of multi-objective RL algorithms with nonlinear MPC algorithms to address the PEG problem.