Insensitive Mechanism-Based Nonlinear Model Predictive Guidance for UAVs Intercepting Maneuvering Targets with Input Constraints

Abstract

This paper proposed an innovative guidance strategy, denoted as NMPC-IM, which integrates the Insensitive Mechanism (IM) with Nonlinear Model Predictive Control (NMPC) for Unmanned Aerial Vehicle (UAV) pursuit-evasion scenarios, with the aim of effectively intercepting maneuvering targets with consideration of input constraints while minimizing average energy expenditure. Firstly, the basic principle of IM is proposed, and it is transformed into an additional cost function in NMPC. Secondly, in order to estimate the states of maneuvering target, a fixed-time sliding mode disturbance observer is developed. Thirdly, the UAV’s interception task is formulated into a comprehensive Quadratic Programming (QP) problem, and the NMPC-IM guidance strategy is presented, which is then improved by the adjustment of parameters and determination of maximum input. Finally, numerical simulations are carried out to validate the effectiveness of the proposed method, and the simulation results show that the NMPC-IM guidance strategy can decrease average energy expenditure by mitigating the impact of the target’s maneuverability, optimizing the UAV’s trajectory during the interception process.

1. Introduction

With the rapid advancement of Unmanned Aerial Vehicle (UAVs) technology, its role has evolved from singular tracking and reconnaissance to multiple types of mobility tasks [1,2,3]. This includes tasks such as airborne interception and multi-UAV cooperative tracking. This transformation signifies that UAVs are transitioning from traditional reconnaissance platforms to multi-tasking platforms, with the enhancement of their intelligence and autonomy levels becoming a critical factor in the future [4].

The UAV interception task primarily refers to the UAV employing tactical maneuvers to bring the target within its fire range. However, there are constraints of control input, and differences in maneuverability between the UAV and the target. When the UAV is confronted with a more maneuverable target, how it can accomplish the interception task under input constraints while minimizing average energy expenditure is a question that merits investigation. In this paper, it is simplified to a point-to-point guidance problem under conditional constraints.

A number of methods have been proposed to solve the UAV pursuit-evasion game. Ref. [5], considering obstacle avoidance and input constraint, proposed a collaborative strategy combining the Apollonius circle algorithm and the geometric algorithm. Ref. [6] introduced a Reinforcement Learning (RL)-based control strategy for the UAV pursuit-evasion game, ensuring equal formation and stability against disturbances, with a novel Euler-based proof and Hamilton–Jacobi–Isaacs (HJI) equation solution for the Nash equilibrium. Ref. [7] presents a robust swarm strategy for UAVs in urban pursuit-evasion scenarios, integrating predictive decision-making and demonstrating enhanced performance under normal and damaged conditions. Ref. [8] established a mathematical model of UAV pursuit and evasion countermeasures based on the DDPG (deep deterministic policy gradient) algorithm, and designed countermaneuver strategies for evading UAVs. Ref. [9] introduced an autonomous exploration and mapping system utilizing heterogeneous UAVs and UGVs in GPS-denied environments, featuring a two-layered exploration strategy and an optimized next-view planning framework for efficient 3D environment perception and navigation. Ref. [10] proposed a robust adaptive event-triggered control strategy for USV-UAV cooperative formation systems to enhance maritime parallel search efficiency. Ref. [11] proposed a varying cells strategy for UAV collision avoidance that integrates aerodynamic constraints, offering more flexible maneuvers and overcoming the limitations of fixed cell methods.

Model Predictive Control (MPC) utilizes a model of the system to predict its future behavior over a certain horizon, optimizing the control inputs to achieve desired objectives while satisfying constraints. This predictive capability makes MPC particularly attractive for the design of guidance laws, which require precise trajectory planning and control under various dynamic conditions and constraints. This solves the problem of violating constraint of lateral acceleration when using conventional proportional navigation [12].

There have been some applications of MPC to the guidance problem. Ref. [13] proposed a guidance law using Robust MPC optimized by a neural network. Ref. [14] proposed an MPC-based cooperative guidance approach, which launches a salvo attack for multi-missile engagement against a stationary target. Ref. [15] introduced an LPMPC-based suboptimal guidance law for 3D interception of ballistic missiles, and achieved high accuracy in impact angle and miss distance. Similarly, ref. [16] presented an LPMPC cooperative guidance law for multi-missile simultaneous attacks with impact angle constraints. It converts a nonlinear control problem into a solvable Quadratic Programming task. Ref. [17] proposed a hybrid missile guidance algorithm for planar engagement scenarios, which combined Nonlinear MPC with a collision cone approach and ensured that the impact angle lies in a predefined range. More recently, ref. [18] presented an MPC-based guidance system for drone interception with visual tracking and impact angle constraints. It optimizes reduced maneuvering at interception, tested across scenarios with both maneuvering and non-maneuvering targets.

In order to solve the problem of input constraints for the UAV interception task, we incorporate IM into the traditional NMPC guidance model and integrate it into an algorithm for the interception task with some additional operations. The primary contributions of this paper are as follows:

1.

We propose the basic concept of the Insensitive Mechanism (IM), which is expressed by concrete formulations through a UAV interception scenario. This concept is then transformed into an additional cost function within the NMPC framework.

2.

To estimate the target states while reducing the required information acquisition, we design a fixed-time sliding mode disturbance observer to obtain information on the target’s maneuverability.

3.

The UAV interception task is then formulated into a comprehensive Quadratic Programming (QP) problem with the consideration of soft constraints. Consequently, a guidance algorithm that integrates NMPC with IM is presented.

4.

To refine the algorithm, we adjust the parameters and qualitatively analyze the reasons. For the reason that NMPC considers more about the predicted position, we design a determination of maximum input to correct the trajectory in a timely manner.

The structure of the paper is as follows: Section 2 presents the engagement kinematics of UAVs and the basic theories of NMPC in guidance. Section 3 introduces the concept and description of IM, designs a fixed-time sliding mode disturbance observer, and proposes the enhanced NMPC algorithm with IM. In Section 4, the effectiveness and robustness of the proposed NMPC-IM algorithm are validated through simulation, and conclusions are drawn in Section 5.

2. Preliminaries and Problem Statement

To concretely represent the UAV interception, we assume that a UAV intercepts a target. The target is more maneuverable due to its smaller size. The task of the UAV is to minimize the average energy expenditure within the constraints to catch up with the target and destroy it once it reaches the attack range.

As shown in Figure 1, the target may be inside the cone of the UAV at the next moment or outside the cone, the latter meaning that the UAV is constrained by control input and cannot be pointed at the target in time.

In this section, we mainly introduce engagement kinematics of UAVs and NMPC in guidance.

2.1. Engagement Kinematics of UAVs

Assuming that the motion between the UAV and the target is in the same plane, a motion schematic is shown in Figure 2, where Oxy is the inertial frame of reference. Since we consider the UAV and target as particles during the interception and ignore their dynamics, we only need to consider the relative motion of the two points in the local inertial coordinate system.

In this figure, $V_M$ and $V_T$ represent the speed of the UAV and the target, respectively, which remain constant during the engagement. ${\alpha }_M$ and ${\alpha }_T$ represent their heading angles. The quantities of $a_M$ and $a_T$ represent the accelerations of M and T, respectively. r and $\theta$ represent the range and line-of-sight (LOS) angle to the T, respectively. The relationships between the above parameters are represented by the following differential equation:

$$\left\{\begin{array}{cc}\hfill \dot{r}& \hfill =V_r=V_{T}cos\left({\alpha }_T-\theta \right)-V_{M}cos\left({\alpha }_M-\theta \right)\hfill \\ \hfill r\dot{\theta }& \hfill =V_{\theta }=V_{T}sin\left({\alpha }_T-\theta \right)-V_{M}sin\left({\alpha }_M-\theta \right)\hfill \end{array}\right.$$

(1)

In order to establish the state equations describing the motion between M and T, the state vector and control input are defined as $\mathit{x}$ $={\left[\begin{array}{ccccc}r& V_r& \theta & V_{\theta }& {\alpha }_M\end{array}\right]}^T$, $u=a_M=V_M{\dot{\alpha }}_M$, where $V_r$ and $V_{\theta }$ represent the relative velocity components along and perpendicular to the LOS. The kinematics of the engagement can be expressed in the following:

$$\left[\begin{array}{c}\dot{r}\\ {\dot{V}}_r\\ \dot{\theta }\\ {\dot{V}}_{\theta }\\ {\dot{\alpha }}_M\end{array}\right]=\left[\begin{array}{c}{V}_r\\ V_{\theta }^2/r\\ V_{\theta }/r\\ -V_r{V}_{\theta }/r\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ sin\left({\alpha }_M-\theta \right)\\ 0\\ -cos\left({\alpha }_M-\theta \right)\\ 1/V_M\end{array}\right]u+\left[\begin{array}{c}0\\ a_{T}sin\left(\theta -{\alpha }_T\right)\\ 0\\ a_{T}cos\left(\theta -{\alpha }_T\right)\\ 0\end{array}\right]$$

(2)

And the compact form is:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right)+\mathit{g}\left(\mathit{x}\right)u+\mathit{d}$$

(3)

where $u\in \mathbb{R}$, $\mathit{x}\in {\mathbb{R}}^5$, $\mathit{f}:{\mathbb{R}}^5\to {\mathbb{R}}^5$, $\mathit{g}:{\mathbb{R}}^5\to {\mathbb{R}}^5$, $\mathit{d}\in {\mathbb{R}}^5$. The maneuvrability of the target is represented by $\mathit{d}$, which is treated as an unknown bounded disturbance in many articles [17,19]. In this paper, $\mathit{d}$ is considered to be measurable or estimatable, and is used to make improvements to the results.

2.2. Nonlinear Model Predictive Control

MPC typically requires transforming differential equations into discrete form. According to Equation (3), and without the disturbances, the discrete-time form using Euler discretization is expressed in the following:

$$\mathit{x}\left(k+1\right)={\mathit{f}}_d\left(\mathit{x}\left(k\right)\right)+{\mathit{g}}_d\left(\mathit{x}\left(k\right)\right)u\left(k\right)$$

(4)

We use p and c to represent the prediction and control horizons, respectively. The increment in the control input is given by $\mathsf{\Delta }u\left(k\right)=u\left(k\right)-u\left(k-1\right)$. Then, the “referenced predictive form” [20] is

$$\mathit{x}\left(k+j\left|k\right.\right)={\mathit{f}}_d\left(\mathit{x}\left(k+j-1\left|k\right.\right)\right)+{\mathit{g}}_d\left(\mathit{x}\left(k+j-1\left|k\right.\right)\right)\left[u\left(k+j-2\left|k\right.\right)+\mathsf{\Delta }u\left(k+j-1\left|k\right.\right)\right]$$

(5)

where $j=1,\cdots ,p$, and $\left(k+j\left|k\right.\right)$ means that the current time step is k and the distance from the current time step is j. With that, the standard prediction form of Equation (5), required for the NMPC formulations, is given by

$${\mathit{X}}_k={\mathit{F}}_k+{\mathit{g}}_{k}u\left(k-1\left|k-1\right.\right)+{\mathit{G}}_k\mathsf{\Delta }{\mathit{U}}_k$$

(6)

where $u\left(k-1\left|k-1\right.\right)$ is the control input of the last time step, ${\mathit{X}}_k\in {\mathbb{R}}^{5p\times 1}$, ${\mathit{F}}_k\in {\mathbb{R}}^{5p\times 1}$, ${\mathit{g}}_k\in {\mathbb{R}}^{5p\times 1}$, ${\mathit{G}}_k\in {\mathbb{R}}^{5p\times c}$, and $\mathsf{\Delta }{\mathit{U}}_k\in {\mathbb{R}}^{c\times 1}$, $\mathsf{\Delta }{\mathit{U}}_k={\left[\begin{array}{ccc}\mathsf{\Delta }u\left(k\left|k\right.\right)& \cdots & \mathsf{\Delta }u\left(k+c-1\left|k\right.\right)\end{array}\right]}^T$. Since the computation of $\mathit{x}\left(k+j\left|k\right.\right)\left(j=1,\cdots ,p-1\right)$ requires $\mathsf{\Delta }u\left(k+i\left|k\right.\right)\left(i=0,\cdots ,c-1\right)$, which needs to be solved, the predicted states $\mathit{x}\left(k+j\left|k\right.\right)$ cannot be obtained at the current step k, nor can ${\mathit{F}}_k$, ${\mathit{g}}_k$, ${\mathit{G}}_k$. It is a good method to use the predicted states of the previous time step to calculate the above matrics [13]. Therefore, Equation (6) can be reformulated as

$${\mathit{X}}_k={\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}u\left(k-1\left|k-1\right.\right)+{\mathit{G}}_{k-1}\mathsf{\Delta }{\mathit{U}}_k$$

(7)

where ${\mathit{F}}_{k-1}=\left[\begin{array}{c}{\mathit{f}}_d\left(\mathit{x}\left(k-1\left|k-1\right.\right)\right)\\ {\mathit{f}}_d\left(\mathit{x}\left(k\left|k-1\right.\right)\right)\\ ⋮\\ {\mathit{f}}_d\left(\mathit{x}\left(k+c-2\left|k-1\right.\right)\right)\\ ⋮\\ {\mathit{f}}_d\left(\mathit{x}\left(k+p-2\left|k-1\right.\right)\right)\end{array}\right]$, ${\mathit{g}}_{k-1}=\left[\begin{array}{c}{\mathit{g}}_d\left(\mathit{x}\left(k-1\left|k-1\right.\right)\right)\\ {\mathit{g}}_d\left(\mathit{x}\left(k\left|k\right.-1\right)\right)\\ ⋮\\ {\mathit{g}}_d\left(\mathit{x}\left(k+c-2\left|k-1\right.\right)\right)\\ ⋮\\ {\mathit{g}}_d\left(\mathit{x}\left(k+p-2\left|k-1\right.\right)\right)\end{array}\right]$,

${\mathit{G}}_{k-1}=\left[\begin{array}{cccc}{\mathit{g}}_d\left(\mathit{x}\left(k-1\left|k\right.-1\right)\right)& {\mathbf{0}}_{5\times 1}& \cdots & {\mathbf{0}}_{5\times 1}\\ {\mathit{g}}_d\left(\mathit{x}\left(k\left|k-1\right.\right)\right)& {\mathit{g}}_d\left(\mathit{x}\left(k\left|k-1\right.\right)\right)& \cdots & {\mathbf{0}}_{5\times 1}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{g}}_d\left(\mathit{x}\left(k+c-2\left|k-1\right.\right)\right)& {\mathit{g}}_d\left(\mathit{x}\left(k+c-2\left|k-1\right.\right)\right)& \cdots & {\mathit{g}}_d\left(\mathit{x}\left(k+c-2\left|k-1\right.\right)\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{g}}_d\left(\mathit{x}\left(k+p-2\left|k-1\right.\right)\right)& {\mathit{g}}_d\left(\mathit{x}\left(k+p-2\left|k-1\right.\right)\right)& \cdots & {\mathit{g}}_d\left(\mathit{x}\left(k+p-2\left|k-1\right.\right)\right)\end{array}\right].$

In order to design the cost function, it is necessary to define a more appropriate output such as $V_{\theta }$. $V_{\theta }=r\dot{\theta }$, which contains both the position information between the M and T as well as the LOS rate. Therefore, the output is expressed as $y=\mathit{Cx}$, where $\mathit{C}=\left[\begin{array}{ccccc}0& 0& 0& 1& 0\end{array}\right]$. Similarly, $y_d$ denotes the desired output, and usually takes a value like $y_d=c_y\sqrt{r}, c_y⩾0$ [21]. Additionally, the magnitude of the control increments needs to be minimized. Based on the above conditions, the cost function, for time step k, is designed as follows:

$$\begin{array}{c}\hfill J_k={\left({\mathit{Y}}_k-{\mathit{Y}}_{dk}\right)}^T\mathit{Q}\left({\mathit{Y}}_k-{\mathit{Y}}_{dk}\right)+\mathsf{\Delta }{{\mathit{U}}_k}^T\mathit{R}\mathsf{\Delta }{\mathit{U}}_k\end{array}$$

(8)

where ${\mathit{Y}}_k={\left[\begin{array}{ccc}y\left(k+1\left|k\right.\right)& \cdots & y\left(k+p\left|k\right.\right)\end{array}\right]}^T$, ${\mathit{Y}}_{dk}={\left[\begin{array}{ccc}{y}_d\left(k+1\left|k\right.\right)& \cdots & y_d\left(k+p\left|k\right.\right)\end{array}\right]}^T$, $\mathit{Q}=q{\mathit{I}}_p$, $\mathit{R}=r{\mathit{I}}_c$, ${\mathit{Y}}_k=\overline{\mathit{C}}{\mathit{X}}_k=\left({\mathit{I}}_p\otimes \mathit{C}\right){\mathit{X}}_k$. “⊗” denotes the Kronecker product of matrices.

With $q>0$, $r>0$, which represent the weighting factor of output error and control increment, respectively. Therefore, the cost function $J_k$ in Equation (8) can be expressed as

$$\begin{array}{ll}{J}_k& ={\left(\overline{\mathit{C}}{\mathit{X}}_k-{\mathit{Y}}_{dk}\right)}^T\mathit{Q}\left(\overline{\mathit{C}}{\mathit{X}}_k-{\mathit{Y}}_{dk}\right)+\mathsf{\Delta }{{\mathit{U}}_k}^T\mathit{R}\mathsf{\Delta }{\mathit{U}}_k\\ & =J_{ck}+\mathsf{\Delta }{{\mathit{U}}_k}^T{\mathit{H}}_k\mathsf{\Delta }{\mathit{U}}_k+{\mathit{f}}_k^T\mathsf{\Delta }{\mathit{U}}_k\end{array}$$

(9)

where ${\mathit{H}}_k={\mathit{G}}_{k-1}^T{\overline{\mathit{C}}}^T\mathit{Q}\overline{\mathit{C}}{\mathit{G}}_{k-1}+\mathit{R}$, ${\mathit{f}}_k^T=2{\left(\overline{\mathit{C}}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)-{\mathit{Y}}_{dk}\right)}^T\mathit{Q}\overline{\mathit{C}}{\mathit{G}}_{k-1}$, $J_{ck} =$${\left(\overline{\mathit{C}}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)-{\mathit{Y}}_{dk}\right)}^T\mathit{Q}\left(\overline{\mathit{C}}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)-{\mathit{Y}}_{dk}\right)$. Note that the later optimization problem is to find $\mathsf{\Delta }{\mathit{U}}_k$ to minimize $J_k$, while the change of $\mathsf{\Delta }{\mathit{U}}_k$ does not affect $J_{ck}$, which can be disregarded.

Like most MPC problems, some constraints on the amount of control and control increments are required. They are constrained as

$${\mathit{U}}_{min}⩽{\mathit{U}}_k⩽{\mathit{U}}_{max}$$

(10)

$${\mathit{U}}_k={\mathit{U}}_{k-1}+{\mathit{I}}_{1t}\mathsf{\Delta }{\mathit{U}}_k$$

(11)

$$\mathsf{\Delta }{\mathit{U}}_{min}⩽\mathsf{\Delta }{\mathit{U}}_k⩽\mathsf{\Delta }{\mathit{U}}_{max}$$

(12)

where ${\mathit{U}}_{k-1}=u\left(k-1\left|k-1\right.\right){\mathit{1}}_c$, ${\mathit{I}}_{1t}={\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \cdots & 1\end{array}\right]}_{c\times c}$.

Additionally, $V_r$ is required to remain negative for a successful interception, for the reason that $V_r$ represents the relative velocity components along the LOS. Only when $V_r<0$ does the distance between the UAV and the target gradually decrease so that interception can be achieved [17]. The terminal state in the prediction is constrained to be negative, for instance, less than a negative number $v_{rd}$. We choose $v_{rd}=-sign\left(V_{r0}\right) \times 0.1V_{r0}$ here, where $V_{r0}$ represents the initial value of $V_r$. The terminal constraint is given by

$${\mathit{EX}}_k⩽{\mathit{V}}_{rd}$$

(13)

where $\mathit{E}={\left[\begin{array}{ccc}{\mathbf{0}}_{1\times 5}& \cdots & {\mathbf{0}}_{1\times 5}\\ ⋮& \ddots & ⋮\\ {\mathbf{0}}_{1\times 5}& \cdots & \mathit{e}\end{array}\right]}_{p\times 5p}$, ${\mathit{V}}_{rd}={\left[\begin{array}{c}0\\ ⋮\\ 0\\ v_{rd}\end{array}\right]}_{p\times 1}$, $\mathit{e}=\left[\begin{array}{ccccc}0& 1& 0& 0& 0\end{array}\right]$. Substituting ${\mathit{X}}_k$ from Equation (7) in Equation (13), we obtain

$${\mathit{EG}}_{k-1}\mathsf{\Delta }{\mathit{U}}_k⩽{\mathit{V}}_{rd}-\mathit{E}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)$$

(14)

Considering the above constraints, a compact form can be expressed as

$${\mathit{A}}_k\mathsf{\Delta }{\mathit{U}}_k⩽{\mathit{b}}_k$$

(15)

where ${\mathit{A}}_k=\left[\begin{array}{c}{\mathit{I}}_c\\ -{\mathit{I}}_c\\ {\mathit{I}}_{1t}\\ -{\mathit{I}}_{1t}\\ {\mathit{EG}}_{k-1}\end{array}\right]$, ${\mathit{b}}_k=\left[\begin{array}{c}\mathsf{\Delta }{\mathit{U}}_{max}\\ -\mathsf{\Delta }{\mathit{U}}_{min}\\ {\mathit{U}}_{max}-{\mathit{U}}_{k-1}\\ -{\mathit{U}}_{min}+{\mathit{U}}_{k-1}\\ {\mathit{V}}_{rd}-\mathit{E}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)\end{array}\right]$. With the above formulations of the cost function and constraints, the NMPC in guidance can be described as a QP problem:

$$\begin{array}{c}{\displaystyle \underset{\mathsf{\Delta }{\mathit{U}}_k}{min}}\mathsf{\Delta }{{\mathit{U}}_k}^T{\mathit{H}}_k\mathsf{\Delta }{\mathit{U}}_k+{\mathit{f}}_k^T\mathsf{\Delta }{\mathit{U}}_k\\ s.t.{\mathit{A}}_k\mathsf{\Delta }{\mathit{U}}_k⩽{\mathit{b}}_k\end{array}$$

(16)

Assuming that the solution to the QP in Equation (16) at time step k is given by $\mathsf{\Delta }{\mathit{U}}_k^{\ast }$, then ${\mathit{U}}_k$ can be calculated by Equation (11). The first element of ${\mathit{U}}_k$ is $u\left(k\left|k\right.\right)$, which is used as the control input at time step k. Using ${\mathit{U}}_k$ and the measurement of state at the current time step $\mathit{x}\left(k\left|k\right.\right)$, we can compute ${\mathit{X}}_s\left(k\right)={\left[\begin{array}{cccc}{\mathit{x}}^T\left(k\left|k\right.\right)& {\mathit{x}}^T\left(k+1\left|k\right.\right)& \cdots & {\mathit{x}}^T\left(k+p-1\left|k\right.\right)\end{array}\right]}^T$ by Equation (5). ${\mathit{X}}_s\left(k\right)$ is then used at the next time step to calculate the matrices ${\mathit{F}}_{k-1}$, ${\mathit{g}}_{k-1}$, ${\mathit{G}}_{k-1}$. Further, ${\mathit{U}}_k$ is reformulated as ${\mathit{U}}_k=u\left(k\left|k\right.\right){\mathbf{1}}_c$. After the above algorithm, a closed loop control is accomplished.

Remark 1.

At the begining step of the simulation, $u\left(k-1\left|k-1\right.\right)=0$, ${\mathit{U}}_{k-1}={\mathbf{0}}_{c\times 1}$ and ${\mathit{X}}_s\left(k-1\right)={\mathbf{1}}_p\otimes {\mathit{x}}_0$, where ${\mathit{x}}_0$ is the initial state.

2.3. Problem Statement

In order to guide the UAV under input constraints, we designed the basic guidance law by transforming the guidance conditions into a cost function and the constraints into inequality constraints using NMPC. Furthermore, it is our intention to attenuate the oscillatory variation of the control inputs throughout the guidance in order to reduce the impact of the target maneuvering, so that the average energy expenditure is as small as possible. In order to quantify this, we set two parameters as follows:

$$Su={\int }_0^{t_f}\left|u\right|dt$$

(17)

$$su=Su/t_f$$

(18)

where $t_f$ represents the impact time and also the time the simulation stops.

$Su$ in Equation (17) represents accumulated control losses, $su$ in Equation (18) represents average energy expenditure during guidance.

3. Insensitive Mechanism-Based Nonlinear Model Predictive Guidance

In this section, we first introduce the Insensitive Mechanism (IM) in NMPC, then design a fixed-time sliding mode disturbance observer, and finally elaborate a combination algorithm of NMPC and IM used for guidance. In order to make the algorithm more complete, the parameters are adjusted and a determination condition is added as to whether the maximum input is required or not. The general structure of what is designed in this section is shown in Figure 3.

M and T provide the initial position information, and the information related to LOS is calculated based on Equation (1). The maneuvering states of target are obtained after the disturbance observer. The control inputs to the UAV are calculated by the NMPC-IM algorithm. Meanwhile, Equation (41) is utilized for Maximum Input Determination, and if the conditions are satisfied, the input calculated by NMPC-IM will be substituted.

3.1. Design of Insensitive Mechanism

When using NMPC, as well as other methods of guidance on maneuvering targets, if T is bending forward, M may behave sensitively to the movement of T and thus make unnecessary steering control. Especially when M is far away from T, the steering control may produce wasteful lateral overload. Due to this, the guidance method may need to be less sensitive.

One simple idea is that if we can measure or estimate the maneuvering state of T, we might be able to make a prediction of its position some time in the future, choosing to strike directly toward the predicted position without or with little regard for the current position of T [22,23].

A schema of the basic idea of the IM is shown in Figure 4. $M_k$ and $T_k$ represent the position of M and T at time step k, respectively. After a p-step prediction, the position of T is $T_{k+p}$. The blue dotted line represents the predicted trajectory of T, the red solid line is M’s guidance to the current position of T, and the green solid line is M’s guidance to the predicted position of T. In the IM, we would rather choose the green one than choose the red one.

We first assuming that the maneuvering state of T is known to be $a_T$, the position and the heading angle of the target can be measured. In order to represent the insensitivity of M to the maneuverability of T, it is assumed that M remains stationary while the position of T is predicted. According to Equation (1), the kinematic relationship between M and T is simplified as

$$\left\{\begin{array}{cc}\hfill \dot{r}& \hfill =V_{T}cos\left({\alpha }_T-\theta \right)\hfill \\ \hfill r\dot{\theta }& \hfill =V_{T}sin\left({\alpha }_T-\theta \right)\hfill \\ \hfill {\dot{\alpha }}_T& =a_T/V_T\hfill \end{array}\right.$$

(19)

In fact, Equation (19) is not a description of the objective world motion, but a hypothesis required at the moment of calculating the control inputs to the UAV. In this instant, the UAV assumes that it remains stationary and predicts the position of the target in the future, which allows the UAV not to be limited to the target’s current true position, but to reduce unnecessary inputs by anticipating the future. Equation (1), on the other hand, is an objective description of the relative motion of objects that exists in the objective world and follows objective laws. Therefore, it can be considered that Equation (19) is a subjective hypothesis required for prediction and is where the basic idea of IM lies, and it does not conflict with Equation (1).

At the time step k, the predicted position of T after p steps is calculated by Equation (19). The step-by-step calculation of the Insensitive Mechanism is shown in Figure 5.

The red arrows represent the positional shift of the target at each step during the prediction time period, while the missile remains stationary. For $j=1,\cdots ,p$, $r_s\left(k+j\left|k\right.\right)$, and ${\theta }_s\left(k+j\left|k\right.\right)$ represent the range and LOS angle from $M_k$ to the $T_{k+j}$, respectively. Similarly, the heading angle of T can be expressed as ${\alpha }_{Ts}\left(k+j\left|k\right.\right)$, which is given by linear prediction with known $a_T$. Using Equation (19) and the measurement of M and T’s state at time step k, we can compute all the above parameters in the prediction time domain step by step. Once the position of $T_{k+p}$ is obtained, all the M needs to do is to keep the heading angle ${\theta }_s\left(k+p\left|k\right.\right)$, which is abbreviated as ${\theta }_{sd}$.

However, in actual motion, M does not remain stationary. Here, we design an ideal trajectory based on (1) and ${\theta }_{sd}$. Before this, a new output is expressed as ${\mathit{y}}_s={\mathit{C}}_s\mathit{x}$, where ${\mathit{C}}_s=\left[\begin{array}{c}\begin{array}{ccccccccccc}1& 0& 0& 0& 0& ;& 0& 0& k_{\theta }& 0& 0\end{array}\end{array}\right]$. In this case, ${\mathit{y}}_s={\left[\begin{array}{cc}r& k_{\theta }\theta \end{array}\right]}^T$. Similarly, ${\mathit{y}}_{sd}$ denotes the desired output, where ${\mathit{y}}_{sd}\left(k+j\left|k\right.\right)={\left[\begin{array}{cc}{r}_{sd}\left(k+j\left|k\right.\right)& k_{\theta }{\theta }_{sd}\left(k+j\left|k\right.\right)\end{array}\right]}^T$. Generally speaking, the range of the variables can be used as a reference for the value of $k_{\theta }$. For instance, r ranges from 0 to 200, and $\theta$ ranges from 0 to 6.28. In order to balance the range of variation of the two variables, $k_{\theta }$ is supposed to be 32. However, it is unlikely that the two variables will change in the same way, i.e., error values between expected and actual may not even be of the same order of magnitude in each step. Since we design the algorithm for ${\theta }_{sd}$ in IM, we expect its corresponding error value to be smaller. In order for this error term in the QP problem to have a larger impact on the results, a larger $k_{\theta }$ is required. After testing the number in different orders of magnitude, we choose $k_{\theta }=5000$. The calculations of $r_{sd}\left(k+j\left|k\right.\right)$ and ${\theta }_{sd}\left(k+j\left|k\right.\right)$ are displayed as follows

$$\begin{array}{l}forj=1\\ r_{sd}\left(k+1\left|k\right.\right)=r\left(k\right)+\mathsf{\Delta }t{V}_{T}cos\left({\alpha }_T\left(k\right)-\theta \left(k\right)\right)-\mathsf{\Delta }t{V}_{M}cos\left({\theta }_{sd}-\theta \left(k\right)\right)\\ {\theta }_{sd}\left(k+1\left|k\right.\right)=\theta \left(k\right)+\left(\mathsf{\Delta }t/r\left(k\right)\right)V_{T}sin\left({\alpha }_T\left(k\right)-\theta \left(k\right)\right)-\left(\mathsf{\Delta }t/r\left(k\right)\right)V_{M}sin\left({\theta }_{sd}-\theta \left(k\right)\right)\\ forj=2,\cdots ,p\\ r_{sd}\left(k+j\left|k\right.\right)=r_{sd}\left(k+j-1\left|k\right.\right)+\mathsf{\Delta }t{V}_{T}cos\left({\alpha }_{Ts}\left(k+j-1\left|k\right.\right)-{\theta }_{sd}\left(k+j-1\left|k\right.\right)\right)\\ -\mathsf{\Delta }t{V}_{M}cos\left({\theta }_{sd}-{\theta }_{sd}\left(k+j-1\left|k\right.\right)\right)\\ {\theta }_{sd}\left(k+j\left|k\right.\right)={\theta }_{sd}\left(k+j-1\left|k\right.\right)\\ +\left(\mathsf{\Delta }t/r_{sd}\left(k+j-1\left|k\right.\right)\right)V_{T}sin\left({\alpha }_{Ts}\left(k+j-1\left|k\right.\right)-{\theta }_{sd}\left(k+j-1\left|k\right.\right)\right)\\ -\left(\mathsf{\Delta }t/r_{sd}\left(k+j-1\left|k\right.\right)\right)V_{M}sin\left({\theta }_{sd}-{\theta }_{sd}\left(k+j-1\left|k\right.\right)\right)\end{array}$$

(20)

Similar to Equation (9), the cost function, designed for IM, at time step k, is as follows:

$$\begin{array}{ll}{J}_{sk}& ={\left({\mathit{Y}}_{sk}-{\mathit{Y}}_{sd}\right)}^T\mathit{S}\left({\mathit{Y}}_{sk}-{\mathit{Y}}_{sd}\right)\\ & =J_{sck}+\mathsf{\Delta }{{\mathit{U}}_k}^T{\mathit{H}}_{sk}\mathsf{\Delta }{\mathit{U}}_k+{\mathit{f}}_{sk}^T\mathsf{\Delta }{\mathit{U}}_k\end{array}$$

(21)

where ${\mathit{Y}}_{sk}={\left[\begin{array}{ccc}{\mathit{y}}_s^T\left(k+1\left|k\right.\right)& \cdots & {\mathit{y}}_s^T\left(k+p\left|k\right.\right)\end{array}\right]}^T$, ${\mathit{Y}}_{sd}={\left[\begin{array}{ccc}{\mathit{y}}_{sd}^T\left(k+1\left|k\right.\right)& \cdots & {\mathit{y}}_{sd}^T\left(k+p\left|k\right.\right)\end{array}\right]}^T$, $\mathit{S}=s{\mathit{I}}_{2p}$, ${\mathit{Y}}_{sk}=\overline{{\mathit{C}}_s}{\mathit{X}}_k=\left({\mathit{I}}_p\otimes {\mathit{C}}_s\right){\mathit{X}}_k$.

${\mathit{H}}_{sk}={\mathit{G}}_{k-1}^T{\overline{{\mathit{C}}_{\mathit{s}}}}^T\mathit{S}\overline{{\mathit{C}}_{\mathit{s}}}{\mathit{G}}_{k-1}$, ${\mathit{f}}_{sk}^T=2{\left(\overline{{\mathit{C}}_{\mathit{s}}}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)-{\mathit{Y}}_{sd}\right)}^T\mathit{S}\overline{{\mathit{C}}_{\mathit{s}}}{\mathit{G}}_{k-1}$,

$J_{sck}={\left(\overline{{\mathit{C}}_{\mathit{s}}}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)-{\mathit{Y}}_{sd}\right)}^T\mathit{S}\left(\overline{{\mathit{C}}_{\mathit{s}}}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)-{\mathit{Y}}_{sd}\right)$

with $s>0$, which represents the weighting factor of the output error in IM.

It is noteworthy that Equations (9) and (21) have the same structure, which allows the addition of IM to NMPC to produce certain effects

Remark 2.

Although the target maneuvering state assumed here is known, it can still be set to a condition that cannot be measured, which could be estimated by using the disturbance observer later. Further, a time series prediction of $a_T$ can be used, for example, a neural network. The more accurate the prediction of the target’s movement, the more effective IM is.

Moreover, the concept of IM is not limited to this paper. Scenarios such as a robotic arm gripping a flexible moving object under constraints, or autonomous driving with more soothing control inputs for the comfort of the passengers, are examples of scenarios where IM has a promising application.

3.2. NMPC-IM Guidance Law with Fixed-Time Sliding Mode Observer

In the previous sections, we derived the prediction equation, cost function, and control constraints in NMPC and then transformed the optimization problem into a QP problem. In addition, the conditions of IM are transformed into another cost function to be added into NMPC. We combine the two into a QP problem and add some additional conditions to the composition of the algorithm. To address the case where the maneuvering state of the target is unknown, we introduce a fixed-time sliding mode observer to estimate the state of the target, and eventually give a complete proof of feasibility.

3.2.1. Full Description of the QP Problem

In many cases, when the MPC problem is transformed into a QP problem, it is usually necessary to use a slack variable to solve the unsolvable problem due to hard constraints. To soften the terminal constraint for $V_r$ in Equation (14), we relax the inequality with

$${\mathit{EG}}_{k-1}\mathsf{\Delta }{\mathit{U}}_k+{\left[\begin{array}{cc}{\mathbf{0}}_{1\times \left(p-1\right)}& -{\gamma }_k\end{array}\right]}^T⩽{\mathit{V}}_{rd}-\mathit{E}\left({\mathit{F}}_{k-1}+{\mathit{g}}_{k-1}\right)$$

(22)

where ${\gamma }_k⩾0$ is the slack variable. With the addition of ${\gamma }_k$, and the cost function in IM, the QP in Equation (16) is modified as

$$\begin{array}{c}{\displaystyle \underset{{\mathit{\zeta }}_k}{min}}{\mathit{\zeta }}_k^T{\tilde{\mathit{H}}}_k{\mathit{\zeta }}_k+{\tilde{\mathit{f}}}_k^T{\mathit{\zeta }}_k\\ s.t.{\tilde{\mathit{A}}}_k{\mathit{\zeta }}_k⩽{\tilde{\mathit{b}}}_k\end{array}$$

(23)

where ${\mathit{\zeta }}_k={\left[\begin{array}{cc}\mathsf{\Delta }{\mathit{U}}_k^T& {\gamma }_k\end{array}\right]}^T\in {\mathbb{R}}^{c+1}$, ${\tilde{\mathit{H}}}_k=\mathrm{diag}\left({\mathit{H}}_k+{\mathit{H}}_{sk},r_1\right)$, ${\tilde{\mathit{f}}}_k=\left[\begin{array}{c}{\mathit{f}}_k+{\mathit{f}}_{sk}\\ r_2\end{array}\right]$, ${\tilde{\mathit{A}}}_k={\left[\begin{array}{cc}{\mathit{A}}_k& \tilde{\mathbf{0}}\\ {\mathbf{0}}_{1\times c}& -1\end{array}\right]}_{\left(4c+p+1\right)\times \left(c+1\right)}$, ${\tilde{\mathit{b}}}_k={\left[\begin{array}{c}{\mathit{b}}_k\\ 0\end{array}\right]}_{\left(4c+p+1\right)\times 1}$.

with $r_1,r_1>0$, $\tilde{\mathbf{0}}={\left[\begin{array}{c}{\mathbf{0}}_{c\times 1}\\ {\mathbf{0}}_{c\times 1}\\ {\mathbf{0}}_{c\times 1}\\ {\mathbf{0}}_{c\times 1}\\ \overline{\mathbf{0}}\end{array}\right]}_{\left(4c+p\right)\times 1}$, $\overline{\mathbf{0}}=\left[\begin{array}{c}{\mathbf{0}}_{\left(p-1\right)\times 1}\\ -1\end{array}\right]$.

The combination of NMPC and IM belongs to a linear superposition; however, with a fixed linear relationship, it is certainly more difficult to fulfill the demand. For example, when the maneuverability of T is relatively large, we would like to consider IM more.

We choose ${\dot{a}}_T$ to represent the maneuverability of T; the corresponding weighting factors in the cost function of NMPC and IM are adjusted, respectively.

$$\left\{\begin{array}{l}q=0.1\times (1-{\displaystyle \frac{2}{\pi }}\mathrm{arc}tan\left(\left|{\dot{a}}_T\right|\right))\\ r=1\\ s=10\times \left(5+10\left|{\dot{a}}_T\right|\right)\end{array}\right.$$

(24)

According to Equation (24), it can be seen that q and s vary monotonically decreasing and monotonically increasing with respect to $\left|{\dot{a}}_T\right|$, respectively. Additionally, the effect of target maneuvering on the result of prediction is considered. In the IM, the predicted target is considered to have constant $\left|{\dot{a}}_T\right|$, which means the result of prediction could be better. In contrast, when $\left|{\dot{a}}_T\right|$ changes greatly, the result could be worse, and we hope that p could be smaller in order to avoid a larger error. Therefore, p should be inversely proportional to $\left|{\dot{a}}_T\right|$. In addition, the NMPC-based guidance law may focus on the predicted position of T, resulting in an excessive amount of eventual off-targeting. Thus, once $r<r_0$, a threshold, we expect p to decrease as r decreases, which can expressed as follows:

$$p=\left\{\begin{array}{c}\begin{array}{cc}\left[\left(1-{\displaystyle \frac{2}{\pi }}\mathrm{arc}tan\left(\left|{\dot{a}}_T\right|\right)\right)\left(p_{max}-p_{min}\right)+p_{min}\right],& r>r_0\end{array}\\ \begin{array}{cc}\left[{\displaystyle \frac{r}{r_0}}\left(1-{\displaystyle \frac{2}{\pi }}\mathrm{arc}tan\left(\left|{\dot{a}}_T\right|\right)\right)\left(p_{max}-p_{min}\right)+p_{min}\right],& r⩽r_0\end{array}\end{array}\right.$$

(25)

where “$[\ast ]$” represents the Gauss mark; $p_{max}>p_{min}>0$ are positive integers.

Remark 3.

Because p is changeable at each step, ${\mathit{X}}_s\left(k-1\right)\in {\mathbb{R}}^{5p\times 1}$ changes too. While calculating ${\mathit{F}}_{k-1}$, ${\mathit{g}}_{k-1}$, ${\mathit{G}}_{k-1}$ using ${\mathit{X}}_s\left(k-1\right)$, we need to change the scale of ${\mathit{X}}_s\left(k-1\right)$ first. Actually, ${\mathit{X}}_s\left(k-1\right)\in {\mathbb{R}}^{5p_{k-1}\times 1}$, where $p_{k-1}$ represents p the previous time step. After computing $p_k$ at this time step by Equation (25), change ${\mathit{X}}_s\left(k-1\right)$ as

$${\mathit{X}}_s\left(k-1\right)=\left\{\begin{array}{c}{\left[\begin{array}{c}\mathit{x}\left(k-1\left|k-1\right.\right)\\ ⋮\\ \mathit{x}\left(k+p_{k-1}-2\left|k-1\right.\right)\\ ⋮\\ \mathit{x}\left(k+p_{k-1}-2\left|k-1\right.\right)\end{array}\right]}_{5p_k\times 1},p_{k-1}<p_k\\ {\left[\begin{array}{c}\mathit{x}\left(k-1\left|k-1\right.\right)\\ \mathit{x}\left(k\left|k-1\right.\right)\\ ⋮\\ \mathit{x}\left(k+p_k-2\left|k-1\right.\right)\end{array}\right]}_{5p_k\times 1},p_{k-1}⩾p_k\end{array}\right.$$

(26)

After adjusting the weighting factors and the prediction time domain, the IM in NMPC can work better.

This is a complete description of the QP problem in NMPC-IM.

3.2.2. Fixed-Time Sliding Mode Observer for Acceleration Estimation of Maneuvering Target

Sliding Mode Observer (SMO) is a type of nonlinear observer commonly used for estimating the state of a system. It operates by constructing an observation error system and designing a switching control law achieving an accurate estimation of the system’s state. There are also successful examples of combining observer and MPC [24,25].

To ensure that the state of the system converges to the origin in a finite time T and the convergence time is minimally affected by the initial state, ref. [26] proposed a control law by improving the super-twisting algorithm.

Assuming that the maneuvering states of T are unknown and bounded. We let $a_{Tr}=a_{T}sin(\theta -{\alpha }_T)$ and $a_{T\theta }=a_{T}cos(\theta -{\alpha }_T)$ represent the maneuvering states of T, and $\left|a_{Tr}\right|, \left|a_{T\theta }\right|⩽L$, $\left|{\dot{a}}_{Tr}\right|,\left|{\dot{a}}_{T\theta }\right|⩽H$. According to Equation (2), there exist two first-order systems: ${\dot{V}}_r={\overline{u}}_1+a_{Tr}$, where ${\overline{u}}_1=V_{\theta }^2/r+sin\left({\alpha }_M-\theta \right)u$ and ${\dot{V}}_{\theta }={\overline{u}}_2+a_{T\theta }$, where ${\overline{u}}_2=-V_r{V}_{\theta }/r-cos\left({\alpha }_M-\theta \right)u$.

With the addition of the designed fixed-time sliding mode disturbance observer, a theorem about the feasibility of the QP in Equation (23) is given as follows:

Theorem 1.

Observer (27) can observe the actual value of the maneuvering states $a_{Tr}$ and $a_{T\theta }$, respectively. The errors converge to zero in a fixed time $T_f$, which is independent of the initial states. The QP in Equation (23) is a Strictly Convex QP (SCQP) for all k, and the errors will not effect the feasibility of the solution.

$$\left\{\begin{array}{ll}\dot{x}& =u+d\\ {\sigma }_1& =z_1-x\\ {\sigma }_2& =z_2-d\\ {\dot{z}}_1& =-k_1{\left|{\sigma }_1\right|}^{p}sign\left({\sigma }_1\right)-k_2{\left|{\sigma }_1\right|}^{q}sign\left({\sigma }_1\right)+z_2+u\\ {\dot{z}}_2& =-k_{3}sign\left({\sigma }_1\right)\end{array}\right.$$

(27)

$$T_f⩽M\left({\displaystyle \frac{1}{k_1\left(p-1\right){\epsilon }^{p-1}}}+{\displaystyle \frac{{\epsilon }^{1-q}}{k_2\left(1-q\right)}}\right){\displaystyle \frac{1}{m\left(1-Mh\left(k_2\right)/k_2\right)}}$$

(28)

where $k_1,k_2>0$, $k_3>H$, $p>1$, $0<q<1$, $\epsilon >0$, $h\left(k_2\right)=1/k_2+{\left(2e/m{k}_2\right)}^{1/3}$, $M=k_3+H$, $m=k_3-H$.

Proof of Theorem 1.

We use the estimation of $a_{Tr}$ as an example. For the purpose of the proof, the initial time of the system is considered to be 0 and the value of the sliding mode surface at the initial moment is ${\sigma }_1\left(0\right)$, and ${\sigma }_2\left(0\right)=0$ [27]. The time derivatives of ${\sigma }_1$ and ${\sigma }_2$ are

$$\left\{\begin{array}{l}{\dot{\sigma }}_1={\dot{z}}_1-{\dot{V}}_r=-k_1{\left|{\sigma }_1\right|}^{p}sign\left({\sigma }_1\right)-k_2{\left|{\sigma }_1\right|}^{q}sign\left({\sigma }_1\right)+z_2-a_{Tr}\\ {\dot{\sigma }}_2={\dot{z}}_2-{\dot{a}}_{Tr}=-k_{3}sign\left({\sigma }_1\right)-{\dot{a}}_{Tr}\end{array}\right.$$

(29)

Given $\epsilon >0$, when $\left|{\sigma }_1\left(0\right)\right|>\epsilon$, according to Equation (29), we have

$${\displaystyle \frac{d\left|{\sigma }_1\left(t\right)\right|}{dt}}⩽-k_1{\left|{\sigma }_1\left(t\right)\right|}^p$$

(30)

It is easy to analyze that when $t>0$, the signs of ${\sigma }_1$ and ${\sigma }_2$ are opposite, i.e., the signs of ${\sigma }_1$ remain unchanged, so Equation (30) can be expressed as

$${\left|{\sigma }_1\left(t\right)\right|}^{-p}d\left|{\sigma }_1\left(t\right)\right|⩽-k_{1}dt$$

(31)

Integrating both sides of the inequality yields:

$${\displaystyle \frac{{\left|{\sigma }_1\left(t\right)\right|}^{1-p}}{1-p}}-{\displaystyle \frac{{\left|{\sigma }_1\left(0\right)\right|}^{1-p}}{1-p}}⩽-k_1\left(t-0\right)$$

(32)

$${\displaystyle \frac{{\left|{\sigma }_1\left(t\right)\right|}^{1-p}}{1-p}}⩽-k_1\left(t-0\right)+{\displaystyle \frac{{\left|{\sigma }_1\left(0\right)\right|}^{1-p}}{1-p}}⩽-k_{1}t$$

(33)

Multiply both sides simultaneously by $(p-1)$:

$$-{\left|{\sigma }_1\left(t\right)\right|}^{1-p}⩽-k_1\left(p-1\right)t$$

(34)

$${\left|{\sigma }_1\left(t\right)\right|}^{p-1}⩽{\displaystyle \frac{1}{k_1\left(p-1\right)t}}$$

(35)

After analysis, it can be obtained that $\left|{\sigma }_1\left(t\right)\right|$ gradually decreases to $\epsilon$ in $\left(0,T_1\right)$, where $T_1⩽1/\left(k_1\left(p-1\right){\epsilon }^{p-1}\right)$. The signs of ${\sigma }_2$ and ${\sigma }_1$ are opposite and $k_3>H$, $\left|{\sigma }_2\left(t\right)\right|$ increases in $\left(0,T_1\right)$.

In the time period $t>T_1$, $\left|{\sigma }_1\left(t\right)\right|$ decreases from $\epsilon$ to 0. According to Equation (29), we have

$${\displaystyle \frac{d\left|{\sigma }_1\left(t\right)\right|}{dt}}⩽-k_2{\left|{\sigma }_1\left(t\right)\right|}^q$$

(36)

Solving this inequality gives

$${\left|{\sigma }_1\left(t\right)\right|}^{1-q}⩽-k_2\left(1-q\right)\left(t-T_1\right)+{\epsilon }^{1-q}$$

(37)

In this case, $\left|{\sigma }_1\left(t\right)\right|$ gradually decreases to 0 in $\left(T_1,T_2\right)$, where $T_2⩽{\epsilon }^{1-q}/\left(k_2\left(1-q\right)\right)$. $\left|{\sigma }_2\left(t\right)\right|$ increases in $\left(T_1,T_2\right)$, and $\left|{\sigma }_2\left(T_2\right)\right|<M{T}_2=M\left({\displaystyle \frac{1}{k_1\left(p-1\right){\epsilon }^{p-1}}}+{\displaystyle \frac{{\epsilon }^{1-q}}{k_2\left(1-q\right)}}\right)$.

In the time period $t>T_2$, the trajectory of the system state starts from $\left(0,{\sigma }_2\left(T_2\right)\right)$ and gradually converges to the origin. Since the parameters of the observer satisfy the inequality $M\left(h\left(k_2\right)/k_2\right)<1$, it can be obtained that the convergence time satisfies the inequality $T_f⩽{\sigma }_2\left(T_2\right)/\left(m\left(1-Mh\left(k_2\right)/k_2\right)\right)$ according to the proof in ref. [26].

In summary, the sliding mode surfaces of the observer ${\sigma }_1\left(t\right)$ and ${\sigma }_2\left(t\right)$ can converge at a fixed time $T_f$, where $T_f$ satisfies

$$T_f⩽M\left({\displaystyle \frac{1}{k_1\left(p-1\right){\epsilon }^{p-1}}}+{\displaystyle \frac{{\epsilon }^{1-q}}{k_2\left(1-q\right)}}\right){\displaystyle \frac{1}{m\left(1-Mh\left(k_2\right)/k_2\right)}}$$

(38)

We use ${\widehat{a}}_{Tr}$ and ${\widehat{a}}_{T\theta }$ to represent the estimated maneuvering states, and they will converge to the real value in $T_f$ from the above proof.

The feasible set for the optimization problem in Equation (23) is a polytope for all k. Hence, the feasible set is convex for all k [28]. The Hessian matrix ${\mathit{H}}_k+{\mathit{H}}_{sk}$ can be rewritten as ${\mathit{H}}_k+{\mathit{H}}_{sk}={\mathit{G}}_{k-1}^T{\overline{\mathit{C}}}^T\mathit{Q}\overline{\mathit{C}}{\mathit{G}}_{k-1}+\mathit{R}+{\mathit{G}}_{k-1}^T{\overline{{\mathit{C}}_{\mathit{s}}}}^T\mathit{S}\overline{{\mathit{C}}_{\mathit{s}}}{\mathit{G}}_{k-1}$. The maneuvering states have been abandoned in Equation (4), so ${\mathit{G}}_{k-1}$ will not contain ${\widehat{a}}_{Tr}$ or ${\widehat{a}}_{T\theta }$ according to Equation (7). Obviously, $\overline{\mathit{C}}$ and $\overline{{\mathit{C}}_{\mathit{s}}}$ are also not related to ${\widehat{a}}_{Tr}$, ${\widehat{a}}_{T\theta }$. As for the parameters q, r, s in Equation (24), $a_T$ and $\left|{\dot{a}}_T\right|$ are given by

$$a_T^2={\widehat{a}}_{Tr}^2+{\widehat{a}}_{T\theta }^2$$

(39)

$$\left|{\dot{a}}_T\right|={\displaystyle \frac{\left|{\widehat{a}}_{Tr}{\dot{\widehat{a}}}_{Tr}+{\widehat{a}}_{T\theta }{\dot{\widehat{a}}}_{T\theta }\right|}{\sqrt{{\widehat{a}}_{Tr}^2+{\widehat{a}}_{T\theta }^2}}}$$

(40)

Clearly, q, r, s remain larger than zero, the corresponding $\mathit{Q}$, $\mathit{R}$, $\mathit{S}$ are positive definite. Therefore, ${\mathit{G}}_{k-1}^T{\overline{\mathit{C}}}^T\mathit{Q}\overline{\mathit{C}}{\mathit{G}}_{k-1}$ and ${\mathit{G}}_{k-1}^T{\overline{{\mathit{C}}_{\mathit{s}}}}^T\mathit{S}\overline{{\mathit{C}}_{\mathit{s}}}{\mathit{G}}_{k-1}$ are positive semidefinite (at least). Thus, ${\mathit{H}}_k+{\mathit{H}}_{sk}$ is positive definite and the cost function $J_k+J_{sk}$ is strictly convex for all k. In this way, due to the positive definiteness of ${\mathit{H}}_k+{\mathit{H}}_{sk}$ and optional $r_1>0$, the Hessian ${\tilde{\mathit{H}}}_k$ is positive definite for all k. Thus, the cost function ${\mathit{\zeta }}_k^T{\tilde{\mathit{H}}}_k{\mathit{\zeta }}_k+{\tilde{\mathit{f}}}_k^T{\mathit{\zeta }}_k$ is strictly convex for all k. With that, the corresponding QP in Equation (23) is an SCQP for all k.

Due to the strict convexity property, the optimal solution for the QP in Equation (23) can be efficiently obtained using existing convex optimization tools and solvers [17], which means the QP is feasible even with the estimated states of T. □

It is worth stating that although the error in the estimation does not affect the feasibility of the QP problem, it affects the prediction of the target’s direction and determines whether the target can be intercepted quickly and accurately. Since the estimated states will converge in a fixed time, which can usually be limited to a few seconds, errors have little effect on the overall guidance process. Unlike general control processes, the guidance process does not need to deliberately follow a particular trajectory, but only needs to ensure that the relative distance to the target decreases gradually. In Equation (13), $V_r$ is limited to less than 0, which means that r is decreasing, satisfying the basic requirement that the UAV approaches the target gradually.

3.2.3. Maximum Input Determination

In contemplating an exceedingly extreme scenario, ${\alpha }_M$, the heading angle of M has deviated significantly from the LOS angle $\theta$ to the extent that even with the application of maximum control input, it is not feasible to nullify the deviation angle for the remainder of the trajectory. Despite this limitation, it remains imperative to employ the maximum control input in an attempt to effect a timely correction to the trajectory. Considering constraints on control increments, the determination of maximum input is designed as

$$\begin{array}{l}let\left({\alpha }_M-\theta \right)\in \left(-\pi ,\pi \right),t_{go}={\displaystyle \frac{k_{t}r}{V_M}}\\ if\left(sign\left({\alpha }_M-\theta \right)\left({\alpha }_M-\theta -{\displaystyle \frac{\left|u_{max}\right|t_{go}}{V_M}}sign\left({\alpha }_M-\theta \right)\right)\right)>0\\ u=u_{pre}-sign\left({\alpha }_M-\theta \right)\left|\mathsf{\Delta }{u}_{max}\right|\end{array}$$

(41)

where $u_{pre}$ represents the last control input, and $k_t>0$ is an influencing factor that determines how early the effect will take place.

4. Simulation Results

In this section, we present simulation results for the NMPC-IM guidance algorithm. Firstly, we discuss the selection of the desired output $y_d$ in NMPC-IM. The NMPC-IM guidance algorithm without SMO and Maximum Input Determination can be called pure NMPC-IM, which is compared with the PNG to verify the effectiveness of the algorithm. We then analyzed the parameter selection for Maximum Input Determination. The NMPC-IM algorithm is compared with PNG and NMPC and conclusions are drawn. Finally, assuming that the maneuvering states of the target are unknown, SMO is added to estimate them and obtain a satisfactory guidance effect.

The experiments were conducted on a Dell computer equipped with an Intel(R) Core(TM) i7-8565U CPU, operating at a base frequency of 1.80 GHz. The system is furnished with 8 GB of RAM and runs on a 64-bit operating system. For the simulation and computational tasks, Python 3.7 was utilized as the primary programming language. Moreover, we use the Euler method for numerical solutions.

4.1. Experimental Settings and Datasets

We consider the maneuvering target with time-varying acceleration. The basic simulation conditions are shown in

Table 1. Simulation conditions.

Variables	Value	Variables	Value
$\mathsf{\Delta }t$ 1	0.1 s	$r_b$2	10 m
$V_M$	14 m/s	$V_T$	10 m/s
$u_{max}$	2 g	$u_{min}$	−2 g
$\mathsf{\Delta }{u}_{max}$	1 g	$\mathsf{\Delta }{u}_{min}$	−1 g
$p\left(NMPC\right)$	50	$c\left(NMPC\right)$	25
$q\left(NMPC\right)$	0.1	$r\left(NMPC\right)$	1
$r_1$ 3	100	$r_2$	200
$p(NMPC-IM)$	Equation (25)	$c(NMPC-IM)$	[p/2]
$p_{max}$	50	$p_{min}$	3
$q,r,s(NMPC-IM)$	Equation (24)	$K\left(PNG\right)$	3
$a_{M0}$ 4	0 g	$r_0$ 5	50 m
${\dot{a}}_T$ 6	$\left(a_T\left(k\right)-a_T\left(k-1\right)\right)/\mathsf{\Delta }t$	${\dot{a}}_T\left(0\right)$	0

¹ Sample time. ² The simulation will stop once $r<r_b$. ³ Slack variable. ⁴ The initial control input, which can also be the acceleration of M. ⁵ A threshold used in Equation (25). ⁶ Euler discretization.

. Additionally, PNG is set as a comparison test, and the proportional coefficient is K.

4.2. Selection of the Desired Output $y_d$

Because there is a determination of attack range, r will not converge to 0. With the protection of determination conditions, and after experiments, it was found that it is more appropriate to use $y_d=1/\sqrt{r}$ rather than $y_d=\sqrt{r}$ for the algorithm. To verify this statement, we set the initial range and LOS angle to T as $r_0=200$ m and ${\theta }_0={45}^{\circ }$, respectively. T moves with acceleration as $a_T=-10sin\left(2t\right)$ g, and the initial heading angle is ${180}^{\circ }$. The initial heading angle of M is set to four sets of experiments ${\alpha }_{M0}=0^{\circ },{90}^{\circ },{180}^{\circ },{270}^{\circ }$. The results are shown in Figure 6.

According to the results, we can find that when choosing $y_d=1/\sqrt{r}$, the UAV intercepts the target faster than choosing $y_d=\sqrt{r}$. Especially in the cases of ${90}^{\circ }$ and ${270}^{\circ }$, the average energy expenditure is much smaller when choosing $y_d=1/\sqrt{r}$. Thus, it seems that the previous choice of $y_d$ is relatively correct.

4.3. Compare pure NMPC-IM with PNG

In order to better validate the performance of pure NMPC-IM, we set up two different scenarios with the target acceleration set to $a_T=20sin\left(t\right)$ g and $a_T=15sin\left(1.5t\right)$ g, respectively, based on the conditions $r_0=200$ m, ${\theta }_0={45}^{\circ }$, and the initial heading angle of T is $0^{\circ }$. The initial heading angles of M are set to four sets of experiments ${\alpha }_{M0}=0^{\circ },{90}^{\circ },{180}^{\circ },{270}^{\circ }$, and ${\alpha }_{M0}=0^{\circ },{30}^{\circ },{60}^{\circ },{90}^{\circ }$, respectively. The results are shown in Figure 7 and

Table 2. Impact time and average energy expenditure in different guidance laws ($a_T=20sin\left(t\right)$ g).

Guidance Laws	$0^{\circ }$	${90}^{\circ }$	${180}^{\circ }$	${270}^{\circ }$
$PNG$	$t_f=16.0$ s	$t_f=15.4$ s	$t_f=30.6$ s	$t_f=32.8$ s
	$su$ 1 = 3.1628 g	$su=2.7524$ g	$su=2.5487$ g	$su=3.3881$ g
$NMPC-IM$	$t_f=16.8$ s	$t_f=68.0$ s	$t_f=32.4$ s	$t_f=41.3$ s
	$su=1.6464$ g	$su=1.4575$ g	$su=1.5242$ g	$su=1.6239$ g

¹ Refer to Equation (18).

, and Figure 8 and

Table 3. Impact time and average energy expenditure in different guidance laws ($a_T=15sin\left(1.5t\right)$ g).

Guidance Laws	$0^{\circ }$	${30}^{\circ }$	${60}^{\circ }$	${90}^{\circ }$
$PNG$	$t_f=33.2$ s	$t_f=31.1$ s	$t_f=30.8$ s	$t_f=31.5$ s
	$su=3.2696$ g	$su=3.3183$ g	$su=3.2152$ g	$su=2.7818$ g
$NMPC-IM$	$t_f=36.2$ s	$t_f=33.3$ s	$t_f=30.0$ s	$t_f=33.5$ s
	$su=1.2263$ g	$su=0.9766$ g	$su=0.4707$ g	$su=1.2486$ g

, respectively.

From the results, it can be found that the impact time of PNG is shorter than that of NMPC-IM in most cases, but the average energy expenditure is larger. It is worth noting that in the ${90}^{\circ }$ counterpart case in Figure 7a, the guidance trajectory of the NMPC-IM re-intercepts the target after making a large circle, resulting in a much longer impact time than normal. This is due to the fact that the UAV’s miss distance on its first approach to the target is too large to intercept the target, resulting in a missed intercept. As for the reason, it could be that the control inputs are constrained, and it is difficult to increase the inputs without limit to make adjustments when faced with such a situation, as in the case of the PNG. The Maximum Input Determination proposed before is intended to solve this problem.

4.4. Effectiveness Validation of Maximum Input Determination

The effectiveness of Maximum Input Determination is validated first. We set the initial range and LOS angle to T as $r_0=200$ m and ${\theta }_0={120}^{\circ }$, respectively. The initial heading angle of M is ${\alpha }_{M0}={45}^{\circ }$. T moves with acceleration as $a_T=10sin\left(t\right)$ g, and the initial heading angle is $0^{\circ }$. After approximate testing, a range of more pronounced results due to changes in $k_t$ was identified. Within that range, a series of representative points were selected so that the results were distinguishable from each other. $k_t$ is chosen sequentially as $0.2,0.5,0.9,1.5,2.5$, and the results are shown in Figure 9.

We use colored spots to mark the portion of Maximum Input Determination, and it is obvious that when $k_t$ is relatively small, the maximum input is used at the very beginning of the guidance in order to head to T, but the trajectory is correspondingly more curved. As $k_t$ gradually increases, the Maximum Input Determination effects occur later. In general, at the end stage of guidance, since both M and T are dynamic, it is easier to produce angular deviation and the Maximum Input Determination effects occur more easily.

4.5. Comparison of NMPC-IM with PNG and NMPC

To compare NMPC-IM with PNG and NMPC, we designed two sets of experiments with different maneuvering states of T for a more complete validation of the algorithm.

We first set T moving with acceleration as $a_T=10sin\left(1.5t\right)$ g, and the initial heading angle as $0^{\circ }$. $r_0=150$ m and ${\theta }_0={45}^{\circ }$, respectively. The initial heading angle of M is ${\alpha }_{M0}={90}^{\circ }$. $k_t$ is chosen as $0.9$, for a more balanced effect as shown in Figure 9, i.e., not determined too early or too late. Additionally, we use the computation time to represent the real time it takes for the computer to finish running. The results are shown in Figure 10 and

Table 4. Impact time and average energy expenditure in different guidance laws (sin).

Guidance Laws	${\mathit{t}}_{\mathit{f}}$ (s)	$\mathit{su}$ (g)	Computation Time 1 (s)
$PNG$	30.3	2.6829	0.2784
$NMPC$	30.7	1.5450	3.3304
NMPC-IM	31.8	1.2876	1.5531

¹ After five repetitions of the experiment and taking the average.

.

Let $q=0$, which means completely disregarding the cost function in the original NMPC, to show the IM individually.

As shown in Figure 11, it can be seen that the IM is slow with respect to the target’s maneuverability, and will not respond very quickly. Comparing with Figure 10, we can find that IM works mainly by attenuating the implication effect of frequent oscillations of control inputs due to frequent maneuvers of the target (e.g., sometimes left, sometimes right), so that unnecessary energy loss can be reduced by certain prediction.

We then set T moving with acceleration as $a_T=20exp(-0.1t)$ g, and the initial heading angle as $0^{\circ }$. $r_0=180$ m and ${\theta }_0={45}^{\circ }$, respectively. The initial heading angle of M is ${\alpha }_{M0}={90}^{\circ }$. $k_t$ is chosen as $0.9$, for a more balanced effect as shown in Figure 9, i.e., not determined too early or too late. The results are shown in Figure 12 and

Table 5. Impact time and average energy expenditure in different guidance laws (exp).

Guidance Laws	${\mathit{t}}_{\mathit{f}}$ (s)	$\mathit{su}$ (g)	Computation Time (s)
$PNG$	15.0	3.0351	0.2098
$NMPC$	14.7	1.5293	1.1526
NMPC-IM	14.6	1.4304	0.9277

.

It can be seen from the above results that the addition of IM does improve the trajectory guided by NMPC. Although the trajectory of NMPC looks straighter than that of NMPC-IM in Figure 10a, the trajectory is still more curved from the details due to the fact that its control inputs have been oscillating between poles. More importantly, NMPC-IM reduces average energy expenditure during the UAV interception. In general, for $su$, NMPC-IM < NMPC < PNG, because unnecessary bends as well as control waste can be reduced with prediction. In terms of computational performance as shown in Table 3 and Table 4, NMPC and NMPC-IM are considerably slower than PNG, but NMPC-IM has some improvement compared to NMPC. We believe that it may be the inclusion of IM that makes the QP problem easier to solve, as well as the Maximum Input Determination to correct the trajectory in time and optimize the subsequent computation. All in all, the addition of IM is an improvement to the NMPC guidance method in terms of trajectory optimization and expenditure of average energy.

4.6. Further Testing of The NMPC-IM Algorithm

In order to check the universality and robustness of the algorithm, we consider the guidance effect in the case of different initial heading angles of M and with the addition of disturbance observer, respectively.

4.6.1. Different Initial Heading Angles of M

The acceleration of T is set as $a_T=20exp(-0.2t)$ g, with initial heading angle $0^{\circ }$. $r_0=250$ m and ${\theta }_0={45}^{\circ }$, respectively. The initial heading angle of M is ${\alpha }_{M0}={90}^{\circ }$. $k_t$ is chosen as $0.9$, and the results are shown in Figure 13.

With different initial heading angles, M always intercepts T, though $t_f$ differs.

4.6.2. Add Disturbance Observer

Previously, we assumed that the maneuvering states of T are known, including $a_T$ and ${\alpha }_T$. In this simulation, the maneuvering states of T are unknown; the fixed-time sliding mode observer in Equation (27) is used to estimate them. Take the same conditions as Figure 10 except $\mathsf{\Delta }t=0.02$ s, and set $k_1=5$, $k_2=5$, $k_3=15$, $p=1.5$, $q=0.5$; the results are shown in Figure 14.

Based on the comparison of the estimated and true values, it can be observed that the disturbance observer has a relatively good estimate of the maneuvering state of T, converging to the vicinity of the true value in a very short time, despite the inevitable presence of buffeting. By obtaining the approximate maneuvering information of T, under the guidance of the NMPC-IM algorithm, M is indeed able to accomplish the task of pursuing T, and the average energy expenditure remains smaller than NMPC. It is worth noting that, due to the buffeting, the estimation of the state of T may sometimes be diametrically opposed to the true value, but this does not immediately have an irreversible impact on the control inputs, further illustrating the robustness of NMPC-IM.

5. Conclusions

In this paper, a guidance method based on NMPC is developed, on the basis of which IM is proposed. IM is formulated as an additional cost function for incorporation into the NMPC, which is finally transformed into a complete QP problem. Meanwhile, the parameters are adjusted and optimized to make IM play a better role in NMPC. We add a determination condition, which can be regulated, for the timely correction of the guidance of NMPC-IM. Moreover, a fixed-time sliding mode disturbance observer is designed to reduce the information required and verify the robustness of this algorithm. It has been shown by simulation experiments that the inclusion of IM reduces the average energy expenditure by minimizing the impact of the target’s maneuverability on guidance. This holds positive significance for the UAV’s trajectory optimization during the interception of a more maneuverable target under constraints. The parameters used in this paper require further qualitative analysis and optimization through alternative methods to achieve superior results. Additionally, predicting the target’s maneuverability remains a significant challenge, and the current predictions are quite limited. Accurately predicting the target’s state in the future would represent a significant optimization.