Terminal Impact Time Control Cooperative Guidance Law for UAVs under Time-Varying Velocity

Abstract

Aiming at the problem that multiple Unmanned Aerial Vehicles (UAVs) attack the stationary target cooperatively under time-varying velocity, the cooperative guidance law with finite time convergence on two-dimensional plan and the three-dimensional cooperative guidance laws with impact time constraint are designed separately in this paper. Firstly, based on the relative motion equation between UAV and target on two-dimensional plane, the time cooperative guidance model of multiple UAVs is established. Then based on the consistency theory and graph theory, a distributed time cooperative guidance law is designed, which can ensure that the impact time of all UAVs can be quickly consistent in a limited time. Next, the cooperative guidance problem is expanded from two-dimensional plane to three-dimensional space, the motion model between UAV and target in three-dimensional space is established and the expression of time-to-go estimation under time-varying velocity is derived. Finally, according to whether there is the communication among UAVs under the condition of time-varying velocity, a multiple UAVs three-dimensional cooperative guidance law based on desired impact time and a multiple UAVs three-dimensional cooperative guidance law based on coordination variables are designed, respectively. The simulation results show that the cooperative guidance law with finite time convergence on two-dimensional plan and the three-dimensional cooperative guidance law with impact time constraint proposed in this paper are effective, which can both realize the saturation attack under the time-varying velocity.

1. Introduction

Multiple UAVs impact time control cooperative guidance (simultaneous attack) is a difficult problem in the current research. According to whether there is communication among formation members, the guidance methods of impact time control are divided into two categories in [1]: The first one is the impact time control method, which assigns the same designated time to all members before launching, and each member completes the attack task independently; the second one is the cooperative guidance method, in which the members exchange information through the communication topology network, adjust their movements according to the cooperative guidance algorithm, and finally reach the agreement of impact time.

The method of simultaneous attack by designing impact time control guidance law appeared earlier. In [2,3], an impact time control item was added to the proportional guidance law in order to attack the target with a designated time. In reference [4], the impact time control guidance law with proportional coefficient of 3 in reference [2] was extended to arbitrary value, and the accuracy of impact time control was improved by a high-order time-to-go estimation formula. In [5], a nonlinear sliding mode control method was adopted, and the sliding mode surface was weighted by the non-zero function of the relative distance and the ideal time-to-go to achieve the impact time requirements. By simplifying the relationship between the time of flight and velocity, an impact time control guidance law based on integral sliding mode control was proposed under varying velocity in [6], in which, only needing the range of the velocity, the impact time control can be achieved. However, the gravity and induced resistance of the member were not considered in the derivation process of the guidance law, which brought some limitations in practical application. In [7], an impact time control guidance law based on field-of-view shaping under varying velocity was proposed. The guidance law used an online iterative algorithm to predict the remaining average velocity, and then updated the adaptive gain to control the flight time. The time-to-go estimation algorithm was based on proportional guidance, which was completed by online iteration. This algorithm can improve the accuracy and velocity of time prediction.

In recent years, with the development of optimization theory, some optimization approaches such as Reinforcement Learning, Model predictive control (MPC), Model predictive static programming (MPSP) and Convex optimization have been applied to the guidance of UAVs. A training framework for the impact time control guidance law based on reinforcement learning theory was proposed in [8,9], which was robust to parameter uncertainty. In order to simplify the design of guidance law and improve the robustness of guidance law, reference [10] designed the terminal guidance law of interceptor based on meta-reinforcement learning. In [11], the memory pool generation method of traditional depth Q network (DQN) was improved, and a mid-course penetration control model of ballistic missile based on Markov decision process was proposed. In reference [12], Taylor expansion was performed on the model used for prediction on MPC to obtain an approximate linear model, and predictive control was applied to the linearized model. In references [13,14,15], considering the movement of the target in the prediction model, the tracking and landing of the moving target by unmanned aerial vehicle is realized by MPC. Reference [16] studied the suboptimal guidance method with terminal angle constraint based on MPSP, and gave the strategy of selecting the initial value of the control quantity based on the extended proportional guidance method; the ground stationary and maneuvering targets were also tested. In reference [17], with the minimum control energy as the performance index, a nonlinear suboptimal guidance law was derived based on MPSP method, which could satisfy both terminal position and terminal angle constraints. Based on convex optimization, a fast guidance law optimization algorithm with impact time and falling angle constraint was proposed in reference [18].

When the impact angle constraint is considered, the terminal impact angle control guidance laws are designed in references [19,20,21]. The generalized form of energy minimization optimal guidance law was derived in reference [19], and on this basis, an optimal impact angle control guidance law without delay was proposed. Reference [20] designed a controller with double-loop structure, based on which a three-dimensional impact angle control guidance law was proposed. Reference [21] expanded the circular navigation guidance law and realized the control of terminal impact angle.

The above guidance laws need to design the desired impact time in advance, and there is no information exchange during the flight of the members. In view of this, a time cooperative guidance architecture based on the “leader-followers” mode was designed in [22]. The coordination variable was the time-to-go of leader, and for followers, a method for calculating the rate of change of the line-of-sight angle was designed, so as to meet the requirement of the followers following the leader and attacking the target at the same time. In [23], the desired time-to-go was directly set as the average of the time-to-go of each member, thereby a hybrid guidance law satisfying both impact time and impact angle was designed.

However, the battlefield environment is relatively complex, and the centralized communication network adopted in [22,23] is difficult to effectively guarantee, so the distributed network topology is often used. When considering the target maneuvering in a two-dimensional plane, a distributed cooperative guidance law was designed in [24,25]. The guidance law consists of two parts: one is the local guidance law based on augmented proportional navigation; the other is the decentralized coordination strategy based on the principle of network synchronization. When considering the dynamic characteristics of the UAVs, the “leader-followers” mode was adopted in [26], the coordination strategy based on the consistency algorithm and the cooperative guidance law based on the guidance and control integration model were designed. In [27,28,29], a terminal cooperative guidance law under directed topology was proposed. The guidance law consists of two parts: the line of sight guidance law based on the multi-agent control theory and the line of sight normal guidance law based on secondary system stability and sliding mode theory.

According to existing literature, both the impact time control and cooperative guidance methods are obtained under the assumption that the velocity is constant or even adjustable [30,31]. However, in the actual process, the UAV is subject to air resistance and the velocity changes significantly. Therefore, the above research results are not applicable to the time coordination problem when the velocity of the unpowered vehicle changes with the flight state. At present, there are few open literatures on the impact time control problem of unpowered aircraft with varying velocity: In [32], based on the reference trajectory designed in the altitude-velocity profile, the range-to-go and time-to-go were predicted online and maneuvered laterally to correct the trajectory. A double-layer guidance structure was adopted in [33], based on the prediction-correction guidance law, the time-to-go was predicted online by neural network, and the time error was corrected by dynamically adjusting the line-of-sight angle corridor. Considering that the rate of change of velocity is a quadratic function of velocity, Reference [34] proposed a three-dimensional impact time cooperative guidance law based on desired impact time.

In order to solve the problem of multiple UAVs attack the stationary target cooperatively under time-varying velocity, the cooperative guidance law with finite time convergence on two-dimensional plan and the three-dimensional cooperative guidance laws with impact time constraint are designed separately in this paper. The main contributions of this paper are as follows:

(1)

The cooperative guidance method on two-dimensional plane proposed in this paper takes into account three aspects, namely, the finite time consistency of impact time, the time-varying velocity and the distributed communication topology. In addition, this paper adopts the communication topology without the leader, and each member in the cluster has the same status and function. Even if some UAVs in the cluster are destroyed, the impact time of the remaining clusters can still be quickly consistent.

(2)

Compared with reference [34], when the rate of change of velocity is a quadratic function of velocity, a three-dimensional impact time cooperative guidance law based on coordination variables is proposed in this paper, which does not need to set the desired impact time in advance as the impact time is determined by negotiation among UAVs.

(3)

When the rate of change of velocity is the first order function of velocity, the expression of time-to-go estimation is derived firstly, then according to whether there is the communication among UAVs, a multiple UAVs three-dimensional cooperative guidance law based on desired impact time and a multiple UAVs three-dimensional cooperative guidance law based on coordination variable are designed respectively, which can both realize the saturation attack.

The rest of this paper is arranged as follows: In the second part, multiple UAVs time cooperative guidance model is established on two-dimensional plane, and then a cooperative guidance law with finite time convergence under time-varying velocity is designed. In the third part, the cooperative guidance problem is expanded from two-dimensional plane to three-dimensional space, the motion model between UAV and target in three-dimensional space is established and the expression of time-to-go estimation under time-varying velocity is derived. The three-dimensional cooperative guidance law based on the desired impact time and the three-dimensional cooperative guidance law based on the coordination variables are designed respectively under time-varying velocity. Several numerical simulation examples are provided in the fourth part. The final part gives the conclusion of this study.

2. Cooperative Guidance Law with Finite Time Convergence under Time-Varying Velocity on Two-Dimensional Plan

2.1. Mathematical Model of benween UAV and Target on Two-Dimensional Plan

In order to establish the terminal motion Equation of the UAV under time-varying velocity, the following assumptions are made firstly.

Assumption 1.

The UAV and the target are regarded as mass points.

Assumption 2.

The target is considered static.

In Figure 1, $M$ and $T$ represent the missile and the target, respectively, and the mathematical model of motion between them can be expressed as

$$\dot{r}=-V_M\cos {\sigma }_M$$

(1)

$$\dot{q}=-\frac{V_M\sin {\sigma }_M}{r}$$

(2)

$${\dot{\gamma }}_M=a_M/V_M$$

(3)

$${\dot{V}}_M=-\left(D+mg\sin {\gamma }_M\right)/m$$

(4)

$${\sigma }_M={\gamma }_M-q$$

(5)

where, $V_M$ is the UAV velocity, $r$ is the range-to-go. ${\gamma }_M$ is the flight path angle, $q$ is the line-of-sight angle and ${\sigma }_M$ is the leading angle, respectively. $a_M$ is the acceleration, which is normal to the UAV velocity vector. $L$ and $D$ represent the aerodynamic lift and drag respectively. $m$ is the mass of UAV and $g$ is the gravitational acceleration.

The aerodynamic drag $D$ can be expressed as

$$D=C_D\rho V_M^2{S}_{ref}/2$$

(6)

where, $C_D$ is the drag coefficient, $\rho$ is atmospheric density and $S_{ref}$ is the aerodynamic reference area of the UAV.

It can be known from [35] that when $C_D$ is considered as constant $C_{D0}$ and $\rho$ is also considered as constant, the change rate of velocity can be expressed as

$${\dot{V}}_M=-K_M{V}_M^2$$

(7)

where, $K_M=C_{D0}\rho S_{ref}/2m$ is the coefficient of ${\dot{V}}_M$.

When the influence of aerodynamic drag and gravity on UAV velocity is comprehensively considered, it is pointed out in [36] that the range of velocity decrease is directly proportional to the current velocity, and the change rate of velocity can be expressed as

$${\dot{V}}_M=-K_{M1}{V}_M^{}$$

(8)

where, $K_{M1}$ is the deceleration control parameter, and the parameter size varies from 0.001 to 0.125.

The time-to-go estimation method in reference [34] is adopted in this paper, which can be expressed as

$$t_{go}=\frac{r}{V_M}\left(1+N^{\prime }{\sigma }_M^2\right){\left(\frac{{\sigma }_M}{\sin {\sigma }_M}\right)}^{\frac{1}{N-1}}$$

(9)

Defining

$$N^{\prime }=\frac{2-N}{6\left(N-1\right)\left(2N-1\right)}$$

(10)

where, $N$ is the proportional guidance coefficient, all the following $N$ have the same meaning.

2.2. Design of Cooperative Guidance Law with Finite Time Convergence

Terminal strike time $t_f{}_i$ of the $ith$ UAV can be expressed as

$$t_f{}_i=t_{goi}+t$$

(11)

where, $t_{goi}$ is the actual time-to-go of the $ith$ UAV.

It can be obtained from Equation (11)

$$t_f{}_i-t_f{}_j=t_{goi}-t_{goj}$$

(12)

As can be seen from Equation (12), when $t_{goi}$ of all UAVs is consistent, cooperative attack can be realized.

It can be obtained from Equation (1)

$$\left|{\sigma }_i\right|=\arccos \left(\frac{-{\dot{r}}_i}{V_i}\right)$$

(13)

Substitute Equation (13) into Equation (9)

$${\widehat{t}}_{go}=\frac{r_i}{V_i}\left(1+N^{\prime }{(-{\dot{r}}_i/V_i)}^2\right){\left(\frac{-{\dot{r}}_i/V_i}{\sin (-{\dot{r}}_i/V_i)}\right)}^{\frac{1}{N-1}}$$

(14)

Defining

$$x_i=\frac{r_i}{V_i},\hspace{0.17em}\hspace{0.17em}{v}_i=\frac{{\dot{r}}_i}{V_i}$$

(15)

It can be obtained from Equations (14) and (15)

$${\widehat{t}}_{go}=x_i\left(1+N^{\prime }{(-v_i)}^2\right){\left(\frac{-v_i}{\sin (-v_i)}\right)}^{\frac{1}{N-1}}$$

(16)

It can be seen from Equation (16) that the aim of cooperative guidance law with finite time convergence can be transformed into: Control $x_i$ and $v_i$ to achieve consistency in finite time, which is

$$\underset{t\to T}{\lim }||x_i-x_j||=0,\hspace{0.17em}\hspace{0.17em}\underset{t\to T}{\lim }||v_i-v_j||=0,\hspace{0.17em}\hspace{0.17em}i,j=1,2,\cdots n,\hspace{0.17em}T>0$$

(17)

It can be obtained from Equation (15)

$${\dot{x}}_i=v_i$$

(18)

It can be obtained from Equations (1), (5) and (15)

$$\begin{array}{ll}{\dot{v}}_i& ={\dot{\sigma }}_i\sin {\sigma }_i=\left({\dot{\gamma }}_i-{\dot{q}}_i\right)\sin {\sigma }_i\\ & =\frac{-a_i\sin {\sigma }_i}{V_i}+\frac{V_i{\sin }^2{\sigma }_i}{r_i}\end{array}$$

(19)

Substitute Equation (2) into Equation (19)

$${\dot{v}}_i=-\frac{a_i{r}_i{\dot{q}}_i{}_{}}{V_i{}^2}+\frac{r_i{\dot{q}}_i{}^2}{V_i}$$

(20)

Combining Equations (18) and (20), the second-order time cooperative guidance model for multiple UAVs is as follows:

$$\{\begin{array}{l}{\dot{x}}_i=v_i\\ {\dot{v}}_i=-\frac{a_i{r}_i{\dot{q}}_i}{V_i{}^2}+\frac{r_i{\dot{q}}_i{}^2}{V_i}\end{array}$$

(21)

The communication topology among multiple UAVs can be represented by undirected graph $G(A)=(\upsilon ,\xi ,A)$, in which $\upsilon$ describes the set of nodes, $\xi$ represents the connection between two nodes, matrix $A=[a_{ij}]\in R^{n\times n}$ represents the weight coefficient matrix. If information can be exchanged between the ith UAV and the jth UAV, $a_{ij}=1$; otherwise, $a_{ij}=0$. In particular, $a_{ii}=0,\hspace{0.17em}\hspace{0.17em}i=1,2,\cdots ,n$. Since $G(A)$ is an undirected graph, there is $a_{ij}=a_{ji}$. If there is at least one path between any two nodes in $G(A)$, the whole graph is connected. The Laplace matrix corresponding to undirected graph $G(A)$ among multiple UAVs is defined as $L=[l_{ij}]\in R^{n\times n}$, Where the elements of the matrix are

$$l_{ij}=\{\begin{array}{l}{\displaystyle \sum _{k=1,k\ne i}^n{a}_{ik}},j=i\hfill \\ -a_{ij},j\ne i\hfill \end{array}$$

(22)

The cooperative guidance law with finite time convergence is given in the form of Theorem 1.

Theorem 1. 

For the multiple UAVs system, if the communication topology$G(A)$is undirected and connected, the distributed time cooperative guidance law as shown in Equation (23) is designed to make the system states$x_i$and$v_i$converge to the same in a limited time.

$$a_i=\frac{V_i^2}{r_i{\dot{q}}_i}\left\{\begin{array}{l}-\alpha [{\displaystyle \sum _{j=1}^n{l}_{ij}{x}_j}+sig{({\displaystyle \sum _{j=1}^n{l}_{ij}{x}_j})}^k]-\beta [{\displaystyle \sum _{j=1}^n{l}_{ij}{v}_j}+sig{({\displaystyle \sum _{j=1}^n{l}_{ij}{v}_j})}^k]\\ -\gamma [sign(x_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{x}_j})+sign(v_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{v}_j})]\end{array}\right\}+V_i{\dot{q}}_i$$

(23)

where, $\alpha >0,\hspace{0.17em}\beta >0,\hspace{0.17em}\gamma >0$ are the coupling strengths,$sig{(x)}^k={\left|x\right|}^{k}sign(x)$, $sign(x)$ is a symbolic function, $\left|x\right|$ is the absolute value of $x$and$0<k<1$. ${\displaystyle \sum _{j=1}^n{l}_{ij}{x}_j}={\displaystyle \sum _{j=1}^n{a}_{ij}(x_i(t)-x_j(t))}$ and ${\displaystyle \sum _{j=1}^n{l}_{ij}{v}_j}={\displaystyle \sum _{j=1}^n{a}_{ij}(v_i(t)-v_j(t))}$represent relative displacement information and relative velocity information, respectively. Different variations of$V_i$are shown in Equations (7) and (8) respectively.

2.3. Proof of Theorem 1

First, the following two lemmas are introduced.

Lemma 1. 

In [37]: Let$\overline{G}$be an undirected connected graph, then the corresponding Laplace matrix$L$is a positive semidefinite symmetric matrix, that is$L\ge 0$. In particular, the symmetric matrix has exactly one zero eigenvalue, and the left and right eigenvectors corresponding to the zero eigenvalue are$1_n{}^T$and$1_n$respectively, that is,$L{1}_n=0_n$and$1_n{}^{T}L=0_n{}^T$hold at the same time. The$n$eigenvalues of matrix$L$are nonnegative real numbers, which are recorded as

$$0={\lambda }_1(L)\le {\lambda }_2(L)\le \cdots \le {\lambda }_n(L)$$

(24)

where, ${\lambda }_2(L)$ is the algebraic connectivity of $\overline{G}$, the following formula holds and satisfies ${\xi }^{T}L\xi \ge {\xi }^T{\lambda }_2\xi$, $\xi$ is a positive vector.

$${\lambda }_2(L)=\underset{\xi \ne 0,\hspace{0.17em}{1}^T\xi =0}{\min }\frac{{\xi }^{T}L\xi }{\left|\right|\xi |{|}^2}$$

(25)

Lemma 2.

In [38]: For system

$$\dot{x}=f(x),\hspace{0.17em}f(0)=0,\hspace{0.17em}x={\left(x_1,x_2,\dots ,x_n\right)}^T\in R^n$$

(26)

in which $f(x)={\left[f_1(x),\dots ,f_n(x)\right]}^T$. Assuming that there is a continuous function $V(x):R^n\to [0,\infty )$ and real numbers $k>0(\beta >0),\hspace{1em}\alpha \in (0,1)$that satisfy the following Equation

$$\dot{V}(x)+k{(V(x))}^{\alpha }\le 0$$

(27)

or

$$\dot{V}(x)+kV(x)+\beta {(V(x))}^{\alpha }\le 0,x\in R^n$$

(28)

The system is considered as finite time stable. In addition, let $T(x_0)$ represent the adjustment time function, then the finite convergence time$T(x_0)$satisfies

$$T\left(x_0\right)\le \frac{V^{1-\alpha }\left(x_0\right)}{k(1-\alpha )}$$

(29)

or

$$T\left(x_0\right)\le \frac{1}{k(1-\alpha )}\ln \frac{k{V}^{1-\alpha }\left(x_0\right)+\beta }{\beta },x\in R^n$$

(30)

and for all$t>T(x_0)$, there is $V(x)=0$.

Proof. 

First, substitute Equation (23) into Equation (21) to obtain

$$\{\begin{array}{l}{\dot{x}}_i=v_i\\ {\dot{v}}_i=u_i=-\alpha [{\displaystyle \sum _{j=1}^n{l}_{ij}{x}_j}+sig{({\displaystyle \sum _{j=1}^n{l}_{ij}{x}_j})}^k]-\beta [{\displaystyle \sum _{j=1}^n{l}_{ij}{v}_j}+sig{({\displaystyle \sum _{j=1}^n{l}_{ij}{v}_j})}^k]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\gamma [sign(x_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{x}_j})+sign(v_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{v}_j})],\hspace{0.17em}\hspace{0.17em}i=1,2,\cdots ,n\end{array}$$

(31)

Let ${\overline{x}}_i=x_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{x}_j},\hspace{0.17em}{\overline{v}}_i=v_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{v}_j},\hspace{0.17em}\hspace{0.17em}i=1,2,\cdots ,n$ represent the relative position error and relative speed error of each UAV respectively. Then the error system Equation is

$$\{\begin{array}{l}{\dot{\overline{x}}}_i={\overline{v}}_i\\ {\dot{\overline{v}}}_i=u_i-\frac{1}{n}{\displaystyle \sum _{j=1}^n{u}_j},\hspace{0.17em}\hspace{0.17em}i=1,\cdots ,n\end{array}$$

(32)

Let $\overline{x}=[{\overline{x}}_1{}^T,{\overline{x}}_2{}^T,\cdots ,{\overline{x}}_n{}^T]=Mx$, $\overline{v}=[{\overline{v}}_1{}^T,{\overline{v}}_2{}^T,\cdots ,{\overline{v}}_n{}^T]=Mv$, in which $M=I-\frac{1}{n}\overline{N}$, $\overline{N}$ is an all 1 square matrix of order $n$.

Combine Equations (31) and (32) to obtain

$$\{\begin{array}{l}\dot{\overline{x}}=\overline{v}\\ \dot{\overline{v}}=-M\dot{v}=-Mu\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=-\alpha M[Lx+sig{(Lx)}^k]-\beta M[Lv+sig{(Lv)}^k]-\gamma M[sign(Mx)+sign(Mv)]\end{array}$$

(33)

According to $M=I-\frac{1}{n}\overline{N}$ and Lemma 1

$$LM=L=ML$$

(34)

Substitute Equation (34) into Equation (33)

$$\{\begin{array}{l}\dot{\overline{x}}=\overline{v}\\ \dot{\overline{v}}=-\alpha [L\overline{x}+sig{(L\overline{x})}^k]-\beta [L\overline{v}+sig{(L\overline{v})}^k]-\gamma M[sign(\overline{x})+sign(\overline{v})]\end{array}$$

(35)

Let $\overline{\zeta }={({\overline{x}}^T,{\overline{v}}^T)}^T$, then Equation (35) can be re-expressed as

$${\overline{\zeta }}^{\prime }=-B_1\overline{\zeta }-B_{2}sig{(\overline{L}\overline{\zeta })}^k-\gamma B_{3}sign(\overline{\zeta })$$

(36)

where

$$B_1=\left[\begin{array}{cc}{0}_n& -I_n\\ \alpha L& \beta L\end{array}\right]$$

(37)

$$B_2=\left[\begin{array}{cc}{0}_n& 0_n\\ \alpha I_n& \beta I_n\end{array}\right]$$

(38)

$$B_3=\left[\begin{array}{cc}{0}_n& 0_n\\ M& M\end{array}\right]$$

(39)

$$\overline{L}=\left[\begin{array}{cc}{L}_n& 0_n\\ 0_n& L_n\end{array}\right]$$

(40)

The following mainly proves Equation (36). The proof of Equation (36) is transformed into the proof of Theorem 2. □

Theorem 2.

When the control gain satisfies$\alpha /{\beta }^2<{\lambda }_2(L),\hspace{0.17em}0<k<1$, the undirected connected topology Equation (31) can achieve finite time consistency, that is, there are

$$x_i\to \frac{1}{n}{\displaystyle \sum _{j=1}^n{x}_j},\hspace{0.17em}{v}_i\to \frac{1}{n}{\displaystyle \sum _{j=1}^n{v}_j}$$

(41)

Proof. 

Since matrix M has a simple eigenvalue of 0, the corresponding right eigenvector is $1_n$, and 1 is the other $n-1$ multiple eigenvalues. Then if $\overline{x}=0$ and $\overline{v}=0$ hold, if and only if $x_1=x_2=\cdots =x_n$ and $v_1=v_2=\cdots =v_n$. Therefore, when the error variables $\overline{x}$ and $\overline{v}$ converge to zero in finite time, the finite time consistent control of the system can be obtained.

The following Lyapunov function is constructed

$$V_1(t)=\frac{1}{2}{\overline{\zeta }}^T{P}_1\overline{\zeta }$$

(42)

where, $P_1=\left[\begin{array}{cc}\frac{\alpha \beta }{n}(L+L^T)& \frac{\alpha }{n}{I}_n\\ \frac{\alpha }{n}{I}_n& \frac{\beta }{n}{I}_n\end{array}\right]$.

Expand Equation (42)

$$\begin{array}{l}{V}_1(t)=\frac{1}{2}{\overline{\zeta }}^T{P}_1\overline{\zeta }\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{\alpha \beta }{2n}{\overline{x}}^T(L+L^T)\overline{x}+\frac{\alpha }{n}{\overline{x}}^T{I}_n\overline{v}+\frac{\beta }{2n}{\overline{v}}^T{I}_n\overline{v}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\ge \frac{1}{2}{\overline{\zeta }}^T{Q}_1\overline{\zeta }\end{array}$$

(43)

where, $Q_1=\left[\begin{array}{cc}\frac{2\alpha \beta }{n}{\lambda }_2(L)I_n& \frac{\alpha }{n}{I}_n\\ \frac{\alpha }{n}{I}_n& \frac{\beta }{n}{I}_n\end{array}\right]$.

If the control parameters satisfy $\beta >0$ and ${\lambda }_2(L)>\alpha /2{\beta }^2$, according to Lemma 1, $Q_1>0$ can be obtained. It can be concluded from Equation (43) that Lyapunov function $V_1(t)\ge 0$, and if $V_1(t)=0$, if and only if $\overline{\zeta }=0_{2n}$.

Find the derivative of Lyapunov function $V_1(t)$ and substitute the control protocol (36) to obtain

$$\begin{array}{l}{V}_1^{\prime }=-\frac{1}{2}{\overline{\zeta }}^T(P_1{B}_1+B_1{}^T{P}_1)\overline{\zeta }-\frac{1}{2}{\overline{\zeta }}^T(P_1{B}_2+B_2{}^T{P}_1)sig{(\overline{L}\overline{\zeta })}^k\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{1}{2}{\overline{\zeta }}^T[\gamma (P_1{B}_2+B_2{}^T{P}_1)]sign(\overline{\zeta })\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=-\frac{1}{2}{\overline{\zeta }}^T\left[\begin{array}{cc}\frac{{\alpha }^2}{n}(L+L^T)& 0_n\\ 0_n& \frac{{\beta }^2}{n}(L+L^T)-\frac{2\alpha }{n}{I}_n\end{array}\right]\overline{\zeta }\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\overline{\zeta }}^T\left[\begin{array}{cc}\frac{{\alpha }^2}{n}{I}_n& \frac{\alpha \beta }{n}{I}_n\\ \frac{\alpha \beta }{n}{I}_n& \frac{{\beta }^2}{n}{I}_n\end{array}\right]sig{(\overline{L}\overline{\zeta })}^k\hspace{0.17em}-\frac{1}{2}{\overline{\zeta }}^T\left[\begin{array}{cc}\frac{2\alpha }{n}\gamma M& \frac{\alpha +\beta }{n}\gamma M\\ \frac{\alpha +\beta }{n}\gamma M& \frac{2\beta }{n}\gamma M\end{array}\right]sign(\overline{\zeta })\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\le -{\overline{\zeta }}^T{P}_2\overline{\zeta }-{\overline{\zeta }}^T{P}_3{\overline{L}}^{k}sig{(\overline{\zeta })}^k-\frac{1}{2}\gamma {\overline{\zeta }}^T{P}_{4}Msign(\overline{\zeta })\end{array}$$

(44)

According to the known condition $\left|\right|{\delta }_i\left|\right|\le \gamma <\infty$

$$V_1^{\prime }\le -{\overline{\zeta }}^T{P}_2\overline{\zeta }-{\overline{\zeta }}^T{P}_3{\overline{L}}^{k}sig{(\overline{\zeta })}^k$$

(45)

where

$$P_2=\left[\begin{array}{cc}\frac{{\alpha }^2{\lambda }_2(L)}{n}& 0\\ 0& \frac{{\beta }^2{\lambda }_2(L)-\alpha }{n}\end{array}\right]$$

(46)

$$P_3=\left[\begin{array}{cc}\frac{{\alpha }^2}{n}& \frac{\alpha \beta }{n}\\ \frac{\alpha \beta }{n}& \frac{{\beta }^2}{n}\end{array}\right]$$

(47)

$$P_4=\left[\begin{array}{cc}\frac{2\alpha }{n}& \frac{\alpha +\beta }{n}\\ \frac{\alpha +\beta }{n}& \frac{2\beta }{n}\end{array}\right]$$

(48)

From Lemma 1 to obtain

$$\begin{array}{l}{V}_1(t)=\frac{1}{2}{\overline{\zeta }}^T{P}_1\overline{\zeta }\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{\alpha \beta }{2n}{\overline{x}}^T(L+L^T)\overline{x}+\frac{\alpha }{n}{\overline{x}}^T{I}_n\overline{v}+\frac{\beta }{2n}{\overline{v}}^T{I}_n\overline{v}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\le \frac{\alpha \beta }{n}{\lambda }_{\max }(L){\overline{x}}^T{I}_n\overline{x}+\frac{\alpha }{n}{\overline{x}}^T{I}_n\overline{v}+\frac{\beta }{2n}{\overline{v}}^T{I}_n\overline{v}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\le {\overline{\zeta }}^T{P}_5\overline{\zeta }\end{array}$$

(49)

where

$$P_5=\left[\begin{array}{cc}\alpha \beta {\lambda }_{\max }(L)& \frac{\alpha }{2n}\\ \frac{\alpha }{2n}& \frac{\beta }{2n}\end{array}\right]$$

(50)

From Equation (50), the following Equation holds

$$V_1(t)\le {\lambda }_{\max }(P_5){\overline{\zeta }}^T\overline{\zeta }$$

(51)

where, take the parameter $m_3={\lambda }_{\max }(P_5)$.

From Equation (45) to obtain

$$V_1^{\prime }\le -{\lambda }_{\min }(P_2){\overline{\zeta }}^T\overline{\zeta }-{\lambda }_{\min }(P_3){\lambda }_2{}^k(L){\overline{\zeta }}^{T}sig{(\overline{\zeta })}^k$$

(52)

where, take the parameter $m_1={\lambda }_{\min }(P_2)$, $m_2={\lambda }_{\min }(P_3){\lambda }_2{}^k(L)$.

Combining Equations (51) and (52), we can get

$$V_1^{\prime }(t)\le -\frac{m_1}{m_3}{V}_1(t)-\frac{m_2}{m_3}{V}_1{(t)}^{\frac{1+k}{2}}$$

(53)

To summarize, Equation (53) satisfies the conditions of Lemma 2, that is, it can be obtained in finite time

$$\overline{\zeta }={({\overline{x}}^T,{\overline{v}}^T)}^T={(0,0)}^T$$

(54)

That is, in finite time

$$x_i\to \frac{1}{n}{\displaystyle \sum _{j=1}^n{x}_j},\hspace{0.17em}{v}_i\to \frac{1}{n}{\displaystyle \sum _{j=1}^n{v}_j},\hspace{0.17em}\hspace{0.17em}i=1,2,\cdots ,n$$

(55)

The specific convergence time is

$$T=\frac{2m_3}{m_1(1-k)}\ln \frac{m_1{V}^{1-k}(0)+m_2}{m_2}$$

(56)

To summarize, the cooperative guidance law expressed by Equation (23) with finite time convergence has been proved. □

3. Three-Dimensional Cooperative Guidance Law with Impact Time Constraint under Time-Varying Velocity

3.1. Mathematical Model of UAV in Three Dimensional Space

In this section, the cooperative guidance problem is expanded from two-dimensional plane to three-dimensional space, and the relative motion model between UAV and target is shown in Figure 2. The symbols in this section that are the same as those in the Section 2.1 have the same meanings.

Then the kinematics Equation of the UAV can be expressed as

$$\begin{array}{l}\dot{x}=V_M\cos {\gamma }_M\cos {\psi }_M\\ \dot{y}=V_M\sin {\gamma }_M\\ \dot{z}=V_M\cos {\gamma }_M\sin {\psi }_M\end{array}$$

(57)

The dynamic Equation of the UAV can be expressed as

$$\begin{array}{l}{\dot{V}}_M=-\frac{D}{m}-g\sin {\gamma }_M\\ {\dot{\gamma }}_M=\frac{L\cos \nu }{m{V}_M}-\frac{g\cos {\gamma }_M}{V_M}\\ {\dot{\psi }}_M=\frac{L\sin \nu }{m{V}_M\cos {\gamma }_M}\end{array}$$

(58)

where, $\nu$ is the bank angle.

The change rate of $r$, ${\psi }_L$ and $q$ can be expressed as

$$\begin{array}{l}\dot{r}=V_M\sin {\gamma }_M\sin q-V_M\cos {\gamma }_M\cos \left({\psi }_M-{\psi }_L\right)\cos q\\ {\dot{\psi }}_L=\frac{V_M\cos q\sin \left({\psi }_M-{\psi }_L\right)}{r\cos q}\\ \dot{q}=\frac{V_M\cos {\gamma }_M\sin \left({\psi }_M-{\psi }_L\right)\sin q-V_M\sin {\gamma }_M\cos q}{r}\end{array}$$

(59)

When the rate of change of velocity is the first order function of velocity, $a_z$ and $a_y$ are the acceleration of the pitch channel and the acceleration of the yaw channel respectively, then the dynamic equation of UAV can be simplified as

$$\begin{array}{l}{\dot{V}}_M=-K_{M1}{V}_M^{}\\ {\dot{\gamma }}_M=a_z/V_M\\ {\dot{\psi }}_M=a_y/V_M\end{array}$$

(60)

Similarly, when the rate of change of velocity is the quadratic function of velocity, the dynamic Equation is adopted in reference [34], which is

$$\begin{array}{l}{\dot{V}}_M=-K_M{V}_M^2\\ {\dot{\gamma }}_M=a_z/V_M\\ {\dot{\psi }}_M=a_y/V_M\end{array}$$

(61)

3.2. Time-to-Go Estimation

When the rate of change of velocity is the first order function of velocity, integrating Equation (8) to get

$$V_M\left(t\right)=V_{M0}{e}^{-K_{M1}t}$$

(62)

where, $V_{M0}$ is the initial velocity of the UAV.

From reference [34] we can get

$${\dot{\sigma }}_M=-\frac{\left(N-1\right)V_M\left(t\right)\sin {\sigma }_M}{r_0{\left(\frac{\sin {\sigma }_M}{\sin {\sigma }_{M0}}\right)}^{\frac{1}{N-1}}}$$

(63)

where, $r_0$ is initial range-to-go.

Substituting Equation (62) into Equation (63), yields

$${\dot{\sigma }}_M=K{e}^{-K_{M1}t}{\left(\sin {\sigma }_M\right)}^{\frac{N-2}{N-1}}$$

(64)

where, $K=-\frac{\left(N-1\right)V_{M0}}{r_0}{\left(\sin {\sigma }_{M0}\right)}^{\frac{1}{N-1}}$.

From Equation (64), we can get

$$e^{-K_{M1}t}dt=\frac{1}{K}{\left(\sin {\sigma }_M\right)}^{\frac{2-N}{N-1}}d{\sigma }_M$$

(65)

By integrating Equation (65), we can get

$$-\frac{1}{K_{M1}}{e}^{-K_{M1}t}=\frac{1}{K}{\displaystyle {\int }_{{\sigma }_{M0}}^{{\sigma }_M}{\left(\sin {\sigma }_M\right)}^{\frac{2-N}{N-1}}d{\sigma }_M}$$

(66)

From reference [34], we can get

$$\frac{1}{K}{\displaystyle {\int }_{{\sigma }_{M0}}^{{\sigma }_M}{\left(\sin {\sigma }_M\right)}^{\frac{2-N}{N-1}}d{\sigma }_M}\approx \frac{1}{K}{\displaystyle {\int }_{{\sigma }_{M0}}^{{\sigma }_M}\left({\sigma }_M^{\frac{2-N}{N-1}}+\frac{2-N}{N-1}\frac{{\sigma }_M^{\frac{N}{N-1}}}{6}\right)d{\sigma }_M}$$

(67)

Equation (66) can be further simplified to obtain

$$t=\frac{\ln \left(-K_{M1}\left(\frac{1}{K}{\displaystyle {\int }_{{\sigma }_{M0}}^{{\sigma }_M}\left({\sigma }_M^{\frac{2-N}{N-1}}+\frac{N-2}{N-1}\frac{{\sigma }_M^{\frac{N}{N-1}}}{6}\right)d{\sigma }_M}\right)+e^{-K_{M1}{t}_0}\right)}{K_{M1}}$$

(68)

where, $t_0$ is the initial time and ${\sigma }_{M0}$ is the initial leading angle.

When $r$ is zero, ${\sigma }_M$ is zero, so Equation (68) can be re-expressed as

$$t_{go}=\frac{\ln \left(-K_{M1}\left(\frac{r}{V_M}\left(1+\frac{2-N}{6\left(N-1\right)\left(2N-1\right)}{\sigma }_{{}_M}^2\right){\left(\frac{{\sigma }_M}{\sin {\sigma }_M}\right)}^{\frac{1}{N-1}}\right)+1\right)}{K_{M1}}$$

(69)

Then Equation (69) can be re-expressed as

$$t_{go}=\frac{\ln \left(-K_{M1}\left(\frac{r}{V_M}\left(1+N^{\prime }{\sigma }_{{}_M}^2\right){\left(\frac{{\sigma }_M}{\sin {\sigma }_M}\right)}^{\frac{1}{N-1}}\right)+1\right)}{K_{M1}}$$

(70)

where, $N^{\prime }$ has the same meaning as in Equation (10).

It is also pointed out in [34] that when the rate of change of velocity is the quadratic function of velocity, the time-to-go can be expressed as

$$t_{go}=\frac{e^{K_M{V}_M\left(\frac{r}{V_M}\left(1+N^{\prime }{\sigma }_{{}_M}^2\right){\left(\frac{{\sigma }_M}{\sin {\sigma }_M}\right)}^{\frac{1}{N-1}}\right)}-1}{K_M{V}_M}$$

(71)

3.3. Three-Dimensional Cooperative Guidance Law Based on Desired Impact Time

In this section, the main objective is to design the three-dimensional cooperative guidance law based on desired impact time when the rate of change of velocity is the first-order function of velocity.

The guidance law adopts the form in reference [34], which is as follows:

$$\begin{array}{l}{\dot{\gamma }}_M=N\dot{q}\\ {\dot{\psi }}_M=N{\dot{\psi }}_L\left(\frac{3}{2}-\frac{1}{2}\sqrt{1+\frac{240V_{Mxz}^5}{{\left(N{V}_{Mxz}{\dot{\psi }}_L\right)}^2{r}_{xz}^3}{e}_t}\right)\end{array}$$

(72)

where, $V_{Mxz}$ and $r_{xz}$ are the projection of velocity and the projection of the distance in the XZ plane, respectively. Where, $V_{Mxz}=V_M\cos {\gamma }_M$, $r_{xz}=r\cos q$. $e_t$ is the impact time error, which is

$$e_t=t_d-t-t_{go}$$

(73)

where, $t_d,\hspace{0.17em}t,\hspace{0.17em}{t}_{go}$ are the desired impact time, the current time and the estimation of the time-to-go in the yaw channel, respectively.

According to Equation (70), the $t_{go}$ in yaw channel can be expressed as

$$t_{go}=\frac{\ln \left(-K_{Mxz}\left(\frac{r_{xz}}{V_{Mxz}}\left(1+N^{\prime }{\eta }_{{}_M}^2\right){\left(\frac{{\eta }_M}{\sin {\eta }_M}\right)}^{\frac{1}{N-1}}\right)+1\right)}{K_{Mxz}}$$

(74)

where, ${\eta }_M$ is the leading angle in the yaw channel, which meets ${\eta }_M={\psi }_L-{\psi }_M$. $K_{Mxz}$ is the coefficient of rate of velocity change in the XZ plane, and which satisfies

$${\dot{V}}_{Mxz}=-K_{Mxz}{V}_{Mxz}^{}$$

(75)

Derive both sides of $V_{Mxz}=V_M\cos {\gamma }_M$ with respect to time

$${\dot{V}}_{Mxz}={\dot{V}}_M\cos {\gamma }_M-V_M{\dot{\gamma }}_M\sin {\gamma }_M$$

(76)

By Equations (75) and (76), $K_{Mxz}$ can be expressed as

$$K_{Mxz}=-\left({\dot{V}}_M\cos {\gamma }_M-V_M{\dot{\gamma }}_M\sin {\gamma }_M\right)/V_{Mxz}^{}$$

(77)

Up to now, the parameters in Equation (72) have clear meanings.

3.4. Three-Dimensional Cooperative Guidance Law Based on Coordinated Variables

In Section 3.3, the three-dimensional cooperative guidance law based on desired impact time needs to set the desired impact time for each UAV, there is no information exchange among members. For this purpose, the three-dimensional cooperative guidance law based on coordinated variables is designed in this section.

The desired impact time $\xi$ is selected as the coordination variable. Since the impact time constraint is controlled only through the yaw channel, the control energy consumption of the pitch channel is not considered, and the yaw channel acceleration can be expressed as

$$a_y=N{V}_{Mxz}{\dot{\psi }}_L\left(\frac{3}{2}-\frac{1}{2}\sqrt{1+\frac{240V_{Mxz}^5}{{\left(N{V}_{Mxz}{\dot{\psi }}_L\right)}^2{r}_{xz}^3}\left(\xi -t-t_{go}\right)}\right)$$

(78)

Using Taylor expansion series, Equation (78) can be re-expressed as

$$\begin{array}{ll}{a}_y& =N{V}_{Mxz}{\dot{\psi }}_L-\frac{60V_{Mxz}^5}{N{V}_{Mxz}{\dot{\psi }}_L{r}_{xz}^3}\left(\xi -t-t_{go}\right)\\ & =a_{y1}-a_{y2}\left(\xi -t-t_{go}\right)\end{array}$$

(79)

where, $a_{y1}=N{V}_{Mxz}{\dot{\psi }}_L,a_{y2}=\frac{60V_{Mxz}^5}{N{V}_{Mxz}{\dot{\psi }}_L{r}_{xz}^3}$.

The cost function of UAV $i\left(i=1,2,\dots ,n\right)$ is taken as

$$J_i\left(\xi \right)=a_{y,i}^2$$

(80)

The total cost function of formation can be expressed as the sum of control energy of each UAV when the control energy consumption of the pitching channel is not considered.

$$J_t\left(\xi \right)={\displaystyle \sum _{i=1}^n{a}_{y,i}^2}$$

(81)

Take the impact time that minimizes the total energy consumption of the formation as the desired impact time, namely

$${\xi }^{\ast }=\arg \min J_t\left(\xi \right)$$

(82)

It can be obtained from Equations (79) and (81)

$$J_t\left(\xi \right)={\displaystyle \sum _{i=1}^n{\left(a_{y1,i}-a_{y2,i}\left(\xi -t-t_{go,i}\right)\right)}^2}$$

(83)

Find the partial derivative of Equation (83) with respect to $\xi$ can be obtained

$$\begin{array}{ll}\frac{\partial J_t\left(\xi \right)}{\partial \xi }& =2{\displaystyle \sum _{i=1}^n{a}_{y,i}^{}}\frac{\partial a_{y,i}^{}}{\partial \xi }\\ & =2{\displaystyle \sum _{i=1}^n\left(-a_{y1,i}^{}{a}_{y2,i}^{}+a_{y2,i}^2\left(\xi -t-t_{go,i}\right)\right)}\end{array}$$

(84)

According to Equation (82), it can be obtained

$${\xi }^{*}=\frac{{\displaystyle \sum _{i=1}^n\left(a_{y2,i}^2\left(t+t_{go,i}\right)+a_{y1,i}{a}_{y2,i}\right)}}{{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2}}$$

(85)

Let $\delta ={\displaystyle \sum _{i=1}^n{a}_{y1,i}{a}_{y2,i}}/{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2}$, then Equation (85) can be re-expressed as

$${\xi }^{*}=\frac{{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2\left(t+t_{go,i}\right)}}{{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2}}+\delta$$

(86)

It can be seen from Equation (86) that ${\xi }^{\ast }$ consists of two parts. When the UAV is close to the target, compared with the first part, $\delta$ is very small, not in an order of magnitude. So $\delta$ can be ignored and the desired suboptimal solution of impact time can be obtained.

$${\xi }^{†}=\frac{{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2\left(t+t_{go,i}\right)}}{{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2}}$$

(87)

Let the weight of each UAV be $w_i$, and $w_i=a_{y2,i}^2/{\displaystyle \sum _{i=1}^n{a}_{y2,i}^2}$, then Equation (87) can be re-expressed as

$${\xi }^{†}={\displaystyle \sum _{i=1}^n{w}_i\left(t+t_{go,i}\right)}$$

(88)

To summarize, taking the desired impact time as the coordination variable, the cooperative guidance law designed based on the coordination variable can be expressed as

$$\begin{array}{c}{\dot{\gamma }}_M=N\dot{q}\hfill \\ {\dot{\psi }}_M=N{\dot{\psi }}_L-\frac{60V_{Mxz}^5}{N{V}_{{}_{Mxz}}^2{\dot{\psi }}_L{r}_{xz}^3}\left(\xi -t-t_{go}\right)\hfill \end{array}$$

(89)

4. Numerical SIMULATION

4.1. Performance Verification of Cooperative Guidance Law with Finite Time Convergence under Time-Varying Velocity on Two-Dimensional Plan

The effectiveness of the cooperative guidance law with finite time convergence under time-varying velocity on two-dimensional plan proposed in this paper is verified by numerical simulations. The simulation situation is as follows: in the horizontal plane, aiming at the stationary target, four UAVs fly cooperatively and finally hit the target. The available overload limit for UAVs is 10 g, and the parameters of cooperative guidance law are $\alpha =0.025,\hspace{0.17em}\beta =0.1,\hspace{0.17em}\gamma =0.5,\hspace{0.17em}k=0.2$.

The communication topology of the four UAVs is undirected and connected, which is as shown in Figure 3, and the corresponding weight coefficient matrix is shown in the following formula.

$$\mathrm{A}=\left[\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 1& 0\\ 0& 1& 0& 1\\ 0& 0& 1& 0\end{array}\right]$$

(90)

The position of the target is $(12,0)\hspace{0.17em}\mathrm{km}$, the initial conditions of the four UAVs are shown in

Table 1. Initial parameters of the four UAVs.

UAVs	Initial Position (m)	Velocity (m/s)	Initial Flight Path Angle (deg)
M1	(0, 1900)	242	−8
M2	(0, 400)	238	−3
M3	(0, −700)	240	3
M4	(0, −1100)	241	5

.

4.1.1. When the Rate of Change of Velocity Is the First Order Function of Velocity

The miss distance and strike time of four UAVs are given in

Table 2. Initial parameters of the four UAVs.

UAVs	Miss Distance (m)	Strike Time (s)
M1	0.45	61.03
M2	1.42	61.03
M3	0.78	61.03
M4	1.97	61.03

. It can be seen from the table that the maximum miss distance is 1.97 m and the maximum strike time deviation is less than 0.01 s, which reflects the effectiveness of the guidance law when the rate of change of velocity is the first order function of velocity.

The simulation results are shown in Figure 4. It can be seen from Figure 4a that all four UAVs can gradually fly to the target. As can be seen from Figure 4b, the distance between each UAV and the target is getting closer and closer, and finally decreases to 0, which indicates that the four UAVs can hit the target accurately. As can be seen from Figure 4c, the $\dot{r}/V_m$ of the four UAVs are inconsistent at the initial time, but reach the same after 8 s, which indicates that the cooperative guidance law can make the $\dot{r}/V_m$ reach the consistency quickly, and also means that the time-to-go of four UAVs reach the consistency within 8 s. It can be seen from Figure 4d that the time-to-go of the four UAVs can quickly converge to the consistency within 8 s, which shows that the cooperative guidance law makes the time-to-go reach the consistency quickly, and also shows the correctness of the analysis in Figure 4c. It can be seen from Figure 4e that the velocity of the four UAVs decrease linearly with time, which shows the effectiveness of the cooperative guidance law when the rate of change of velocity is the first order function of velocity. It can be seen from Figure 4f that the overload amplitude of the four UAVs is relatively large in the first 3.5 s. According to the analysis of Figure 4c, this is because the deviation of $\dot{r}/V_m$ is relatively large in the initial 3.5 s, and a large overload is required to make the $\dot{r}/V_m$ converge to the consistency quickly, which reflects the advantages of fast convergence and strong robustness of cooperative guidance law with finite time. In addition, although the overload command is large in initial 3.5 s, it is still within the available overload range of [−10 g, 10 g], so the algorithm has certain engineering application potential. There is a certain sudden change in the overload command around 10 s, which is caused by switching the cooperative guidance law to the proportional guidance law. After 10 s, under the action of proportional guidance law, the overload command amplitude of each UAV gradually decreases and converges to 0 in the final stage.

4.1.2. When the Rate of Change of Velocity Is the Quadratic Function of Velocity

When the rate of change of velocity is the quadratic function of velocity, as can be seen from Figure 5a,b, four UAVs gradually fly to the target and hit the target at the same time in 53.27 s, which proves the effectiveness of the cooperative guidance law. It can be seen from Figure 5e that the velocity of each UAV decreases more slowly than that in Figure 4e. Therefore, when the rate of change of velocity is the quadratic function of velocity, the UAVs can attack the target in a shorter time.

To summarize, the cooperative guidance law with finite time convergence under time-varying velocity on two-dimensional plan proposed in this paper can ensure that multiple UAVs reach the target at the same time, the miss distance and strike time deviation meet the requirements, and the required overload is within the available overload range, which shows the effectiveness of this method.

4.2. Performance Verification of Three-Dimensional Cooperative Guidance Law under Time-Varying Velocity

In order to verify the three-dimensional cooperative guidance law based on desired impact time and the three-dimensional cooperative guidance law based on coordination variables proposed in this paper, the following simulations are carried out. The position of the target is (200, 0, 0) km, and the four UAVs are required to reach the target position at the same time to strike the target. The proportional guidance constants are all 3, and the initial parameters are shown in

Table 3. Simulation parameters under the three-dimensional cooperative guidance law.

UAVs	Initial Position (km)	Velocity (m/s)	Initial Flight Path Angle (deg)	Initial Heading Angle (deg)
M1	(−35, 20, −30)	1950	0	0
M2	(−15, 21, −17)	1850	2	−2
M3	(15, 25, 20)	1750	1	2
M4	(25, 22, 31)	1800	−1	3

.

4.2.1. Simulations of Three-Dimensional Cooperative Guidance Law Based on Desired Impact Time When the Rate of Change of Velocity Is the First-Order Function of Velocity

The desired impact time is set to 170 s, the simulation results are shown in Figure 6.

As can be seen from Figure 6a–e, under the action of the cooperative guidance law, the four UAVs gradually approach and finally strike the target. From Figure 6b–e, it can be seen that the time for the UAVs to reach the target is 170 s, which is the same as the desired impact time, which indicates the effectiveness of the cooperative guidance law. As can be seen from Figure 6b, the impact time errors of the four UAVs gradually decrease, which reach the same in about 80 s, and converge to zero in the terminal time. As can be seen from Figure 6e, the leading angle in the yaw channel increases in the early stage and then decreases gradually until the terminal time converges to 0.

4.2.2. SIMULATIONS of Three-Dimensional Cooperative Guidance Law Based on Coordinated Variables When the Rate of Change of Velocity Is the First-Order Function of Velocity

Considering the mutual communication, the desired impact time is negotiated among UAVs. The simulation results are shown in Figure 7.

Figure 7a–f show the variation curves of three-dimensional flight trajectory, impact time error, flight path angle, heading angle, the leading angle of velocity in yaw channel and the weight with time respectively. As can be seen from Figure 7a, all four UAVs can reach the specified target point. As can be seen from Figure 7b, the impact time errors of the four UAVs converge to 0 at terminal time, and the desired impact time obtained through negotiation is 188.2 s. As can be seen from Figure 7f, M1 has the largest weight, so this member plays a decisive role in the negotiation process. The remaining three UAVs adjust their flight trajectories according to the desired impact time determined by M1, so as to achieve cooperative strike target. As can be seen from Figure 7b, the impact time error of M1 first converges to 0.

4.2.3. Simulations of Three-Dimensional Cooperative Guidance Law Based on Coordinated Variables When the Rate of Change of Velocity Is the Quadratic Function of Velocity

Simulations of three-dimensional cooperative guidance law based on coordinated variables when the rate of change of velocity is the quadratic function of velocity are shown in Figure 8. As can be seen from Figure 8a–f, the four UAVs hit the target at the same time through negotiation. Compared with Figure 7, the flight time is shorter, only 154.1 s, this is because when the velocity change rate is a quadratic function of velocity, the velocity decreases more slowly. As can be seen from Figure 8f, M1 has the largest weight, so this member plays a decisive role in the negotiation process. The remaining three UAVs adjust their flight trajectories according to the desired impact time determined by M1, so as to achieve cooperative strike target. As can be seen from Figure 8b, the impact time error of M1 first converges to 0.

To summarize, the three-dimensional cooperative guidance law with impact time constraint under time-varying velocity proposed in this paper is effective and can make multiple UAVs attack targets cooperatively.

5. Conclusions

This paper has studied the terminal multiple UAVs cooperative guidance problem with impact time control under time-varying velocity. Firstly, on a two-dimensional plane, the cooperative guidance problem is transformed into a second-order consistency control problem, considering the finite time consistency of impact time, the time-varying velocity and the distributed communication topology, a distributed time cooperative guidance law based on graph theory and consistency theory is designed. Then the cooperative guidance problem is expanded from two-dimensional plane to three-dimensional space, compared with reference [34], when the rate of change of velocity is a quadratic function of velocity, a three-dimensional impact time cooperative guidance law based on coordination variables is proposed in this paper, which does not need to set the desired impact time in advance; the impact time is determined by negotiation among UAVs. When the rate of change of velocity is the first order function of velocity, the expression of time-to-go estimation is derived firstly, then according to whether there is the communication among UAVs, a multiple UAVs three-dimensional cooperative guidance law based on desired impact time and a multiple UAVs three-dimensional cooperative guidance law based on coordination variable are designed, respectively. The simulation results show that the cooperative guidance laws designed in this paper can all realize the saturation attack.