A DMPC-based three-dimensional cooperative guidance scheme with impact time and impact angle constraints

Abstract

In this paper, a distributed model predictive control (DMPC)-based three-dimensional (3D) cooperative guidance scheme is proposed for multiple interceptors with controllable thrust and constrained acceleration, to simultaneously attack a maneuvering target at desired terminal angles. The cooperative guidance scheme consists of a fractional power function-based extended state observer (FPESO) for disturbance estimation, and a DMPC controller to generate the control actions. By introducing the assumed state trajectory, the DMPC controller (of each interceptor) can simultaneously employ the zebra optimization algorithm (ZOA) to solve the constrained optimization problem, and generate the real-time guidance commands with no need of any global information. For restricting the uncertain deviation between the assumed state trajectory and the actual one, a compatibility constraint is designed, and the terminal ingredients are developed for coupling and decoupling stage costs, respectively. By applying the above designed constraints and the terminal ingredients, the recursive feasibility of the resulting optimization problem and the closed-loop stability of multi-interceptor system are both ensured. Through numerical simulations, it is verified that the proposed cooperative guidance scheme can achieve cooperative interception with impact time and angle constraints in various scenarios, and has advantageous performance.

Introduction

With the improvement of target maneuverability and the advancement of anti-interception techniques, the traditional single interceptor-based operation is no longer suitable for the complex and ever-changing modern battlefield environment. Owing to the ability to increase the success rate of interception, the salvo attack of multiple interceptors is considered a more effective combat pattern, thus the cooperative guidance laws have received increasing concern.^1–5

Early cooperative guidance laws mainly focused on the simultaneous arrival of interceptors, which was achieved through two main methods. The first method is by designating an impact time, then controlling all the interceptors to hit the target exactly at that time.^6–8 This sort of control methods is called impact time control guidance (ITCG) laws, and has long attracted interests of researchers. However, the inherent weakness of ITCG laws cannot be evaded. As pointed out by Zhao et al., ⁹ due to the absence of communication among interceptors during guidance, the ICTG law cannot be viewed as a genuine cooperative guidance law. Besides, an inappropriate impact time assigned in advance might lead to the failure of interception. For these considerations, extensive research has been conducted on the second method of achieving simultaneous arrival, which is by designing a specific coordination algorithm utilizing the communication information among interceptors.^9–13 As the pioneering research of this way, Zhao et al. ⁹ presented a cooperative guidance scheme which combined the coordination algorithms with local guidance laws, where the coordination algorithm was designed using multi-agent consensus theory, serving as an additional term to local guidance laws such as PNG. Due to the utilization of the local guidance laws, the performance of the above cooperative guidance law is guaranteed. However, since the assumption that the velocities of interceptors are constants are implicit in the above cooperative guidance law, ¹² it may not be applicable in scenarios where the magnitude of velocity of interceptors is variable. Different from the above work, Lyu et al. ¹³ developed a cooperative guidance law with two components in the planar line-of-sight (LOS) frame, in which the component along the direction of LOS guaranteed all interceptors’ impact times achieved consensus asymptotically, while the component along vertical direction of LOS ensured all the interceptors precisely hit the target. This design method is applicable to interceptors with controllable thrust, and adopted by most research on impact angle-constrained guidance law. ¹⁴

Recent cooperative guidance laws focused more on the specific requirements in practical missions, such as impact angle control (or impact angle constraints),^14–19 because cooperative interception at different impact angles is highly important for maximizing the damage effect against high value targets. ¹⁴ Adopting the design method from Lyu, ¹³ Zhang et al. ¹⁴ proposed the cooperative guidance laws in the LOS direction and in the normal direction to the LOS, respectively. Among them, the former was based on the finite-time consensus to achieve the simultaneous attack of all the interceptors, the latter was derived using the super-twisting control algorithm, to make the impact angle of each interceptor converge to the expected value. Song et al. ¹⁸ and Dong et al. ¹⁹ further extended such a design method to three-dimensional (3D) interception scenario, using adaptive super-twisting sliding mode control, integral sliding mode control, fixed-time consensus and other techniques to achieve 3D cooperative interception with impact time and angle constraints.

Expect for the impact angle constraint, the acceleration saturation constraint is another factor that should be considered from the practical viewpoint, since the acceleration that can be provided is restricted by the control ability of actuator. ²⁰ With the natural advantage of handling various constraints, model predictive control (MPC) is a preferred approach in distributed control of multi-agent system,^21–23 and has been applied to multi-interceptor cooperative guidance with acceleration saturation constraint in recent years.^12,24–27 The most significant advantage of MPC over traditional optimal control is that it employs the optimization procedure online to determine the control actions, rather than using a pre-computed control law. This endows MPC with the ability to handle external disturbances and uncertainty. Among the existing works that applied MPC to cooperative guidance problems, Zhao et al. ²⁴ proposed an impact time-constrained cooperative guidance approach using nonlinear MPC, and introduced the assumed (estimated) state trajectory to make interceptors synchronously solve the optimization problems subject to control input constraints. Jiang et al. ¹² employed MPC and convex optimization technique to solve for the optimal guidance gain of each interceptor, thereby achieving the consensus on times-to-go among interceptors under the constraint on the maximum lateral acceleration. Cong et al. ²⁷ used distributed model predictive control (DMPC) to design the cooperative guidance algorithm for follower interceptors under virtual leader-follower cooperative strategy, and introduced the assumed state trajectory. Although the above works provide evidence that MPC (DMPC) can be well used for cooperative guidance with multiple constraints, the following issues have not yet been addressed for most studies on DMPC-based cooperative guidance.

• (a) In several works,^24,25,27 the assumed state trajectory was introduced to make interceptors synchronously compute their guidance commands, however, none of the above studies have considered the uncertain deviation between the assumed state trajectory and the actual one, which might affect the stability of the overall system.

• (b) Although the performance of DMPC has been validated through simulations or experiments,^12,24–27 the recursive feasibility and closed-loop stability of DMPC have not been theoretically proven in any cooperative guidance scheme, and few of studies have formulated the associated constraints to guarantee the above properties.

• (c) Up to now, the research that simultaneously incorporate the acceleration saturation constraint, impact time, and impact angle control into the cooperative guidance law design is still little. Besides, most of studies based on MPC only concerned the cooperative guidance against stationary or low-speed targets.

In view of the above issues, a DMPC-based 3D cooperative guidance scheme is worthy of further investigation, especially considering the huge potential of MPC in integrating and handling multiple objectives and constraints. These objectives and constraints typically reflect the practical requirements of combat, and have great significance in maximizing the success rates of cooperative interception and enhancing the combat capability of interceptors in different engagement scenarios. Therefore, this paper proposes a DMPC-based 3D cooperative guidance scheme for intercepting a maneuvering target, which is able to simultaneously achieve impact time and impact angle control under the acceleration saturation constraints. The main contributions of this paper are highlighted as follows.

• 1) To the best of our knowledge, it is the first time that DMPC is used to solve 3D cooperative guidance problem with impact time and impact angle constraints for intercepting highly maneuvering targets. Different from existing methods for the above problem, we carry out the controller design under a unified framework of DMPC, rather than decomposing the guidance law into several parts.^13,14,18,19 Besides, the acceleration saturation constraint is also taken into our consideration from the practical perspective.

• 2) Compared with studies^24,25,27 which only applied synchronous DMPC to cooperative guidance problems, we further formulate the associated constraints to ensure the closed-loop stability of the multi-interceptor system under DMPC, without utilizing any globally shared reference states.^21–23 Different from the method ²⁶ which integrally generates the guidance commands for all the interceptors, we adopt zebra optimization algorithm (ZOA) to solve for the optimal guidance commands independently for each interceptor.

The remainder of this paper is organized as follows. In section 2, some preliminaries about 3D cooperative guidance are provided, and the nonlinear dynamics for each subsystem are established. In section 3, the detailed design of the distributed cooperative guidance scheme is conducted, and the complete implement process is introduced. In section 4, the properties of the cooperative guidance scheme are analyzed. In section 5, the numerical simulations are carried out under diverse cases, in which the proposed cooperative guidance scheme is compared with other advanced methods. Finally, in section 6, the conclusions are drawn.

Preliminaries

Consider the guidance scenario of a group of interceptors attacking a single maneuvering target. The 3D engagement geometry between the $i$ th interceptor and the target is depicted in Figure 1, where ${\mathrm{M}}_i$ and $\mathrm{T}$ denote the $i$ th interceptor and the target, $i=1,2,\cdots ,N_{\mathrm{M}}$ , ${\mathrm{M}}_i\mathrm{XYZ}$ and ${\mathrm{M}}_i{X}_{\mathrm{L}}{Y}_{\mathrm{L}}{Z}_{\mathrm{L}}$ represent the inertial reference frame and the line-of-sight (LOS) reference frame, respectively. The velocity of ${\mathrm{M}}_i$ is denoted as $v_{{\mathrm{M}}_i}$ . ${\theta }_{{\mathrm{M}}_i}$ and ${\psi }_{{\mathrm{M}}_i}$ represent the pitch and yaw path angles of ${\mathrm{M}}_i$ . The interceptors are assumed to be thrust controllable, that is, the acceleration of ${\mathrm{M}}_i$ can be decomposed into three mutually orthogonal components ${[a_{\mathrm{M}{x}_i},a_{\mathrm{M}{y}_i},a_{\mathrm{M}{z}_i}]}^{\mathrm{T}}$ in inertial frame, and the kinematic equations of ${\mathrm{M}}_i$ can be expressed as equation (1).

$$\{\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{x}}_{{\mathrm{M}}_i}=v_{{\mathrm{M}}_i}\cos {\theta }_{{\mathrm{M}}_i}\cos {\psi }_{{\mathrm{M}}_i}}\\ \multicolumn{1}{c}{{\stackrel{·}{y}}_{{\mathrm{M}}_i}=v_{{\mathrm{M}}_i}\sin {\theta }_{{\mathrm{M}}_i}}\\ \multicolumn{1}{c}{{\stackrel{·}{z}}_{{\mathrm{M}}_i}=-v_{{\mathrm{M}}_i}\cos {\theta }_{{\mathrm{M}}_i}\sin {\psi }_{{\mathrm{M}}_i}}\\ \multicolumn{1}{c}{{\stackrel{·}{v}}_{{\mathrm{M}}_i}=a_{\mathrm{M}{x}_i}}\\ \multicolumn{1}{c}{{\stackrel{·}{\theta }}_{{\mathrm{M}}_i}=a_{\mathrm{M}{y}_i}/v_{{\mathrm{M}}_i}}\\ \multicolumn{1}{c}{{\stackrel{·}{\psi }}_{{\mathrm{M}}_i}=-a_{\mathrm{M}{z}_i}/\left(v_{{\mathrm{M}}_i}\cos {\theta }_{{\mathrm{M}}_i}\right)}\end{array}$$

(1)

Figure 1.

Three dimensional engagement geometry between the interceptor and the target.

Similarly, the motion of the target in inertial frame is obtained as equation (2), the meanings of $v_{\mathrm{T}}$ , ${\theta }_{\mathrm{T}}$ , ${\psi }_{\mathrm{T}}$ and ${[a_{\mathrm{T}x},a_{\mathrm{T}y},a_{\mathrm{T}z}]}^{\mathrm{T}}$ are identical to those of ${\mathrm{M}}_i$ .

$$\{\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{x}}_{\mathrm{T}}=v_{\mathrm{T}}\cos {\theta }_{\mathrm{T}}\cos {\psi }_{\mathrm{T}}}\\ \multicolumn{1}{c}{{\stackrel{·}{y}}_{\mathrm{T}}=v_{\mathrm{T}}\sin {\theta }_{\mathrm{T}}}\\ \multicolumn{1}{c}{{\stackrel{·}{z}}_{\mathrm{T}}=-v_{\mathrm{T}}\cos {\theta }_{\mathrm{T}}\sin {\psi }_{\mathrm{T}}}\\ \multicolumn{1}{c}{{\stackrel{·}{v}}_{\mathrm{T}}=a_{\mathrm{T}x}}\\ \multicolumn{1}{c}{{\stackrel{·}{\theta }}_{\mathrm{T}}=a_{\mathrm{T}y}/v_{\mathrm{T}}}\\ \multicolumn{1}{c}{{\stackrel{·}{\psi }}_{\mathrm{T}}=-a_{\mathrm{T}z}/\left(v_{\mathrm{T}}\cos {\theta }_{\mathrm{T}}\right)}\end{array}$$

(2)

The relative motion between ${\mathrm{M}}_i$ and $\mathrm{T}$ is described using differential equations (3)–(5), where $r_i$ , ${\theta }_{{\mathrm{L}}_i}$ , ${\psi }_{{\mathrm{L}}_i}$ are the relative distance, the elevation, and azimuth angles of the LOS to the inertial reference. ${[a_{\mathrm{M}{r}_i},a_{\mathrm{M}{\theta }_i},a_{\mathrm{M}{\psi }_i}]}^{\mathrm{T}}$ and ${[a_{\mathrm{T}{r}_i},a_{\mathrm{T}{\theta }_i},a_{\mathrm{T}{\psi }_i}]}^{\mathrm{T}}$ represent the acceleration vectors of ${\mathrm{M}}_i$ and T in LOS frame, both of them can be transformed into ${[a_{\mathrm{M}{x}_i},a_{\mathrm{M}{y}_i},a_{\mathrm{M}{z}_i}]}^{\mathrm{T}}$ or ${[a_{\mathrm{T}x},a_{\mathrm{T}y},a_{\mathrm{T}z}]}^{\mathrm{T}}$ utilizing transformation matrices.

$${\stackrel{··}{r}}_i-r_i{\stackrel{·}{\theta }}_{{\mathrm{L}}_i}^2-r_i{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}^2{\cos }^2{\theta }_{{\mathrm{L}}_i}=a_{\mathrm{T}{r}_i}-a_{\mathrm{M}{r}_i}$$

(3)

$$r_i{\stackrel{··}{\theta }}_{{\mathrm{L}}_i}+2{\stackrel{·}{r}}_i{\stackrel{·}{\theta }}_{{\mathrm{L}}_i}+r_i{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}^2\sin {\theta }_{{\mathrm{L}}_i}\cos {\theta }_{{\mathrm{L}}_i}=a_{\mathrm{T}{\theta }_i}-a_{\mathrm{M}{\theta }_i}$$

(4)

$$\begin{array}{c}\multicolumn{1}{c}{-r_i{\stackrel{··}{\psi }}_{{\mathrm{L}}_i}\cos {\theta }_{{\mathrm{L}}_i}-2{\stackrel{·}{r}}_i{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}\cos {\theta }_{{\mathrm{L}}_i}+2r_i{\stackrel{·}{\theta }}_{{\mathrm{L}}_i}{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}\sin {\theta }_{{\mathrm{L}}_i}}\\ \multicolumn{1}{c}{=a_{\mathrm{T}{\psi }_i}-a_{\mathrm{M}{\psi }_i}}\end{array}$$

(5)

Commonly, the guidance problem with impact angle constraints is also deemed as the control problem of line-of-sight angles, ¹⁸ that is, controlling the LOS angles of each interceptor to converge to their desired values at the end of guidance.

$${\theta }_{{\mathrm{L}}_i}\left(t_{{\mathrm{f}}_i}\right)={\theta }_{\mathrm{L}{\mathrm{d}}_i},{\psi }_{{\mathrm{L}}_i}\left(t_{{\mathrm{f}}_i}\right)={\psi }_{\mathrm{L}{\mathrm{d}}_i}$$

(6)

In equation (6), $t_{{\mathrm{f}}_i}$ denotes the terminal time of guidance (i.e. impact time), which is defined as equation (7). ${\theta }_{\mathrm{L}{\mathrm{d}}_i}$ , ${\psi }_{\mathrm{L}{\mathrm{d}}_i}$ are the desired terminal LOS angles.

$$t_{{\mathrm{f}}_i}=t+t_{\mathrm{g}{\mathrm{o}}_i}$$

(7)

In equation (7), $t_{\mathrm{g}{\mathrm{o}}_i}$ represents the remaining time before interception (also known as time-to-go), which can be estimated by

$$t_{\mathrm{g}{\mathrm{o}}_i}=-r_i/{\stackrel{·}{r}}_i$$

(8)

According to equations (7) and (8), coordinating the impact time among interceptors is equivalent to achieving the consensus on the remaining time $t_{\mathrm{g}{\mathrm{o}}_i}$ , so the cooperative guidance problem with impact time and impact angle constraints can be described as below.

• Definition 1: The multi-interceptor system is said to achieve cooperative interception, if the following equations hold for arbitrary $i,j\in \{1,2,\cdots ,N_{\mathrm{M}}\}$ as time $t$ approaches $t_{{\mathrm{f}}_i}$ .

$$t_{\mathrm{g}{\mathrm{o}}_i}=t_{\mathrm{g}{\mathrm{o}}_j}$$

(9)

$${\theta }_{{\mathrm{L}}_i}={\theta }_{\mathrm{L}{\mathrm{d}}_i},{\psi }_{{\mathrm{L}}_i}={\psi }_{\mathrm{L}{\mathrm{d}}_i}$$

(10)

$${\stackrel{·}{\theta }}_{{\mathrm{L}}_i}=0,{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}=0$$

(11)

• Remark 1: According to the parallel approach principle, the interceptor ${\mathrm{M}}_i$ can hit the target with a sufficiently small miss distance as long as the LOS angular rates ${\stackrel{·}{\theta }}_{{\mathrm{L}}_i}$ , ${\stackrel{·}{\psi }}_{{\mathrm{L}}_i}$ converge and stabilize at the vicinity of zero, which provides the necessity for equation (11).

According to definition 1 and equations (3)–(5), the system dynamics of each interceptor are established as follows.

$${\stackrel{·}{\mathit{x}}}_i=f({\mathit{x}}_i,{\mathit{u}}_i,{\mathit{w}}_i)$$

(12)

where

$$\begin{array}{c}\multicolumn{1}{c}{f({\mathit{x}}_i,{\mathit{u}}_i,{\mathit{w}}_i)}\\ \multicolumn{1}{c}{=\left[\begin{array}{c}-x_{i1}/x_{i2}\hfill \\ \multicolumn{1}{c}{-1+x_{i2}^2{x}_{i4}^2+x_{i2}^2{x}_{i6}^2{\cos }^2(x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\\ \multicolumn{1}{c}{x_{i4}}\\ \multicolumn{1}{c}{2x_{i4}/x_{i2}-x_{i6}^2\sin (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})\cos (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\\ \multicolumn{1}{c}{x_{i6}}\\ \multicolumn{1}{c}{2x_{i6}/x_{i2}+2x_{i4}{x}_{i6}\tan (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\end{array}\right]}\\ \multicolumn{1}{c}{+\left({\mathit{B}}_i({\mathit{u}}_i-{\mathit{w}}_i)\right)\otimes \left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\right]}\end{array}$$

(13)

$${\mathit{x}}_i={[r_i,t_{\mathrm{g}{\mathrm{o}}_i},{\theta }_{{\mathrm{L}}_i}-{\theta }_{\mathrm{L}{\mathrm{d}}_i},{\stackrel{·}{\theta }}_{{\mathrm{L}}_i},{\psi }_{{\mathrm{L}}_i}-{\psi }_{\mathrm{L}{\mathrm{d}}_i},{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}]}^{\mathrm{T}}$$

(14)

$${\mathit{B}}_i=\mathrm{diag}\{-x_{i2}^2/x_{i1},-1/x_{i1},1/\left(x_{i1}\cos (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})\right)\}$$

(15)

$${\mathit{u}}_i={[a_{\mathrm{M}{r}_i},a_{\mathrm{M}{\theta }_i},a_{\mathrm{M}{\psi }_i}]}^{\mathrm{T}}$$

(16)

$${\mathit{w}}_i={[a_{\mathrm{T}{r}_i},a_{\mathrm{T}{\theta }_i},a_{\mathrm{T}{\psi }_i}]}^{\mathrm{T}}$$

(17)

In the above equations, ${\mathit{u}}_i$ denotes the control input signal to be designed, which is constrained by

$${\mathit{u}}_i\in {\mathbb{U}}_i$$

(18)

$${\mathbb{U}}_i=\left\{{[u_{i1},u_{i2},u_{i3}]}^{\mathrm{T}}|{\underset{¯}{u}}_{\mathrm{ij}}\le u_{\mathrm{ij}}\le {\overline{u}}_{\mathrm{ij}},j=1,2,3\right\}$$

(19)

${\mathit{w}}_i$ represents the disturbance caused by the target maneuvering, and the nominal system of (12) is given as equation (20), by neglecting the disturbances ${\mathit{w}}_i$ . ⊗ is the operator of Kronecker product.

$${\stackrel{·}{\stackrel{~}{\mathit{x}}}}_i=f({\tilde{\mathit{x}}}_i,{\tilde{\mathit{u}}}_i)$$

(20)

Based on the system formulated above, we will design a guidance scheme which consists of the control inputs ${\tilde{\mathit{u}}}_i$ for the nominal system and the estimated disturbances ${\hat{\mathit{w}}}_i$ to offset the effect induced by disturbances. The proposed scheme is given by

$${\mathit{u}}_i={\tilde{\mathit{u}}}_i+{\hat{\mathit{w}}}_i$$

(21)

In order to facilitate the guidance scheme design, the following assumptions are necessary.

• Assumption 1: Suppose that the system disturbances are second-order differentiable and satisfy $\left|w_{\mathrm{ij}}^{\left[l\right]}\right|\le {\omega }_{\mathrm{ij}}$ , for any $i=1,2,\cdots ,N_{\mathrm{M}}$ , $j=1,2,3$ , and $l=0,1,2$ , where $w_{\mathrm{ij}}^{\left[l\right]}$ is the $l$ -th order differentiative of $w_{\mathrm{ij}}$ with respect to time, ${\omega }_{\mathrm{ij}}$ is a positive constant.

• Assumption 2: Throughout the process of guidance, the relative velocity between the target and the interceptor is always negative, that is, ${\stackrel{·}{r}}_i<0$ . Besides, the relative distance when the interception occurs is a small positive constant, that is, $r_i\left(t_{{\mathrm{f}}_i}\right)>0$ , thus the time-to-go of each interceptor satisfies $t_{\mathrm{g}{\mathrm{o}}_i}>0$ for any $t\le t_{{\mathrm{f}}_i}$ .

• Assumption 3: the relative distance $r_i$ , the relative velocity between the target and the interceptor ${\stackrel{·}{r}}_i$ , LOS angles ${\theta }_{{\mathrm{L}}_i}$ , ${\psi }_{{\mathrm{L}}_i}$ and LOS angular rates ${\stackrel{·}{\theta }}_{{\mathrm{L}}_i}$ , ${\stackrel{·}{\psi }}_{{\mathrm{L}}_i}$ can be acquired by each interceptor in real time.

• Remark 2: the above assumptions are reasonable from the practical viewpoint, because

• (1) For the target, the control actuator requires a process to implement the maneuvering command, thus the variation of the maneuvering accelerations is continuous and smooth.

• (2) When the relative velocity changes from negative to positive, the target either escapes interception or is intercepted, the guidance will immediately terminate no matter which of the above situation occurs, thus the relative velocity is always negative during guidance.

• (3) For the interceptor equipped with an active radar seeker, the aforementioned information can be acquired or calculated in real time.

Design and implementation of the distributed cooperative guidance scheme

For a strongly coupled and highly nonlinear system as shown in (13), it is challenging to develop a distributed control law to achieve the cooperative interception as required in definition 1, especially when the system is subjected to unknown disturbances. However, by adopting the guidance scheme given by (21), we can simultaneously achieve the disturbance rejection and the stabilization of the multi-interceptor system. This motivates us to design a DMPC-based cooperative guidance scheme, as shown in Figure 2.

Figure 2.

The control loop of the cooperative guidance scheme.

The cooperative guidance scheme consists of a fractional power function-based extended state observer (FPESO) for disturbance estimation, and a DMPC controller to generate the control actions. By the virtue of DMPC, we can handle all the constraints involved in the cooperative guidance problem. The focus of controller design is on the development of compatibility constraints and the terminal ingredients, as they are all crucial guarantees for the stability of the multi-interceptor system, and have not yet been designed in relevant research. Besides, we also aim to enhance the real-time performance of the cooperative guidance scheme by efficiently solving the DMPC problem.

Design of the extended state observer

The extended state observer is central to the disturbance rejection control. Through a properly designed ESO, the disturbances can be estimated and canceled in the feedback loop. ²⁸ In the recent studies regarding ESO, Han ²⁹ proposed a nonlinear ESO based on fractional power function (FPESO), Zhao ²⁸ illustrated the relationship between the selection of the gains and the convergence of the aforementioned FPESO, and Zhang ³⁰ designed a FPESO to improve the estimation accuracy for the compound disturbances. These works demonstrate the advantage of FPESO over conventional linear ESO in terms of noise tolerance. Motivated by them, in this subsection, a FPESO is designed to estimate disturbances caused by target maneuver.

First, the states to be observed are selected as follows.

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\chi }}_i={[{\chi }_{i1},{\chi }_{i2},{\chi }_{i3}]}^{\mathrm{T}}={[{\stackrel{·}{r}}_i,r_i{\stackrel{·}{\theta }}_{{\mathrm{L}}_i},-r_i{\stackrel{·}{\psi }}_{{\mathrm{L}}_i}\cos {\theta }_{{\mathrm{L}}_i}]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{={[-x_{i1}/x_{i2},x_{i1}{x}_{i4},-x_{i1}{x}_{i6}\cos (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})]}^{\mathrm{T}}}\end{array}$$

(22)

Differentiating ${\mathbf{\chi }}_i$ with respect to time, yields

$${\stackrel{·}{\chi }}_i={\mathit{G}}_i+{\mathit{w}}_i$$

(23)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{G}}_i=}\\ \multicolumn{1}{c}{\left[\begin{array}{c}\multicolumn{1}{c}{x_{i1}{x}_{i4}^2+x_{i1}{x}_{i6}^2{\cos }^2(x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\\ \multicolumn{1}{c}{x_{i1}{x}_{i4}/x_{i2}-x_{i1}{x}_{i6}^2\sin (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})\cos (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\\ \multicolumn{1}{c}{-x_{i1}{x}_{i6}\cos (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})/x_{i2}-x_{i1}{x}_{i4}{x}_{i6}\sin (x_{i3}+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\end{array}\right]-{\mathit{u}}_i}\end{array}$$

(24)

Then, given assumption 1 and inspired by Zhao, ²⁸ we design the FPESO as follows.

$$\{\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\stackrel{·}{\mathit{z}}}}_{i1}={\mathit{G}}_i+{\mathit{z}}_{2i}+{\beta }_1\mathrm{fal}({\alpha }^2{\mathit{e}}_{i1},{\vartheta }_1,\varphi )/\alpha }\\ \multicolumn{1}{c}{{\stackrel{·}{\stackrel{·}{\mathit{z}}}}_{i2}={\mathit{z}}_{3i}+{\beta }_2\mathrm{fal}({\alpha }^2{\mathit{e}}_{i1},{\vartheta }_2,\varphi )}\\ \multicolumn{1}{c}{{\stackrel{·}{\stackrel{·}{\mathit{z}}}}_{i3}=\alpha {\beta }_3\mathrm{fal}({\alpha }^2{\mathit{e}}_{i1},{\vartheta }_3,\varphi )}\\ \multicolumn{1}{c}{{\mathit{e}}_{i1}={[e_{i1}^{(1)},e_{i1}^{(2)},e_{i1}^{(3)}]}^{\mathrm{T}}={\mathbf{\chi }}_i-{\mathit{z}}_{i1}}\end{array}$$

(25)

In (25), ${\mathit{z}}_{i1}$ , ${\mathit{z}}_{i2}$ , ${\mathit{z}}_{i3}$ are the estimates of ${\mathbf{\chi }}_i$ , ${\mathit{w}}_i$ and ${\stackrel{·}{\mathit{w}}}_i$ , respectively. $\mathrm{fal}(•)$ is the fractional power function, with $\mathrm{fal}({\alpha }^2{\mathit{e}}_{i1},{\vartheta }_j,\varphi )$ , $j=1,2$ being

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{fal}({\alpha }^2{\mathit{e}}_{i1},{\vartheta }_j,\varphi )=\left[\begin{array}{c}\multicolumn{1}{c}{\mathrm{fal}({\alpha }^2{e}_{i1}^{\left(1\right)},{\vartheta }_j,\varphi )}\\ \multicolumn{1}{c}{\mathrm{fal}({\alpha }^2{e}_{i1}^{\left(2\right)},{\vartheta }_j,\varphi )}\\ \multicolumn{1}{c}{\mathrm{fal}({\alpha }^2{e}_{i1}^{\left(3\right)},{\vartheta }_j,\varphi )}\end{array}\right]}\end{array}$$

(26)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{fal}({\alpha }^2{e}_{i1}^{\left(k\right)},{\vartheta }_j,\varphi )=}\\ \multicolumn{1}{c}{\{\begin{array}{c}{\left|{\alpha }^2{e}_{i1}^{\left(k\right)}\right|}^{{\vartheta }_j}\mathrm{sign}\left({\alpha }^2{e}_{i1}^{\left(k\right)}\right),\left|{\alpha }^2{e}_{i1}^{\left(k\right)}\right|>\varphi \hfill \\ \multicolumn{1}{c}{{\alpha }^2{e}_{i1}^{\left(k\right)}/\left({\varphi }^{1-{\vartheta }_j}\right),\mathrm{else}}\end{array},k=1,2,3}\end{array}$$

(27)

parameter $\alpha$ is a constant gain, ${\beta }_1$ , ${\beta }_2$ , ${\beta }_3$ are constants to be chosen so that the following matrix is Hurwitz.

$$\left[\begin{array}{ccc}-{\beta }_1\hfill & 1& 0\\ \multicolumn{1}{c}{-{\beta }_2}& 0& 1\\ \multicolumn{1}{c}{-{\beta }_3}& 0& 0\end{array}\right]$$

The positive constants ${\vartheta }_j\in (0,1)$ are determined by equation ${\vartheta }_j=j\vartheta -(j-1)$ , $j=1,2,3$ , where $\vartheta$ is a positive constant that satisfies $2/3<\vartheta <1$ .

Let the estimation errors be ${\mathit{e}}_{i1}$ , ${\mathit{e}}_{i2}$ , and ${\mathit{e}}_{i3}$ , where ${\mathit{e}}_{i2}={\mathit{w}}_i-{\mathit{z}}_{i2}$ , ${\mathit{e}}_{i3}={\stackrel{·}{\mathit{w}}}_i-{\mathit{z}}_{i3}$ , then the convergence of the FPESO is presented by lemma 1.

• Lemma 1 ²⁸ : Suppose that the disturbances ${\mathit{w}}_i$ and their $l$ -th order derivatives are uniformly bounded, $l=1,2$ . Then, for the FPESO constructed as (25), there exist ${\vartheta }^{*}\in (2/3,1)$ and ${\alpha }^{*}>0$ such that, for any $\vartheta \in [{\vartheta }^{*},1)$ , $\alpha >{\alpha }^{*}$ , any initial ${\mathbf{\chi }}_i\left(t_0\right)$ , ${\mathit{w}}_i\left(t_0\right)$ , ${\stackrel{·}{\mathit{w}}}_i\left(t_0\right)$ , and initial states of FPESO ${\mathit{z}}_{i1}\left(t_0\right)$ , ${\mathit{z}}_{i2}\left(t_0\right)$ , ${\mathit{z}}_{i3}\left(t_0\right)$ , the estimation errors satisfy, for any $t>t_{\alpha }$ , that

$$\left|e_{\mathrm{ij}}^{\left(k\right)}\right|\le \Gamma {\alpha }^{j-4},k=1,2,3,j=1,2,3$$

(28)

where $t_{\alpha }$ is $\alpha$ -dependent and satisfies $\underset{\alpha \to \infty }{lim}{t}_{\alpha }=0$ , $\Gamma$ is a constant independent of $\alpha$ .

• Remark 3: Equation (28) together with $\underset{\alpha \to \infty }{lim}{t}_{\alpha }=0$ implies that, for any $k=1,2,3$ and $j=1,2,3$ , $\underset{\alpha \to \infty }{lim}\underset{t>0}{sup}\left|e_{\mathrm{ij}}^{\left(k\right)}\left(t\right)\right|=0$ . Generally speaking, as long as the gain $\alpha$ is tuned according to the variation speed of the disturbances, the estimates of ${\mathit{w}}_i$ are able to approximate the actual disturbances to an acceptable small error. That is, when $t>t_{\alpha }$ , the impact of estimation errors on the multi-interceptor system is negligible.

Based on lemma 1, the disturbances caused by target maneuvering are completely compensated by implementing control inputs (21) and the FPESO (25). Therefore, the closed-loop stability of the multi-interceptor system is guaranteed so long as the nominal control inputs are appropriately designed. In the next subsection, a cooperative guidance law based on DMPC will be proposed to guarantee the stability of the multi-interceptor system, thereby achieving all the goals in definition 1.

Design of the DMPC controller

As an optimization-based control technique, MPC uses an optimization problem over a given horizon to describe what system behavior is expected and which constraints must be satisfied. ³¹ With the up-to-date information of the system being acquired at each sampling instant, such an optimization problem will be renewed and solved by employing an optimization procedure, and a portion of the optimal solution will be implemented as the real-time control action. This pattern provides MPC with advantages in handling complex constraints and disturbances rejection.

As an emerging control strategy of MPC, DMPC is characterized by lower computational burden and more flexible structure, as it does not require a central agent to compute the control actions for, and communicate with all remaining agents. Among several update strategies of DMPC (including sequential DMPC, ³² iterative DMPC, ³³ and synchronous DMPC^21–23), the synchronous DMPC allows all the agents to synchronously solve the optimization problem and obtain the control inputs, which ensures the real-time performance of the overall system.

Therefore, the synchronous DMPC will be used to design the cooperative guidance scheme in this subsection. By implementing the synchronous DMPC, each interceptor only needs to communicate with its neighbors and independently solve its optimization problem. Since each neighbor is coupled with the interceptor through a coupling cost function, the main focus is to handle the coupling terms by designing associated constraints and terminal ingredients, thereby ensuring the recursive feasibility of the optimization problem and stability of the overall system.

Basic model of DMPC

Considering the control objectives given in definition 1, a finite-horizon optimization control problem (FHOCP) for the nominal system (20) is defined as follows.

• Problem 1: At the time instant $t_k$ , given the system state ${\mathit{x}}_i\left(t_k\right)$ , the (nominal) state trajectories of neighbors ${\tilde{\mathit{x}}}_{i-}\left(\tau |t_k\right)$ over the time interval $\tau \in [t_k,t_k+T]$ , find the optimal control inputs

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)=\underset{{\tilde{\mathit{u}}}_i\left(\tau |t_k\right)}{\arg min}{\int }_{t_k}^{t_k+T}{L}_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-},{\tilde{\mathit{u}}}_i)\mathrm{d}\tau }\\ \multicolumn{1}{c}{+g_i\left({\tilde{\mathit{x}}}_i(t_k+T|t_k)\right)}\end{array}$$

(29)

subject to: system dynamics (30), (31) and control input admissible constraints (32).

$${\tilde{\mathit{x}}}_i\left(t_k|t_k\right)={\mathit{x}}_i\left(t_k\right)$$

(30)

$${\stackrel{·}{\stackrel{~}{\mathit{x}}}}_i\left(\tau |t_k\right)=f({\tilde{\mathit{x}}}_i\left(\tau |t_k\right),{\tilde{\mathit{u}}}_i\left(\tau |t_k\right))$$

(31)

$${\tilde{\mathit{u}}}_i\left(\tau |t_k\right)\in {\mathbb{U}}_i$$

(32)

In which $t_k=t_0+k\delta$ , $k=0,1,\cdots$ , $\delta$ denotes the sampling period of DMPC controller, $T$ is the control horizon satisfying $N_{\mathrm{c}}=T/\delta$ and $N_{\mathrm{c}}\in {\mathbb{Z}}^{+}$ . $L_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-},{\tilde{\mathit{u}}}_i)$ represents the stage cost, which is formulated as (33) according to definition 1, with ${\lambda }_{1i}$ , ${\lambda }_{2i}$ and ${\lambda }_{3i}$ being the preset positive weighting scalars

$$\begin{array}{c}\multicolumn{1}{c}{L_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-},{\tilde{\mathit{u}}}_i)=\underset{\mathrm{consensus}\mathrm{on}\mathrm{time}\mathrm{-}\mathrm{to}\mathrm{-}\mathrm{go}}{\underbrace{\stackrel{\mathrm{coupling}\mathrm{cost}}{\overbrace{{\lambda }_{1i}\left(\sum _{j\in {\mathbb{N}}_i}{({\tilde{x}}_{i2}\left(\tau |t_k\right)-{\tilde{x}}_{j2}\left(\tau |t_k\right))}^2\right)}}}}}\\ \multicolumn{1}{c}{+\stackrel{\mathrm{decoupling}\mathrm{costs}}{\overbrace{\underset{\mathrm{control}\mathrm{of}\mathrm{LOS}\mathrm{angles}}{\underbrace{{\lambda }_{2i}({\tilde{x}}_{i3}^2\left(\tau |t_k\right)+{\tilde{x}}_{i5}^2\left(\tau |t_k\right))}}+\underset{\mathrm{control}\mathrm{of}\mathrm{LOS}\mathrm{angular}\mathrm{rates}}{\underbrace{{\lambda }_{3i}({\tilde{x}}_{i4}^2\left(\tau |t_k\right)+{\tilde{x}}_{i6}^2\left(\tau |t_k\right))}}}}}\end{array}$$

(33)

$g_i\left({\tilde{\mathit{x}}}_i(t_k+T|t_k)\right)$ denotes the terminal cost (to be designed in subsection 3.2.3), which is a differentiable function satisfying $g_i\left(0\right)=0$ and $g_i\left({\tilde{\mathit{x}}}_i\left(t\right)\right)>0$ for any ${\tilde{\mathit{x}}}_i\left(t\right)\ne 0$ . ${\tilde{\mathit{x}}}_{i-}$ represents the collection of the states of neighbors, and ${\mathbb{N}}_i$ is the neighbor index set of ${\mathrm{M}}_i$ .

Assumed state trajectories and compatibility constraints

It is observed in (33), the stage cost function is divided into the coupling cost and decoupling costs. For the decoupling costs, the nominal state trajectory can be obtained according to the system dynamics (31) of each interceptor. However, for the coupling cost, it is not possible to exactly know the actual state trajectory of neighbors until they find their optimal control inputs. Therefore, in order to guarantee the synchronization in computing the control actions, an estimation for ${\tilde{\mathit{x}}}_{i-}\left(\tau |t_k\right)$ (named as assumed state trajectory) is required when solving the optimization problem. To avoid confusion, the definition and notation with regard to the assumed and actual state trajectories are given as follows.

• ${\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)$ : the optimal control input of ${\mathrm{M}}_i$ over time interval $\tau \in [t_k,t_k+T]$ , the first portion of which is implemented as the real-time guidance command.

• ${\tilde{\mathit{u}}}_i^a\left(\tau |t_k\right)$ : the assumed control input of ${\mathrm{M}}_i$ over time interval $\tau \in [t_k,t_k+T]$ , which is constructed as below.

$${\tilde{\mathit{u}}}_i^a\left(\tau |t_k\right)=\{\begin{array}{c}{\tilde{\mathit{u}}}_i^{*}\left(\tau |t_{k-1}\right),\tau \in [t_k,t_{k-1}+T)\hfill \\ \multicolumn{1}{c}{{\tilde{\mathit{u}}}_i^{\kappa }\left({\tilde{\mathit{x}}}_i^{*}\left(\tau |t_{k-1}\right)\right),\tau \in [t_{k-1}+T,t_k+T)}\end{array}$$

(34)

In equation (34), ${\tilde{\mathit{u}}}_i^{\kappa }\left({\tilde{\mathit{x}}}_i^{*}\left(\tau |t_{k-1}\right)\right)$ is generated by the terminal auxiliary controller. The terminal auxiliary controller and the other two terminal ingredients are all essential for ensuring the convergence of the system, and will be designed in subsection 3.2.3.

By implementing ${\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)$ and ${\tilde{\mathit{u}}}_i^a\left(\tau |t_k\right)$ , respectively, two different state trajectories will be obtained, they are

• ${\tilde{\mathit{x}}}_i^{*}\left(\tau |t_k\right)$ : optimal state trajectory, that is, the state trajectory obtained by implementing ${\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)$ over time interval $\tau \in [t_k,t_k+T]$ .

• ${\tilde{\mathit{x}}}_i^a\left(\tau |t_k\right)$ : assumed state trajectory, that is, the state trajectory obtained by implementing ${\tilde{\mathit{u}}}_i^a\left(\tau |t_k\right)$ over time interval $\tau \in [t_k,t_k+T]$ .

Replacing the actual state trajectory with the assumed one, the stage cost is then rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{L_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)={\lambda }_{1i}\left(\sum _{j\in {\mathbb{N}}_i}{({\tilde{x}}_{i2}\left(\tau |t_k\right)-{\tilde{x}}_{j2}^a\left(\tau |t_k\right))}^2\right)}\\ \multicolumn{1}{c}{+{\lambda }_{2i}({\tilde{x}}_{i3}^2\left(\tau |t_k\right)+{\tilde{x}}_{i5}^2\left(\tau |t_k\right))+{\lambda }_{3i}({\tilde{x}}_{i4}^2\left(\tau |t_k\right)+{\tilde{x}}_{i6}^2\left(\tau |t_k\right))}\end{array}$$

(35)

The introduction of the assumed state trajectory guarantees the synchronization of interceptors in computing their control actions, so long as each interceptor transmits its assumed trajectory to the neighbors before they solve their respective optimization problems. However, the uncertain deviation from the actual state trajectory will also be induced, thus a proper constraint should be designed to limit such a deviation and further to guarantee the convergence of the overall system.

According to the stage cost given by (35), the only state variable included in the coupling cost is ${\tilde{x}}_{i2}$ , so the deviation between the actual and assumed state trajectory is defined as

$${\epsilon }_i\left(\tau |t_k\right)={\tilde{x}}_{i2}^a\left(\tau |t_k\right)-{\tilde{x}}_{i2}\left(\tau |t_k\right)$$

(36)

The constraint to limit that deviation, named compatibility constraint, is designed as

$$C_i^{\mathrm{CB}}\left(\tau |t_k\right)=\left|{\epsilon }_i\left(\tau |t_k\right)\right|-{\upsilon }_i\left(t_k\right)\le 0$$

(37)

In (37),

$${\upsilon }_i\left(t_k\right)=\frac{-b_i\left(t_k\right)+\sqrt{b_i^2\left(t_k\right)-4a_i·c_i\left(t_k\right)}}{2a_i}$$

(38)

$$\{\begin{array}{c}{a}_i={\lambda }_{1i}\left|N_i\right|[T+({\gamma }_i-1)\delta ]+2(T-\delta )\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}\hfill \\ \multicolumn{1}{c}{b_i\left(t_k\right)=2(T-\delta )\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}{\phi }_{\mathrm{ji}}\left(t_k\right)}\\ \multicolumn{1}{c}{c_i\left(t_k\right)=-\delta {\gamma }_i{\lambda }_{1i}\sum _{j\in {\mathbb{N}}_i}{\zeta }_{\mathrm{ij}}^2\left(t_k\right)}\end{array}$$

(39)

In (39), $\left|N_i\right|$ is the number of the neighbors of ${\mathrm{M}}_i$ , ${\phi }_{\mathrm{ij}}\left(t_k\right)={\phi }_{\mathrm{ji}}\left(t_k\right)=\underset{\tau \in [t_k,t_k+T)}{max}|{\tilde{x}}_{i2}^a\left(\tau |t_k\right)-{\tilde{x}}_{j2}^a\left(\tau |t_k\right)|$ , ${\zeta }_{\mathrm{ij}}\left(t_k\right)={\zeta }_{\mathrm{ji}}\left(t_k\right)=\underset{\tau \in [t_k,t_k+T)}{min}|{\tilde{x}}_{i2}^a\left(\tau |t_k\right)-{\tilde{x}}_{j2}^a\left(\tau |t_k\right)|$ , and ${\gamma }_i\in (0,1)$ is the scalar given a priori, which affects the convergence speed of the multi-interceptor system.

The compatibility constraints (37) and the coupling cost ${\lambda }_{1i}\left(\sum _{j\in {\mathbb{N}}_i}{({\tilde{x}}_{i2}\left(\tau |t_k\right)-{\tilde{x}}_{j2}^a\left(\tau |t_k\right))}^2\right)$ are both to render the state of ${\mathrm{M}}_i$ close to that of its neighbors, thereby achieving the consensus on the remaining time in a fully distributed manner. In the next subsection, the terminal ingredients corresponding to the decoupling costs and coupling cost will be respectively developed to provide the guarantee for the convergence of LOS angles and angular rates, and to maintain the consensus among the interceptors.

Terminal ingredients

The terminal ingredients (including the terminal auxiliary controller, the terminal cost, and the terminal set constraints) are of great significance in ensuring the recursive feasibility and local convergence of each subsystem, as well as the stability of the overall multi-agent system. However, their development remains challenging for a cooperative guidance problem,^24,25,27 which prompts us to carry out the following design.

1. The terminal ingredients corresponding to the coupling cost. In conventional DMPC problems, the multi-agent systems are required to track a globally shared dynamic reference^21–23 this reference is vital for the development of terminal ingredients. However, for the situation where such a global reference cannot be timely obtained by the interceptors, the more general terminal ingredients are needed. For this consideration, we develop the terminal ingredients by designing a locally shared dynamic reference, this practice equally maintains the consensus of the overall system.

2. The terminal ingredients corresponding to the decoupling costs. The control of LOS angles and their rates can be deemed as the control problem for a nonlinear second-order (or double integrator) system, the DMPC problem of which has been addressed by Gao et al. ²³ However, the terminal region designed by Gao ²³ is given by a terminal equality constraint, which is possibly strict and may lead to an unsolvable optimization problem. Unlike the above work, we develop a less strict terminal set constraint and its corresponding terminal auxiliary controller, by employing which the positive invariance and the recursive feasibility can be guaranteed.

Prior to the development of terminal ingredients, we introduce the following definition, which specifies the requirements that must be satisfied by the terminal ingredients.^21,23

• Definition 2: The terminal set (or terminal region) ${\mathbf{\Omega }}_i\left(t_k\right)$ and the terminal auxiliary controller ${\tilde{\mathit{u}}}_i^{\kappa }(•)$ are such that: for any ${\tilde{\mathit{x}}}_i\left(t|t_k\right)\in {\mathbf{\Omega }}_i\left(t_k\right)$ and $\tau \ge t$ , by implementing ${\tilde{\mathit{u}}}_i\left(\tau \right)={\tilde{\mathit{u}}}_i^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)$ ,

$${\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\in {\mathbf{\Omega }}_i\left(t_k\right)$$

(40)

$${\tilde{\mathit{u}}}_i^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)\in {\mathbb{U}}_i$$

(41)

$$\sum _{i=1}^{N_{\mathrm{M}}}({\stackrel{·}{g}}_i\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)+L_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-},{\tilde{\mathit{u}}}_i))\le 0$$

(42)

In definition 2, equation (40) represents positive invariance, equation (41) represents control input admissibility, and equation (42) represents local convergence.

Based on definition 2, the terminal ingredients corresponding to the coupling cost are designed as follows.

The terminal cost is

$$g_{1i}\left({\tilde{x}}_{i2}\left(\tau |t_k\right)\right)={\rho }_{1i}{\tilde{x}}_{i2}^2\left(\tau |t_k\right)$$

(43)

with the positive scalar ${\rho }_{1i}$ to be given.

The terminal set constraint is given by

$${\tilde{x}}_{i2}\left(\tau |t_k\right)-x_{\mathrm{ir}}\left(\tau |t_k\right)\ge {\Lambda }_{i1}$$

(44)

with the locally dynamic reference $x_{\mathrm{ir}}\left(\tau |t_k\right)$ to be designed and the negative scalar ${\Lambda }_{i1}$ to be given.

The terminal auxiliary controller is designed as

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{u}}_{i1}^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)={\tilde{x}}_{i1}\left(\tau |t_k\right){\tilde{x}}_{i4}^2\left(\tau |t_k\right)+{\tilde{x}}_{i1}\left(\tau |t_k\right){\tilde{x}}_{i6}^2\left(\tau |t_k\right)}\\ \multicolumn{1}{c}{{\cos }^2({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\\ \multicolumn{1}{c}{+k_{1i}{\tilde{x}}_{i1}\left(\tau |t_k\right)({\tilde{x}}_{i2}\left(\tau |t_k\right)-x_{\mathrm{ir}}\left(\tau |t_k\right))/\left({\tilde{x}}_{i2}^2\left(\tau |t_k\right)\right)}\end{array}$$

(45)

with the control gain $k_{1i}$ to be given.

In (44), the locally dynamic reference is designed as

$$\begin{array}{c}{x}_{ir}\left(\tau |t_k\right)=x_{i2}^d\left(t_k\right)+t_k-\tau \\ x_{i2}^d\left(t_k\right)=\min \left\{\left(\left(\sum _{j\in {\mathbb{N}}_i}{x}_{j2}\left(t_k\right)\right)/\left|N_i\right|\right),x_{i2}^d\left(t_{k-1}\right)-\delta \right\}\end{array}$$

(46)

which satisfies $\mathrm{d}{x}_{\mathrm{ir}}\left(\tau |t_k\right)/\mathrm{d}\tau =-1$ .

Let ${\overline{\lambda }}_{1i}={\lambda }_{1i}\left|N_i\right|+\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}$ , then the scalars ${\rho }_{1i}$ , ${\Lambda }_{i1}$ , and $k_{1i}$ should meet

$${\rho }_{1i}>{\overline{\lambda }}_{1i}\left(\underset{\tau ,t_k}{max}{x}_{\mathrm{ir}}\left(\tau |t_k\right)\right)/2$$

(47)

$$k_{1i}>{\overline{\lambda }}_{1i}/\left(2{\rho }_{1i}\right)$$

(48)

$$(-2{\rho }_{1i}+{\overline{\lambda }}_{1i}\underset{\tau ,t_k}{max}{x}_{\mathrm{ir}}\left(\tau |t_k\right))/(2k_{1i}{\rho }_{1i}-{\overline{\lambda }}_{1i})\le {\Lambda }_{i1}<0$$

(49)

Thus, the design of terminal ingredients corresponding to the coupling cost is completed.

Similarly, according to the requirements in definition 2, the terminal ingredients corresponding to the decoupling costs are developed as follows.

The terminal cost is defined as

$$\begin{array}{c}\multicolumn{1}{c}{g_{2i}\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)={\rho }_{2i}{({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\mu }_i{\tilde{x}}_{i4}\left(\tau |t_k\right))}^2}\\ \multicolumn{1}{c}{+{\rho }_{2i}{({\tilde{x}}_{i5}\left(\tau |t_k\right)+{\mu }_i{\tilde{x}}_{i6}\left(\tau |t_k\right))}^2}\\ \multicolumn{1}{c}{+{\rho }_{3i}{\tilde{x}}_{i4}^2\left(\tau |t_k\right)+{\rho }_{3i}{\tilde{x}}_{i6}^2\left(\tau |t_k\right)}\end{array}$$

(50)

with positive scalars ${\rho }_{2i}$ , ${\rho }_{3i}$ to be given, and positive scalar ${\mu }_i$ given a priori.

The terminal set constraints are given by

$${\rho }_{2i}{({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\mu }_i{\tilde{x}}_{i4}\left(\tau |t_k\right))}^2+{\rho }_{3i}{\tilde{x}}_{i4}^2\left(\tau |t_k\right)\le {\Lambda }_{i2}$$

(51)

$${\rho }_{2i}{({\tilde{x}}_{i5}\left(\tau |t_k\right)+{\mu }_i{\tilde{x}}_{i6}\left(\tau |t_k\right))}^2+{\rho }_{3i}{\tilde{x}}_{i6}^2\left(\tau |t_k\right)\le {\Lambda }_{i3}$$

(52)

Since the local convergence of LOS angles and angular rates is guaranteed by applying the terminal cost and the terminal auxiliary controller (see the proof of theorem 1 in subsection 4.1), there are no extra requirements for the selection of ${\Lambda }_{i2}$ , ${\Lambda }_{i3}$ , which will not be specifically discussed in the follows. In applications, considering that the control input admissible set may be small, the positive scalars ${\Lambda }_{i2}$ , ${\Lambda }_{i3}$ can be set slightly larger to get a easily reaching but conservative terminal region, or can be optimized offline via methods proposed by Chen et al. ³⁴ and Yu et al. ³⁵ to obtain a less conservative terminal region.

• Remark 4: Different from the dual-mode MPC strategy, ³⁶ that is, if the states of ${\mathrm{M}}_i$ are inside the terminal region, the terminal auxiliary controller will be applied as the actual control input, our strategy is the same as Dai ²² and Gao, ²³ that is, the control action will be obtained by solving the optimization problem, regardless of whether the states are inside the terminal region or not.

The terminal auxiliary controller is designed as

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{u}}_{i2}^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)={\tilde{x}}_{i1}\left(\tau |t_k\right)(k_{2i}{\tilde{x}}_{i3}\left(\tau |t_k\right)+(({\rho }_{2i}-k_{2i}({\rho }_{3i}+{\rho }_{2i}{\mu }_i^2))/\left({\rho }_{2i}{\mu }_i\right)+2/{\tilde{x}}_{i2}\left(\tau |t_k\right)){\tilde{x}}_{i4}\left(\tau |t_k\right))}\\ \multicolumn{1}{c}{-{\tilde{x}}_{i1}\left(\tau |t_k\right){\tilde{x}}_{i6}^2\left(\tau |t_k\right)\sin ({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})\cos ({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\end{array}$$

(53)

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{u}}_{i3}^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)=-{\tilde{x}}_{i1}\left(\tau |t_k\right)\cos ({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})(({\rho }_{2i}-k_{2i}({\rho }_{3i}+{\rho }_{2i}{\mu }_i^2))/\left({\rho }_{2i}{\mu }_i\right)+2/{\tilde{x}}_{i2}\left(\tau |t_k\right)){\tilde{x}}_{i6}\left(\tau |t_k\right)}\\ \multicolumn{1}{c}{-k_{2i}{\tilde{x}}_{i1}\left(\tau |t_k\right){\tilde{x}}_{i5}\left(\tau |t_k\right)\cos ({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\\ \multicolumn{1}{c}{-2{\tilde{x}}_{i1}\left(\tau |t_k\right){\tilde{x}}_{i4}\left(\tau |t_k\right){\tilde{x}}_{i6}\left(\tau |t_k\right)\sin ({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\end{array}$$

(54)

with the control gain $k_{2i}$ to be given.

The scalars ${\rho }_{2i}$ , ${\rho }_{3i}$ in (50), and $k_{2i}$ in (53), are supposed to satisfy

$${\rho }_{3i}>{\mu }_i({\lambda }_{3i}+{\lambda }_{2i}{\mu }_i^2)/2$$

(55)

$$\begin{array}{c}\multicolumn{1}{c}{{\rho }_{2i}>({\lambda }_{2i}{\rho }_{3i}{\mu }_i^2+\sqrt{{\lambda }_{2i}{\mu }_i{\rho }_{3i}^2(2{\rho }_{3i}-{\lambda }_{3i}{\mu }_i)})/}\\ \multicolumn{1}{c}{({\mu }_i(2{\rho }_{3i}-{\lambda }_{3i}{\mu }_i)-{\lambda }_{2i}{\mu }_i^4)}\end{array}$$

(56)

$${\lambda }_{2i}/\left(2{\mu }_i{\rho }_{2i}\right)<k_{2i}<{\rho }_{2i}(2{\rho }_{3i}-{\lambda }_{3i}{\mu }_i)/\left(2{({\rho }_{2i}{\mu }_i^2+{\rho }_{3i})}^2\right)$$

(57)

where equations (55) and (56) are both to ensure the existence of $k_{2i}$ which satisfies equation (57).

Thus, the design of terminal ingredients corresponding to the decoupling costs is completed.

To sum up, the terminal ingredients for a single interceptor are given as follows.

$$g_i\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)=g_{1i}\left({\tilde{x}}_{i2}\left(\tau |t_k\right)\right)+g_{2i}\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)$$

(58)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\Omega }}_i\left(t_k\right)=}\\ \multicolumn{1}{c}{\left\{{\tilde{\mathit{x}}}_i\left(\tau |t_k\right)|\begin{array}{c}{\tilde{x}}_{i2}\left(\tau |t_k\right)-x_{\mathrm{ir}}\left(\tau |t_k\right)\ge {\Lambda }_{i1}\hfill \\ \multicolumn{1}{c}{{\rho }_{2i}{({\tilde{x}}_{i3}\left(\tau |t_k\right)+{\mu }_i{\tilde{x}}_{i4}\left(\tau |t_k\right))}^2+{\rho }_{3i}{\tilde{x}}_{i4}^2\left(\tau |t_k\right)\le {\Lambda }_{i2}}\\ \multicolumn{1}{c}{{\rho }_{2i}{({\tilde{x}}_{i5}\left(\tau |t_k\right)+{\mu }_i{\tilde{x}}_{i6}\left(\tau |t_k\right))}^2+{\rho }_{3i}{\tilde{x}}_{i6}^2\left(\tau |t_k\right)\le {\Lambda }_{i3}}\end{array}\right\}}\end{array}$$

(59)

$${\tilde{\mathit{u}}}_i^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)={[{\tilde{u}}_{i1}^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right),{\tilde{u}}_{i2}^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right),{\tilde{u}}_{i3}^{\kappa }\left({\tilde{\mathit{x}}}_i\left(\tau |t_k\right)\right)]}^{\mathrm{T}}$$

(60)

For ensuring the existence of the terminal ingredients that satisfy all the requirements in definition 2, the following assumption is desired.

• Assumption 4: given weighting scalars ${\lambda }_{1i}$ , ${\lambda }_{2i}$ , ${\lambda }_{3i}$ (in stage cost (35)) and ${\mu }_i$ (in terminal cost (50)), there exist positive scalars ${\rho }_{1i}$ , ${\rho }_{2i}$ , ${\rho }_{3i}$ satisfying (47), (55) and (56), control gains $k_{1i}$ , $k_{2i}$ satisfying (48) and (57), scalar ${\Lambda }_{i1}$ satisfying (49), and there always exists the terminal auxiliary controller (60) that satisfies (41) for any ${\tilde{\mathit{x}}}_i\left(\tau |t_k\right)$ .

Based upon the above designed terminal ingredients and compatibility constraints, the FHOCP problem for a single interceptor to achieve cooperative interception is defined as

• Problem 2: At the time instant $t_k$ , given the system state ${\mathit{x}}_i\left(t_k\right)$ , the assumed state trajectories of neighbors ${\tilde{\mathit{x}}}_{i-}^a\left(\tau |t_k\right)$ over the time interval $\tau \in [t_k,t_k+T]$ , find the optimal control inputs

$${\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)=\underset{{\tilde{\mathit{u}}}_i\left(\tau |t_k\right)}{\arg min}{J}_i(t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)$$

(61)

subject to: system dynamics (30), (31), control input admissible constraints (32), the compatibility constraint (37) and the terminal set constraint (62).

$${\tilde{\mathit{x}}}_i(t_k+T|t_k)\in {\mathbf{\Omega }}_i\left(t_k\right)$$

(62)

where $J_i(t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)={\int }_{t_k}^{t_k+T}{L}_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)\mathrm{d}\tau +g_i\left({\tilde{\mathit{x}}}_i(t_k+T|t_k)\right)$ , the stage cost $L_i(\tau |t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)$ , terminal cost $g_i\left({\tilde{\mathit{x}}}_i(t_k+T|t_k)\right)$ , and the terminal region ${\mathbf{\Omega }}_i\left(t_k\right)$ are shown in (35), (58), and (59), respectively.

Now, we have completed the design of the distributed cooperative guidance scheme consisting of an extended state observer and a DMPC controller, whose implementation, together with the optimization procedure of problem 2, will be elaborated in the next subsection.

Implementation of the distributed cooperative guidance scheme

Solving procedure based on zebra optimization algorithm

Solving an optimization problem described as problem 2 commonly requires large amounts of integral operations, thus, for improving the computational efficiency, we adopt the Euler’s method to approximate the system trajectory within the control horizon.

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathit{x}}}_i\left(t_{k+l+1}|t_k\right)\approx {\tilde{\mathit{x}}}_i\left(t_{k+l}|t_k\right)+\delta f({\tilde{\mathit{x}}}_i\left(t_{k+l}|t_k\right),{\tilde{\mathit{u}}}_i\left(t_{k+l}|t_k\right)),}\\ \multicolumn{1}{c}{l=0,1,\cdots ,N_{\mathrm{c}}-1}\end{array}$$

(63)

Applying (63), one can obtain the system states and control inputs in the form of sequences, that is,

$${\mathit{X}}_i\left(t_k\right)=[{\tilde{\mathit{x}}}_i\left(t_{k+1}|t_k\right),{\tilde{\mathit{x}}}_i\left(t_{k+2}|t_k\right),\cdots ,{\tilde{\mathit{x}}}_i(t_k+T|t_k)]$$

(64)

$${\mathit{U}}_i\left(t_k\right)=[{\tilde{\mathit{u}}}_i\left(t_k|t_k\right),{\tilde{\mathit{u}}}_i\left(t_{k+1}|t_k\right),\cdots ,{\tilde{\mathit{u}}}_i(t_{k-1}+T|t_k)]$$

(65)

The cost function, is accordingly approximated by

$$\begin{array}{c}\multicolumn{1}{c}{J_i(t_k,{\mathit{X}}_i,{\mathit{X}}_{i-}^a,{\mathit{U}}_i)=\delta \sum _{l=1}^{N_{\mathrm{c}}}{L}_i(t_{k+l}|t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)}\\ \multicolumn{1}{c}{+g_i\left({\tilde{\mathit{x}}}_i(t_k+T|t_k)\right)}\end{array}$$

(66)

All the constraints related to the states are supposed to be approximated in a same manner, such as the compatibility constraint, which will be approximated by (67) and can be further transformed into the constraint directly imposed on the control input ${\tilde{u}}_{i1}\left(t_{k+l}|t_k\right)$ by applying (67) and (68).

$$\left|{\tilde{x}}_{i2}\left(t_{k+l+1}|t_k\right)-{\tilde{x}}_{i2}\left(t_{k+l}|t_k\right)-\delta f_2\left({\tilde{x}}_i\left(t_{k+l}|t_k\right),{\tilde{u}}_{i1}\left(t_{k+l}|t_k\right)\right)\right|-{\upsilon }_i\left(t_k\right)\le 0$$

(67)

$$\begin{array}{c}\multicolumn{1}{c}{f_2({\tilde{\mathit{x}}}_i\left(t_{k+l}|t_k\right),{\tilde{u}}_{i1}\left(t_{k+l}|t_k\right))=-1+{\tilde{x}}_{i2}^2\left(t_{k+l}|t_k\right){\tilde{x}}_{i4}^2\left(t_{k+l}|t_k\right)}\\ \multicolumn{1}{c}{-{\tilde{x}}_{i2}^2\left(t_{k+l}|t_k\right){\tilde{u}}_{i1}\left(t_{k+l}|t_k\right)/{\tilde{x}}_{i1}\left(t_{k+l}|t_k\right)}\\ \multicolumn{1}{c}{+{\tilde{x}}_{i2}^2\left(t_{k+l}|t_k\right){\tilde{x}}_{i6}^2\left(t_{k+l}|t_k\right){\cos }^2({\tilde{x}}_{i3}\left(t_{k+l}|t_k\right)+{\theta }_{\mathrm{L}{\mathrm{d}}_i})}\end{array}$$

(68)

Similar to the compatibility constraint, the terminal set constraints can also be imposed on the control inputs ${\tilde{\mathit{u}}}_i(t_{k-1}+T|t_k)$ through approximation and transformation. For ease of denotation, we use additional admissible sets to represent the transformed constraints.

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathbb{U}}}_i^{\left(\beta \right)}\left(t_{k+l}|t_k\right)=}\\ \multicolumn{1}{c}{\{{\tilde{\mathit{u}}}_i\left(t_{k+l}|t_k\right)|{\tilde{u}}_{\mathrm{ij}}\left(t_{k+l}|t_k\right)\in [{\underset{¯}{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right),{\overline{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)],}\end{array}j=1,2,3\},\beta =1,2,3$$

(69)

where ${\tilde{\mathbb{U}}}_i^{\left(1\right)}\left(t_{k+l}|t_k\right)={\mathbb{U}}_i$ is the original control input admissible set, ${\tilde{\mathbb{U}}}_i^{\left(2\right)}\left(t_{k+l}|t_k\right)$ represents the admissible set for the compatibility constraint, in which only ${\tilde{u}}_{i1}\left(t_{k+l}|t_k\right)$ is constrained, the upper and lower bounds ${\underset{¯}{u}}_{i1}^{\left(2\right)}\left(t_{k+l}|t_k\right)$ , ${\overline{u}}_{i1}^{\left(2\right)}\left(t_{k+l}|t_k\right)$ can be calculated with (67) and (68). ${\tilde{\mathbb{U}}}_i^{\left(3\right)}\left(t_{k+l}|t_k\right)$ represents the admissible set for the terminal set constraints, with ${\underset{¯}{u}}_{\mathrm{ij}}^{\left(3\right)}(t_k+T|t_k)$ , ${\overline{u}}_{\mathrm{ij}}^{\left(3\right)}(t_k+T|t_k)$ , $j=1,2,3$ being calculated in the same way as (67) and (68), by transforming constraints (44), (51), and (52).

In order for implementing the transformed constraints, we define a saturation function given by (70).

$$\begin{array}{c}\multicolumn{1}{c}{\sigma \left({\tilde{u}}_{\mathrm{ij}}\left(t_{k+l}|t_k\right)|{\overline{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right),{\underset{¯}{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)\right)=}\\ \multicolumn{1}{c}{\{\begin{array}{c}{\overline{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right),{\tilde{u}}_{\mathrm{ij}}\left(t_{k+l}|t_k\right)>{\overline{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)\hfill \\ \multicolumn{1}{c}{{\tilde{u}}_{\mathrm{ij}}\left(t_{k+l}|t_k\right),{\underset{¯}{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)\le {\tilde{u}}_{\mathrm{ij}}\left(t_{k+l}|t_k\right)\le {\overline{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)}\\ \multicolumn{1}{c}{{\underset{¯}{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right),\mathrm{else}}\end{array}}\end{array}$$

(70)

Based on the above, an emerging swarm intelligence-based optimization algorithm, ZOA ³⁷ is utilized to solve the processed optimization problem 2. As a new bio-inspired metaheuristic algorithm, ZOA simulates the foraging behavior of zebras and their defense strategy against predators’ attacks, whose performance in solving high-dimensional problem has been validated through comparison with other well-known algorithms.

In ZOA, each zebra corresponds to a candidate solution of the optimization problem, which can be evaluated by an objective function. For the processed optimization problem 2, the position of each zebra is denoted as ${\mathit{Z}}_j\in {\mathbb{R}}_{3\times N_{\mathrm{c}}}$ , $j=1,2,\cdots ,N_{\mathrm{Z}}$ , with the same dimension as the control input sequence ${\mathit{U}}_i\left(t_k\right)$ , $N_{\mathrm{Z}}$ is the number of zebras in population. The objective function value of ${\mathit{Z}}_j$ , according to (66), is given by $F_j=J_i(t_k,{\mathit{X}}_i,{\mathit{X}}_{i-}^a,{\mathit{Z}}_j)$ , and the purpose of running such an algorithm is to find the zebra with the minimum objective function value, by iteratively updating the position of zebras.

The position update of ZOA is divided into two phases in each iteration, imitating the natural behaviors of zebras, they are

• (1) Foraging. In foraging phase, the zebra with the minimum objective function value is considered as the pioneer zebra, who will lead the other zebras toward its position. The position update in this phase is described as

$${\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_1)}={\mathit{Z}}_j+{\mathit{R}}_1\odot ({\mathit{Z}}^{\left(P\right)}-{\varsigma }_1{\mathit{Z}}_j)$$

(71)

$${\mathit{Z}}_j=\{\begin{array}{c}{\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_1)},F_j^{(\mathrm{new},{\mathrm{P}}_1)}<F_j\hfill \\ \multicolumn{1}{c}{{\mathit{Z}}_j,\mathrm{else}}\end{array}$$

(72)

where ${\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_1)}$ is the updated position, with the objective function value being $F_j^{(\mathrm{new},{\mathrm{P}}_1)}$ , ${\mathit{Z}}^{\left(P\right)}$ is the pioneer zebra, ${\mathit{R}}_1\in {\mathbb{R}}_{3\times N_{\mathrm{c}}}$ is a matrix composed of random numbers in interval $[0,1]$ , ${\varsigma }_1$ is a random integer satisfying $\varsigma \in \{1,2\}$ , and ⊙ is the operator of Hadamard product.

• (2) Defense against predators. In this phase, zebras adopt different strategies to resist attacks from different predators. For bigger predators like lions, zebras will escape in random sideways turning movements, for smaller predators such as hyenas, zebras will gather toward the attacked zebra to intimidate predators. The position update in this phase can be modeled as

$${\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_2)}=\{\begin{array}{c}{\mathit{Z}}_j+0.01(1-\eta /K_1){\mathit{R}}_2\odot {\mathit{Z}}_j,P_{\mathrm{s}}\le 0.5\hfill \\ \multicolumn{1}{c}{{\mathit{Z}}_j+{\mathit{R}}_3\odot ({\mathit{Z}}^{\left(\mathrm{A}\right)}-{\varsigma }_2{\mathit{Z}}_j)\mathrm{else}}\end{array}$$

(73)

$${\mathit{Z}}_j=\{\begin{array}{c}{\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_2)},F_j^{(\mathrm{new},{\mathrm{P}}_2)}<F_j\hfill \\ \multicolumn{1}{c}{{\mathit{Z}}_j,\mathrm{else}}\end{array}$$

(74)

where ${\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_2)}$ is the updated position, with the objective function value being $F_j^{(\mathrm{new},{\mathrm{P}}_2)}$ , ${\mathit{Z}}^{\left(A\right)}$ is the attacked zebra randomly selected, ${\mathit{R}}_2\in {\mathbb{R}}_{3\times N_{\mathrm{c}}}$ is a matrix composed of random numbers in interval $[-1,1]$ , ${\mathit{R}}_3$ and ${\varsigma }_2$ are the matrix and the constant with the same requirements as ${\mathit{R}}_1$ and ${\varsigma }_1$ , $\eta$ and $K_1$ respectively represent the current iteration count and the maximum number of iterations, with the latter being set in advance, $P_{\mathrm{s}}$ is a random number in interval $[0,1]$ determining which kind of defense strategy will be adopted.

All the zebras will go through check and correction after update lest they overstep the search area, that is, violate the admissible set constraints given by (69). The procedure of check and correction can be described using the following pseudo-code.

Algorithm 1: Pseudo-code of the position correction for each zebra
1: for $l:=$ 0 to $N_{\mathrm{c}}-1$ do
2:   Calculate ${\underset{¯}{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)$ , for $\beta =1,2,3$
3:   Correct ${\mathit{z}}_{j,l+1}$ using ${\mathit{z}}_{j,l+1}=\sigma \left({\mathit{z}}_{j,l+1}|{\overline{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right),{\underset{¯}{u}}_{\mathrm{ij}}^{\left(\beta \right)}\left(t_{k+l}|t_k\right)\right)$ for $\beta =1,2,3$
4:   Calculate ${\mathit{x}}_i\left(t_{k+l+1}|t_k\right)$ applying (63), with ${\tilde{u}}_{\mathrm{ij}}\left(t_{k+l}|t_k\right)$ being ${\mathit{z}}_{j,l+1}$
5: end for
6: Return corrected ${\mathit{Z}}_j$ and ${\mathit{X}}_i\left(t_k\right)$

And the complete procedure of obtaining ${\mathit{U}}_i^{*}\left(t_k\right)=[{\tilde{\mathit{u}}}_i^{*}\left(t_k|t_k\right),{\tilde{\mathit{u}}}_i^{*}\left(t_{k+1}|t_k\right),\cdots ,{\tilde{\mathit{u}}}_i^{*}(t_{k-1}+T|t_k)]$ is presented as follows.

Algorithm 2: Pseudo-code of ZOA
1: // Initialization
2: Set the parameters of the algorithm: the number of ebras in population $N_{\mathrm{Z}}$ , the maximum number of iterations $K_1$ , and the maximum number of invalid iterations $K_2$ .
3: Randomly generate the initial positions of zebras ${\mathit{Z}}_j$ , $j=1,2,\cdots ,N_{\mathrm{Z}}$ , in the neighborhood of the assumed control input sequence, i.e., $[{\tilde{\mathit{u}}}_i^{*}\left(t_k|t_{k-1}\right),{\tilde{\mathit{u}}}_i^{*}\left(t_{k+1}|t_{k-1}\right),\cdots ,{\tilde{\mathit{u}}}_i^{\kappa }\left({\tilde{\mathit{x}}}_i^{*}(t_{k-1}+T|t_{k-1})\right)]$ .
4: Initialize the current iteration count as $\eta =0$ , and the number of invalid iterations as $\eta \prime =0$ .
5: // Main loop
6: while $\eta <K_1\wedge \eta \prime <K_2$ do
7:  for $j:=1$ to $N_{\mathrm{Z}}$ do
8:   Correct ${\mathit{Z}}_j$ and obtain ${\mathit{X}}_i\left(t_k\right)$ using algorithm 1, calculate $F_j=J_i(t_k,{\mathit{X}}_i,{\mathit{X}}_{i-}^a,{\mathit{Z}}_j)$
9:  end for
10:  Determine the pioneer zebra ${\mathit{Z}}^{\left(P\right)}$ , and randomly select the attacked zebra in this round ${\mathit{Z}}^{\left(A\right)}$ , record the objective function value of ${\mathit{Z}}^{\left(P\right)}$
11:  if the objective function value of ${\mathit{Z}}^{\left(P\right)}$ remains unchanged compared to the last iteration, then
12:   $\eta \prime =\eta \prime +1$
13:  else
14:   $\eta \prime =0$
15:  for $j:=1$ to $N_{\mathrm{Z}}$ do
16:   In foraging phase, calculate and correct ${\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_1)}$ using (71) and algorithm 1, determine the updated ${\mathit{Z}}_j$ with (72).
17:   In the phase of defense against predators, calculate and correct ${\mathit{Z}}_j^{(\mathrm{new},{\mathrm{P}}_2)}$ using (73) and algorithm 1, determine the updated ${\mathit{Z}}_j$ with (74)
18:  end for
19:  Update the current iteration count as $\eta =\eta +1$
20: end while
21: Determine the optimal control input sequence at $t_k$ as   ${\mathit{U}}_i^{*}\left(t_k\right)=[{\tilde{\mathit{u}}}_i^{*}\left(t_k|t_k\right),{\tilde{\mathit{u}}}_i^{*}\left(t_{k+1}|t_k\right),\cdots ,{\tilde{\mathit{u}}}_i^{*}(t_{k-1}+T|t_k)]={\mathit{Z}}^{\left(P\right)}$
22: // Results
23: return ${\mathit{U}}_i^{*}\left(t_k\right)$

Complete implementation of the distributed cooperative guidance scheme

After obtaining ${\mathit{U}}_i^{*}\left(t_k\right)$ , whose first portion ${\tilde{\mathit{u}}}_i^{*}\left(t_k|t_k\right)$ will be implemented as the (nominal) control action, with the remaining portions used to construct the assumed control input sequence. The complete implementation of the proposed distributed cooperative guidance scheme is given as figure 3.

Figure 3.

Flowchart of the proposed DMPC cooperative guidance scheme.

In the next section, an analysis of the distributed cooperative guidance scheme will be conducted to demonstrate that the closed-loop stability of the multi-interceptor system can be ensured by implementing the proposed cooperative guidance scheme, with all the requirements of cooperative interception specified in definition 1 being achieved.

Feasibility and stability analysis of the distributed cooperative guidance scheme

Since the distributed cooperative guidance scheme is designed following the common practice of DMPC, two properties will be illustrated in this section, they are the recursive feasibility of the optimization problem 2 and the closed-loop stability for the overall system. With these two properties, the proposed cooperative guidance scheme is capable of achieving the cooperative interception as required in definition 1.

Recursive feasibility

In this subsection, we will prove that the optimization problem 2 is solvable for each interceptor at each sampling instant $t_k$ , which ensures the recursive feasibility of the proposed cooperative guidance scheme. See Theorem 1 for more details.

• Theorem 1: Let the terminal region and terminal cost be given as in (58) and (59). Suppose: (1) Assumption 4 holds, (2) ${\tilde{\mathit{x}}}_i\left(t|t_k\right)\in {\mathbf{\Omega }}_i\left(t_k\right)$ , (3) problem 2 is feasible at $t_k$ . Then, for any $\tau \ge t$ and $t_k>t_{\alpha }$ , under the terminal auxiliary controller (60),

• (a) all the requirements in definition 2 are satisfied.

• (b) the optimization problem 2 is feasible at every sampling instant after $t_k$ .

whose proof follows.

• Proof: (For the sake of brevity, some “ $\left(\tau \right)$ ” or “ $\left(\tau |t_k\right)$ ” are omitted in the following derivations.)

According to lemma 1, after the convergence of FPESO, the deviation between the actual state trajectory and the nominal state trajectory can be neglected, thus the system will evolve along the following dynamics when $t_k>t_{\alpha }$ .

$${\stackrel{·}{\mathit{x}}}_i=f({\mathit{x}}_i,{\tilde{\mathit{u}}}_i)$$

(75)

For this reason, “ ${\tilde{\mathit{x}}}_i\left(\tau |t_k\right)$ ” and “ ${\tilde{\mathit{x}}}_i\left(t_k\right)$ ” will be rewritten as “ ${\mathit{x}}_i\left(\tau |t_k\right)$ ”and “ ${\mathit{x}}_i\left(t_k\right)$ ” in the follows.

• (a) According to assumption 4, there always exists the terminal auxiliary controller ${\tilde{\mathit{u}}}_i^{\kappa }\left({\mathit{x}}_i\left(\tau |t_k\right)\right)$ that satisfies (41) for any ${\mathit{x}}_i\left(\tau |t_k\right)\in {\mathbf{\Omega }}_i\left(t_k\right)$ , that is, the control input admissibility is guaranteed.

Then, by implementing ${\tilde{u}}_{i1}^{\kappa }\left({\mathit{x}}_i\left(\tau |t_k\right)\right)$ , ${\tilde{u}}_{i2}^{\kappa }\left({\mathit{x}}_i\left(\tau |t_k\right)\right)$ , and ${\tilde{u}}_{i3}^{\kappa }\left({\mathit{x}}_i\left(\tau |t_k\right)\right)$ when $\tau \ge t$ , there are (76)– (78) according to (13), (21), and (75).

$${\stackrel{·}{x}}_{i2}=-k_{1i}(x_{i2}-x_{\mathrm{ir}})-1$$

(76)

$${\stackrel{·}{x}}_{i4}=-k_{2i}{x}_{i3}-({\rho }_{2i}-k_{2i}({\rho }_{3i}+{\rho }_{2i}{\mu }_i^2))x_{i4}/\left({\rho }_{2i}{\mu }_i\right)$$

(77)

$${\stackrel{·}{x}}_{i6}=-k_{2i}{x}_{i5}-({\rho }_{2i}-k_{2i}({\rho }_{3i}+{\rho }_{2i}{\mu }_i^2))x_{i6}/\left({\rho }_{2i}{\mu }_i\right)$$

(78)

Applying (76) and $\mathrm{d}{x}_{\mathrm{ir}}\left(\tau |t_k\right)/\mathrm{d}\tau =-1$ , yields

$$\mathrm{d}(x_{i2}\left(\tau |t_k\right)-x_{\mathrm{ir}}\left(\tau |t_k\right))/\mathrm{d}\tau =-k_{1i}(x_{i2}-x_{\mathrm{ir}})$$

(79)

According to (79), there is $\mathrm{d}(x_{i2}\left(\tau |t_k\right)-x_{\mathrm{ir}}\left(\tau |t_k\right))/\mathrm{d}\tau >0$ for any $x_{i2}\left(t|t_k\right)$ that satisfies $x_{i2}\left(t|t_k\right)-x_{\mathrm{ir}}\left(t|t_k\right)\ge {\Lambda }_{i1}$ and $x_{i2}\left(t|t_k\right)-x_{\mathrm{ir}}\left(t|t_k\right)<0$ , which guarantees the satisfaction of (44) for $\tau =t^{+}$ .

Applying (77) and (57), the constraint (51) is satisfied for $\tau =t^{+}$ , because

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{d}({\rho }_{2i}{(x_{i3}\left(\tau |t_k\right)+{\mu }_i{x}_{i4}\left(\tau |t_k\right))}^2+{\rho }_{3i}{x}_{i4}^2\left(\tau |t_k\right))/\mathrm{d}\tau }\\ \multicolumn{1}{c}{=2{\rho }_{2i}(x_{i3}+{\mu }_i{x}_{i4})(x_{i4}+{\mu }_i{\stackrel{·}{x}}_{i4})+2{\rho }_{3i}{x}_{i4}{\stackrel{·}{x}}_{i4}}\\ \multicolumn{1}{c}{=-2k_{2i}{\mu }_i{\rho }_{2i}{x}_{i3}^2+2(k_{2i}{({\rho }_{3i}+{\rho }_{2i}{\mu }_i^2)}^2-{\rho }_{2i}{\rho }_{3i})x_{i4}^2/\left({\rho }_{2i}{\mu }_i\right)}\\ \multicolumn{1}{c}{<-2k_{2i}{\mu }_i{\rho }_{2i}{x}_{i3}^2-{\lambda }_{3i}{x}_{i4}^2}\\ \multicolumn{1}{c}{<0}\end{array}$$

(80)

Similarly, applying (78) and (57), the constraint (52) is satisfied for $\tau =t^{+}$ , thus, by induction, the condition (40) holds for any $\tau \ge t$ , that is, the positive invariance is guaranteed.

Then, by dividing the stage cost into the coupling cost and the decoupling costs, one can get

$$L_i(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)=L_{1i}(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)+L_{2i}(\tau |t_k,{\mathit{x}}_i,{\tilde{\mathit{u}}}_i)$$

(81)

where

$$L_{1i}(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)={\lambda }_{1i}\left(\sum _{j\in {\mathbb{N}}_i}{(x_{i2}\left(\tau |t_k\right)-x_{j2}\left(\tau |t_k\right))}^2\right)$$

(82)

$$\begin{array}{c}\multicolumn{1}{c}{L_{2i}(\tau |t_k,{\mathit{x}}_i,{\tilde{\mathit{u}}}_i)={\lambda }_{2i}(x_{i3}^2\left(\tau |t_k\right)+x_{i5}^2\left(\tau |t_k\right))}\\ \multicolumn{1}{c}{+{\lambda }_{3i}(x_{i4}^2\left(\tau |t_k\right)+x_{i6}^2\left(\tau |t_k\right))}\end{array}$$

(83)

Applying inequality ${(x_{i2}-x_{j2})}^2<x_{i2}^2+x_{j2}^2$ for $x_{i2},x_{j2}>0$ , yields

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{i=1}^{N_{\mathrm{M}}}{L}_{1i}(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)=\sum _{i=1}^{N_{\mathrm{M}}}{\lambda }_{1i}\sum _{j\in {\mathbb{N}}_i}{(x_{i2}-x_{j2})}^2}\\ \multicolumn{1}{c}{<\sum _{i=1}^{N_{\mathrm{M}}}{\lambda }_{1i}\sum _{j\in {\mathbb{N}}_i}(x_{i2}^2+x_{j2}^2)=\sum _{i=1}^{N_{\mathrm{M}}}({\lambda }_{1i}\left|N_i\right|+\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j})x_{i2}^2}\end{array}$$

(84)

Therefore, the inequality $L_{1i}(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)+{\stackrel{·}{g}}_{1i}\left(x_{i2}\left(\tau |t_k\right)\right)\le 0$ holds if

$$({\lambda }_{1i}\left|N_i\right|+\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j})x_{i2}^2+2{\rho }_{1i}{x}_{i2}{\stackrel{·}{x}}_{i2}\le 0$$

(85)

Substituting (44), (47)– (49), (76) and ${\overline{\lambda }}_{1i}={\lambda }_{1i}\left|N_i\right|+\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}$ into (85), one can obtain

$$\begin{array}{c}\multicolumn{1}{c}{({\lambda }_{1i}\left|N_i\right|+\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j})x_{i2}^2+2{\rho }_{1i}{x}_{i2}{\stackrel{·}{x}}_{i2}}\\ \multicolumn{1}{c}{=({\overline{\lambda }}_{1i}-2k_{1i}{\rho }_{1i})x_{i2}^2+2{\rho }_{1i}(k_{1i}{x}_{\mathrm{ir}}-1)x_{i2}}\\ \multicolumn{1}{c}{=({\overline{\lambda }}_{1i}-2k_{1i}{\rho }_{1i})(x_{i2}-x_{\mathrm{ir}})x_{i2}+({\overline{\lambda }}_{1i}{x}_{\mathrm{ir}}-2{\rho }_{1i})x_{i2}}\\ \multicolumn{1}{c}{\le (({\overline{\lambda }}_{1i}-2k_{1i}{\rho }_{1i}){\Lambda }_{i1}+({\overline{\lambda }}_{1i}\underset{\tau ,t_k}{max}{x}_{\mathrm{ir}}-2{\rho }_{1i}))x_{i2}}\\ \multicolumn{1}{c}{\le 0}\end{array}$$

(86)

Accordingly, inequality $L_{1i}(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)+{\stackrel{·}{g}}_{1i}\left(x_{i2}\left(\tau |t_k\right)\right)\le 0$ holds.

Applying (57) and (80), there is

$$\begin{array}{c}\multicolumn{1}{c}{L_{2i}(\tau |t_k,{\mathit{x}}_i,{\tilde{\mathit{u}}}_i)+{\stackrel{·}{g}}_{2i}\left({\mathit{x}}_i\left(\tau |t_k\right)\right)}\\ \multicolumn{1}{c}{<{\lambda }_{2i}(x_{i3}^2+x_{i5}^2)+{\lambda }_{3i}(x_{i4}^2+x_{i6}^2)-2k_{2i}{\mu }_i{\rho }_{2i}(x_{i3}^2+x_{i5}^2)}\\ \multicolumn{1}{c}{-{\lambda }_{3i}(x_{i4}^2+x_{i6}^2)}\\ \multicolumn{1}{c}{<0}\end{array}$$

(87)

According to $L_{1i}(\tau |t_k,{\mathit{x}}_i,{\mathit{x}}_{i-},{\tilde{\mathit{u}}}_i)+{\stackrel{·}{g}}_{1i}\left(x_{i2}\left(\tau |t_k\right)\right)\le 0$ and (87), condition (42) holds for any $\tau \ge t$ , that is, the local convergence is guaranteed.

In summary, conclusion (a) holds.

• (b) Since optimization problem 2 is feasible at $t_k$ , one can get a candidate solution at $t_{k+1}$ , that is, ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ . Following the property of the control input admissible set, and according to the assumption that there always exists the terminal auxiliary controller ${\tilde{\mathit{u}}}_i^{\kappa }\left({\mathit{x}}_i\left(\tau |t_k\right)\right)$ satisfying (41), the candidate solution ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ meet ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)\in {\mathbb{U}}_i$ .

By implementing ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ over time interval $\tau \in [t_{k+1},t_k+T]$ , the system states at $t_k+T$ , that is, ${\mathit{x}}_i(t_k+T|t_{k+1})$ satisfy ${\mathit{x}}_i(t_k+T|t_{k+1})\in {\mathbf{\Omega }}_i\left(t_k\right)$ , which implies

$$x_{i2}(t_k+T|t_{k+1})-x_{\mathrm{ir}}(t_k+T|t_k)\ge {\Lambda }_{i1}$$

(88)

$$\begin{array}{c}\multicolumn{1}{c}{{\rho }_{2i}{(x_{i3}(t_k+T|t_{k+1})+{\mu }_i{x}_{i4}(t_k+T|t_{k+1}))}^2}\\ \multicolumn{1}{c}{+{\rho }_{3i}{x}_{i4}^2(t_k+T|t_{k+1})\le {\Lambda }_{i2}}\end{array}$$

(89)

$$\begin{array}{c}\multicolumn{1}{c}{{\rho }_{2i}{(x_{i5}(t_k+T|t_{k+1})+{\mu }_i{x}_{i6}(t_k+T|t_{k+1}))}^2}\\ \multicolumn{1}{c}{+{\rho }_{3i}{x}_{i6}^2(t_k+T|t_{k+1})\le {\Lambda }_{i3}}\end{array}$$

(90)

According to (46), there is

$$\begin{array}{c}\multicolumn{1}{c}{x_{\mathrm{ir}}(t_k+T|t_{k+1})=x_{i2}^d\left(t_{k+1}\right)+t_{k+1}-t_k-T=x_{i2}^d\left(t_{k+1}\right)}\\ \multicolumn{1}{c}{+\delta -T\le x_{i2}^d\left(t_k\right)-T=x_{\mathrm{ir}}(t_k+T|t_k)}\end{array}$$

(91)

Applying (91) to (88), yields

$$\begin{array}{c}\multicolumn{1}{c}{x_{i2}(t_k+T|t_{k+1})-x_{\mathrm{ir}}(t_k+T|t_{k+1})\ge x_{i2}(t_k+T|t_{k+1})}\\ \multicolumn{1}{c}{-x_{\mathrm{ir}}(t_k+T|t_k)\ge {\Lambda }_{i1}}\end{array}$$

(92)

Combining (92) with (89), (90), one can obtain ${\mathit{x}}_i(t_k+T|t_{k+1})\in {\mathbf{\Omega }}_i\left(t_{k+1}\right)$

Then, based on condition (40) in definition 2 (which has been proven to hold in conclusion (a)), by implementing ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ over time interval $\tau \in [t_k+T,t_{k+1}+T]$ , there is ${\mathit{x}}_i(t_{k+1}+T|t_{k+1})\in {\mathbf{\Omega }}_i\left(t_{k+1}\right)$ , that is, the terminal set constraint (62) is satisfied for ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ .

Since the state trajectory obtained by implementing ${\tilde{\mathit{u}}}_i^a\left(\tau |t_k\right)$ over time interval $\tau \in [t_k,t_k+T]$ is the same as the assumed state trajectory, the compatibility constraint (37) is naturally satisfied for ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ .

In summary, for the candidate solution ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ , all the constraints of problem 2 at $t_{k+1}$ are satisfied, thus ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ is a feasible solution to the optimization problem 2 at $t_{k+1}$ . By induction, the optimization problem 2 is feasible at every sampling instant after $t_k$ , and conclusion (b) holds.

Thereupon, the proof of theorem 1 is completed. ▪

Closed-loop stability

In this subsection, the multi-interceptor system is proven to be stable under the proposed cooperative guidance scheme, that is, all the interceptors are guaranteed to achieve consensus on impact time, and intercept the target at their desired terminal LOS angles. The details are as follows.

• Theorem 2: For the multi-interceptor system, each interceptor solves the optimization problem 2 and obtains the optimal control input ${\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)$ at each sampling instant $t_k$ . Then, by implementing the real-time guidance command ${\mathit{u}}_i\left(\tau \right)={\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)+{\hat{\mathit{w}}}_i\left(\tau \right)$ , $\tau \in [t_k,t_{k+1}]$ , $t_k>t_{\alpha }$ , all the conditions in definition 1 will be satisfied, that is, the multi-interceptor system is guaranteed to achieve cooperative interception.

Whose proof follows.

• Proof: (Similar to the proof of theorem 1, some “ $\left(\tau \right)$ ” or “ $\left(\tau |t_k\right)$ ” are omitted, and “ ${\tilde{\mathit{x}}}_i\left(\tau |t_k\right)$ ”, “ ${\tilde{\mathit{x}}}_i\left(t_k\right)$ ” are replaced by “ ${\mathit{x}}_i\left(\tau |t_k\right)$ ”, “ ${\mathit{x}}_i\left(t_k\right)$ ” in the following derivations)

For the multi-interceptor system, the Lyapunov function is selected as

$$V\left(t_k\right)=\sum _{i=1}^{N_{\mathrm{M}}}{J}_i(t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^{a*},{\tilde{\mathit{u}}}_i^{*})$$

(93)

which satisfies $V\left(t_k\right)>0$ for any ${\mathit{x}}_i\left(\tau |t_k\right)\ne 0$ , according to (61), (35), and (58).

Since ${\tilde{\mathit{u}}}_i^a\left(\tau |t_{k+1}\right)$ is a feasible solution to the optimization problem 2 at $t_{k+1}$ , one can obtain

$$\begin{array}{c}\multicolumn{1}{c}{V\left(t_{k+1}\right)-V\left(t_k\right)}\\ \multicolumn{1}{c}{\le \sum _{i=1}^{N_{\mathrm{M}}}(J_i(t_{k+1},{\mathit{x}}_i^a,{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^a)-J_i(t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^{a*},{\tilde{\mathit{u}}}_i^{*}))}\\ \multicolumn{1}{c}{=\sum _{i=1}^{N_{\mathrm{M}}}\left(\begin{array}{c}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}({(x_{i2}^{*}\left(\tau |t_k\right)-x_{j2}^{*}\left(\tau |t_k\right))}^2-{(x_{i2}^{*}\left(\tau |t_k\right)-x_{j2}^a\left(\tau |t_k\right))}^2)\mathrm{d}\tau \hfill \\ \multicolumn{1}{c}{-{\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau +{\int }_{t_k+T}^{t_{k+1}+T}{L}_i(\tau |t_k,{\mathit{x}}_i^{\kappa },{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{\kappa })\mathrm{d}\tau }\\ \multicolumn{1}{c}{+g_i\left({\mathit{x}}_i^{\kappa }(t_{k+1}+T|t_k)\right)-g_i\left({\mathit{x}}_i^{*}(t_k+T|t_k)\right)}\end{array}\right)}\end{array}$$

(94)

In which ${\mathit{x}}_i^{\kappa }(•)$ denotes the system states obtained by implementing ${\tilde{\mathit{u}}}_i^{\kappa }(•)$

Then, applying condition (42) in definition 2 (as proven in theorem 1), there is (95) for any $\tau \ge t_k+T$ , by implementing ${\tilde{\mathit{u}}}_i^{\kappa }\left({\mathit{x}}_i^{*}\left(\tau |t_k\right)\right)$ .

$$\sum _{i=1}^{N_{\mathrm{M}}}({\stackrel{·}{g}}_i\left({\mathit{x}}_i^{*}\left(\tau |t_k\right)\right)+L_i(\tau |t_k,{\mathit{x}}_i^{\kappa },{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{\kappa }))<0$$

(95)

Integrating (95) from $t_k+T$ to $t_{k+1}+T$ , yields

$$\sum _{i=1}^{N_{\mathrm{M}}}({\int }_{t_k+T}^{t_{k+1}+T}{L}_i(\tau |t_k,{\mathit{x}}_i^{\kappa },{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{\kappa })\mathrm{d}\tau +g_i\left({\mathit{x}}_i^{\kappa }(t_{k+1}+T|t_k)\right)-g_i\left({\mathit{x}}_i^{*}(t_k+T|t_k)\right))<0$$

(96)

Let ${\epsilon }_i^{*}\left(\tau |t_k\right)=x_{i2}^a\left(\tau |t_k\right)-x_{i2}^{*}\left(\tau |t_k\right)$ . Employing the compatibility constraint (37), one can get

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}({(x_{i2}^{*}\left(\tau |t_k\right)-x_{j2}^{*}\left(\tau |t_k\right))}^2-{(x_{i2}^{*}\left(\tau |t_k\right)-x_{j2}^a\left(\tau |t_k\right))}^2)\mathrm{d}\tau }\\ \multicolumn{1}{c}{=\sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}(2x_{i2}^{*}-x_{j2}^{*}-x_{j2}^a){\epsilon }_j^{*}\mathrm{d}\tau =\sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}(2(x_{i2}^a-x_{j2}^a)-2{\epsilon }_i^{*}+{\epsilon }_j^{*}){\epsilon }_j^{*}\mathrm{d}\tau }\\ \multicolumn{1}{c}{\le \sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}(2{\phi }_{\mathrm{ij}}+2\left|{\epsilon }_i^{*}\right|+\left|{\epsilon }_j^{*}\right|)\left|{\epsilon }_j^{*}\right|\mathrm{d}\tau =\sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}(2{\phi }_{\mathrm{ij}}\left|{\epsilon }_j^{*}\right|+2\left|{\epsilon }_i^{*}\right|\left|{\epsilon }_j^{*}\right|+{\left|{\epsilon }_j^{*}\right|}^2)\mathrm{d}\tau }\\ \multicolumn{1}{c}{\le \sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}(2{\phi }_{\mathrm{ij}}\left|{\epsilon }_j^{*}\right|+{\left|{\epsilon }_i^{*}\right|}^2+2{\left|{\epsilon }_j^{*}\right|}^2)\mathrm{d}\tau =\sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}2{\lambda }_{1j}{\phi }_{\mathrm{ij}}\left|{\epsilon }_i^{*}\right|+{\lambda }_{1i}{\left|{\epsilon }_i^{*}\right|}^2+2{\lambda }_{1j}{\left|{\epsilon }_i^{*}\right|}^2\mathrm{d}\tau }\\ \multicolumn{1}{c}{\le \sum _{i=1}^{N_{\mathrm{M}}}\sum _{j\in {\mathbb{N}}_i}(T-\delta )(2{\lambda }_{1j}{\phi }_{\mathrm{ij}}{\nu }_i+({\lambda }_{1i}+2{\lambda }_{1j}){{\nu }_i}^2)}\\ \multicolumn{1}{c}{=\sum _{i=1}^{N_{\mathrm{M}}}(2(T-\delta )\left(\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}{\phi }_{\mathrm{ij}}\right){\nu }_i+(T-\delta )({\lambda }_{1i}\left|N_i\right|+2\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}){{\nu }_i}^2)}\end{array}$$

(97)

and

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{i=1}^{N_{\mathrm{M}}}{\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau }\\ \multicolumn{1}{c}{\ge \sum _{i=1}^{N_{\mathrm{M}}}{\int }_{t_k}^{t_{k+1}}\left(\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1i}{(x_{i2}^{*}-x_{j2}^a)}^2\right)\mathrm{d}\tau =\sum _{i=1}^{N_{\mathrm{M}}}{\int }_{t_k}^{t_{k+1}}\left(\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1i}{((x_{i2}^a-x_{j2}^a)-{\epsilon }_i^{*})}^2\right)\mathrm{d}\tau }\\ \multicolumn{1}{c}{\ge \sum _{i=1}^{N_{\mathrm{M}}}{\int }_{t_k}^{t_{k+1}}\left(\sum _{j\in {\mathbb{N}}_i}({\lambda }_{1i}{{\zeta }_{\mathrm{ij}}}^2-{\lambda }_{1i}{{\nu }_i}^2)\right)\mathrm{d}\tau =\sum _{i=1}^{N_{\mathrm{M}}}(\delta {\lambda }_{1i}\sum _{j\in {\mathbb{N}}_i}{{\zeta }_{\mathrm{ij}}}^2-\delta {\lambda }_{1i}\left|N_i\right|{{\nu }_i}^2)}\end{array}$$

(98)

Combining (98) with (97), and substituting (38) and (39), there is

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{i=1}^{N_{\mathrm{M}}}(\left(\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}({(x_{i2}^{*}-x_{j2}^{*})}^2-{(x_{i2}^{*}-x_{j2}^a)}^2)\mathrm{d}\tau \right)-{\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau )}\\ \multicolumn{1}{c}{=\sum _{i=1}^{N_{\mathrm{M}}}(\left(\sum _{j\in {\mathbb{N}}_i}{\int }_{t_{k+1}}^{t_k+T}{\lambda }_{1i}({(x_{i2}^{*}-x_{j2}^{*})}^2-{(x_{i2}^{*}-x_{j2}^a)}^2)\mathrm{d}\tau \right)-{\gamma }_i{\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau )}\\ \multicolumn{1}{c}{+\sum _{i=1}^{N_{\mathrm{M}}}\left(({\gamma }_i-1){\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau \right)}\\ \multicolumn{1}{c}{\le \sum _{i=1}^{N_{\mathrm{M}}}(2(T-\delta )\left(\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}{\phi }_{\mathrm{ij}}\right){\nu }_i+(T-\delta )({\lambda }_{1i}\left|N_i\right|+2\sum _{j\in {\mathbb{N}}_i}{\lambda }_{1j}){{\nu }_i}^2+\delta {\gamma }_i{\lambda }_{1i}\left|N_i\right|{{\nu }_i}^2-\delta {\gamma }_i{\lambda }_{1i}\sum _{j\in {\mathbb{N}}_i}{{\zeta }_{\mathrm{ij}}}^2)}\\ \multicolumn{1}{c}{+\sum _{i=1}^{N_{\mathrm{M}}}\left(({\gamma }_i-1){\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau \right)}\\ \multicolumn{1}{c}{=\sum _{i=1}^{N_{\mathrm{M}}}(b_i{\nu }_i+a_i{{\nu }_i}^2+c_i)+\sum _{i=1}^{N_{\mathrm{M}}}\left(({\gamma }_i-1){\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau \right)}\\ \multicolumn{1}{c}{=\sum _{i=1}^{N_{\mathrm{M}}}\left(({\gamma }_i-1){\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau \right)}\end{array}$$

(99)

According to (94), (96) and (99), the following inequality holds.

$$V\left(t_{k+1}\right)-V\left(t_k\right)<\sum _{i=1}^{N_{\mathrm{M}}}\left(({\gamma }_i-1){\int }_{t_k}^{t_{k+1}}{L}_i(\tau |t_k,{\mathit{x}}_i^{*},{\mathit{x}}_{i-}^a,{\tilde{\mathit{u}}}_i^{*})\mathrm{d}\tau \right)<0$$

(100)

From (100), one can conclude that, the Lyapunov function of the multi-interceptor system will exponentially decrease to zero, that is, all the requirements in definition 1 will be satisfied, and the multi-interceptor system is guaranteed to achieve cooperative interception.

The proof is completed.▪

Numerical simulations

In order to verify the effectiveness of the proposed DMPC-based cooperative guidance (DMPC-CG) scheme, numerical simulations are carried out in this section. In different simulation scenarios, a group of 4 interceptors with fixed communication topology are required to attack a maneuvering target which adopts diverse maneuvering strategies. The communication topology among the interceptors is shown in Figure 4, the initial positions and velocities of interceptors and target are given in Table 1, the acceleration magnitude of interceptors is up to 200 m·s⁻².

Figure 4.

Communication topology among interceptors.

Table 1.

Initial states of interceptors.

Index of interceptor	Initial position (m)	Initial velocity (m/s)	Initial pitch path angle (°)	Initial yaw path angle (°)
1	(−7000, −2000, 6000)	400	67.5	0
2	(−6000, 1000, 0)	400	45	−22.5
3	(0, 0, −6000)	400	33.75	−67.5
4	(6000, −3000, −6000)	400	56.25	−90
Target	(15,000, 15,000, 15,000)	-	−30	135

The parameters of the DMPC controller are listed as

$$\begin{array}{c}\multicolumn{1}{c}{{\lambda }_{1i}=1,{\lambda }_{2i}=1,{\lambda }_{3i}=1,{\mu }_i=1,{\rho }_{1i}=40{\overline{\lambda }}_{1i},}\\ \multicolumn{1}{c}{{\rho }_{2i}=4,{\rho }_{3i}=1.5,k_{1i}=2,k_{2i}=0.125,{\Lambda }_{i1}=-0.2515}\end{array}$$

The parameters of the extended state observer are shown as

$$\alpha =50,{\beta }_1=15,{\beta }_2=50,{\beta }_3=15,\vartheta =0.7,\varphi =1$$

The parameters of the Zebra optimization algorithm are given as

$$N_{\mathrm{Z}}=10,K_1=20,K_2=5$$

The sampling period and control horizon of the DMPC controller are set to $\delta =0.1\mathrm{s}$ , $T=0.3\mathrm{s}$ , the sampling period for the system states and ESO is set to 0.001 s. The simulation terminates when the relative velocity becomes positive, or when the time-to-go is <0.001 s.

• Remark 5: provided the time-to-go is <0.001 s, the last sampling for the system states will be conducted after time-to-go to obtain the terminal miss distance and the impact time (terminal time of guidance).

Comparative experiments without external noise

In this subsection, the interceptors are required to perform the interception tasks without external noises. In order to validate the effectiveness of the DPMC cooperative guidance scheme in challenging task settings, the experiments will be carried out in the following scenarios.

• Case 1: In case 1, the initial LOS angles of each interceptor are set to deviate significantly from their expected values, and the interceptors will engage with a target whose velocity is comparable to that of themselves. Associated settings are given as follows.

Desired terminal LOS angles are

$$\begin{array}{c}{\theta }_{\mathrm{L}{\mathrm{d}}_1}=11.25°,{\theta }_{\mathrm{L}{\mathrm{d}}_2}=33.75°{\theta }_{\mathrm{L}{\mathrm{d}}_3}=45°{\theta }_{\mathrm{L}{\mathrm{d}}_4}=22.5°\\ {\psi }_{\mathrm{L}{\mathrm{d}}_1}=-45°,{\psi }_{\mathrm{L}{\mathrm{d}}_2}=-56.25°,{\psi }_{\mathrm{L}{\mathrm{d}}_3}=-22.5°{\psi }_{\mathrm{L}{\mathrm{d}}_4}=-33.75°\end{array}$$

The velocity and acceleration magnitudes of the target in inertial reference frame are

$$v_{\mathrm{T}}\left(t_0\right)=400\mathrm{m}/\mathrm{s},\{\begin{array}{c}{a}_{\mathrm{T}x}\left(t\right)=0\mathrm{m}·{\mathrm{s}}^{-2}\hfill \\ \multicolumn{1}{c}{a_{\mathrm{T}y}\left(t\right)=50\cos (0.5t)\mathrm{m}·{\mathrm{s}}^{-2}}\\ \multicolumn{1}{c}{a_{\mathrm{T}z}\left(t\right)=-t\cos (0.4t)\mathrm{m}·{\mathrm{s}}^{-2}}\end{array}$$

• Case 2: In case 2, the interceptors are required to intercept a target whose velocity is faster than themselves. In order for evading the interception, the target will adopt more drastic maneuvers in the yaw plane. Associated settings are given as follows.

Desired terminal LOS angles are

$$\begin{array}{c}{\theta }_{\mathrm{L}{\mathrm{d}}_1}=45°,{\theta }_{\mathrm{L}{\mathrm{d}}_2}=11.25°{\theta }_{\mathrm{L}{\mathrm{d}}_3}=22.5°{\theta }_{\mathrm{L}{\mathrm{d}}_4}=56.25°\hfill \\ {\psi }_{\mathrm{L}{\mathrm{d}}_1}=0°,{\psi }_{\mathrm{L}{\mathrm{d}}_2}=-22.5°,{\psi }_{\mathrm{L}{\mathrm{d}}_3}=-45°,{\psi }_{\mathrm{L}{\mathrm{d}}_4}=-78.75°\hfill \end{array}$$

The velocity and acceleration magnitudes of the target in inertial reference frame are

$$\begin{array}{c}\multicolumn{1}{c}{v_{\mathrm{T}}\left(t_0\right)=700\mathrm{m}/\mathrm{s},}\\ \multicolumn{1}{c}{\{\begin{array}{c}{a}_{\mathrm{T}x}\left(t\right)=0\mathrm{m}·{\mathrm{s}}^{-2}\hfill \\ \multicolumn{1}{c}{a_{\mathrm{T}y}\left(t\right)=50\cos (0.5t)\mathrm{m}·{\mathrm{s}}^{-2}}\\ \multicolumn{1}{c}{a_{\mathrm{T}z}\left(t\right)=-100\tanh ([(t-4)/10]-2[(t-4)/20]-0.5)\mathrm{m}·{\mathrm{s}}^{-2}}\end{array}}\end{array}$$

where $\tanh (•)$ is the hyperbolic tangent function.

Besides, in order to demonstrate the advantages of the proposed DMPC-CG, the fixed-time cooperative guidance (FxTCG) law ¹⁹ will be used for comparison, and the simulation results for the two cases are shown in Figures 5 and 6 and Tables 2 and 3, respectively.

Figure 5.

Simulation results of case 1. (a) 3D flight trajectories of interceptors and the target (b) variation curves of LOS angles (c) variation curves of times-to-go (d) acceleration commands in LOS frame.

Figure 6.

Simulation results of case 2. (a) 3D flight trajectories of interceptors and the target (b) variation curves of LOS angles (c) variation curves of times-to-go (d) acceleration commands in LOS frame.

Table 2.

Simulation results of case 1.

Method	Index of interceptor	Miss distance (m)	Error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Error of LOS angle ${\psi }_{\mathrm{L}}$ (°)	Impact time (s)
FxTCG	1	0.2253	6.47 × 10−3	9.06 × 10−5	39.3680
	2	0.2266	5.69 × 10−3	1.40 × 10−3	39.3682
	3	0.2295	5.20 × 10−3	0.1210	39.3680
	4	0.2247	6.12 × 10−3	0.0556	39.3682
	average	0.2265	5.87 × 10−3	0.0445	39.3681
DMPC-CG	1	4.74 × 10−3	2.56 × 10−3	9.56 × 10−4	38.4360
	2	6.56 × 10−4	1.05 × 10−3	2.52 × 10−3	38.4360
	3	1.24 × 10−3	8.33 × 10−5	4.07 × 10−3	38.4360
	4	2.40 × 10−3	3.20 × 10−3	2.08 × 10−3	38.4360
	average	2.26 × 10−3	1.72 × 10−3	2.41 × 10−3	38.4360

Table 3.

Simulation results of case 2.

Method	Index of interceptor	Miss distance (m)	Error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Error of LOS angle ${\psi }_{\mathrm{L}}$ (°)	Impact time (s)
FxTCG	1	0.0186	0.0103	0.0356	28.3168
	2	0.0185	0.0118	0.0133	28.3165
	3	0.0191	0.0114	3.42 × 10−3	28.3168
	4	0.0185	0.0105	8.39 × 10−3	28.3165
	average	0.0187	0.0110	0.0152	28.3167
DMPC-CG	1	0.0106	4.16 × 10−3	5.22 × 10−3	27.8487
	2	8.92 × 10−3	3.88 × 10−3	1.11 × 10−3	27.8487
	3	2.07 × 10−3	2.29 × 10−3	5.09 × 10−3	27.8487
	4	0.0460	3.09 × 10−3	6.42 × 10−3	27.8487
	average	0.0169	3.35 × 10−3	4.46 × 10−3	27.8487

It can be observed from Figures 5(a) and 6(a) that all the interceptors eventually hit the target with smooth trajectories. However, in case 1, the interceptors using FxTCG perform more rapid maneuvers, resulting in the strong deflections of the flight trajectories. Similar deflections also occur on the variation curves of the LOS angles, as shown in Figure 5(b). By contrast, the interceptors adopting DMPC-CG have smoother variation curves of LOS angles in both cases owing to the reasonable design of the DMPC controller. From Figures 5(c) and 6(c), it can be seen that all the interceptors reach the consensus on times-to-go using either method, subtle difference lies in the convergence time. Adopting DMPC-CG, the convergence times of consensus errors are about 2 s and 2 s, shorter than those using FxTCG (3 s and 2.5 s). According to Figures 5(d) and 6(d), except for the small fluctuations induced by the computation results of the stochastic optimization algorithm, the acceleration commands under DMPC-CG do not exhibit intense changes during guidance. By comparison, under FxTCG, the lateral and normal acceleration commands ( $u_{i2}$ and $u_{i3}$ ) experience the continuous saturation in the middle stages of guidance, the tangential acceleration commands ( $u_{i1}$ ) exhibit chattering behaviors in the final stage, which has a negative impact on the final guidance effect. Tables 2 and 3 provide a detailed comparison of the guidance performance of the two methods. It can be seen that DMPC-CG outperforms FxTCG in terms of terminal miss distance, convergence errors of LOS angles and consensus errors in ideal cases without external noise.

The next subsection will compare the performance of the two methods in the presence of external noise and disturbances.

Comparative experiments in the presence of external noise

In this subsection, the interceptors will perform tasks in the presence of measurement noises and wind disturbances to test the robustness of the DMPC-CG scheme. Similar to subsection 5.1, the simulations will be conducted in the following scenarios.

• Case 3: In case 3, the measurement noises of the radar seeker are taken into account. As the signal-to-noise ratio of LOS angular rates are much smaller compared with the relative distance or LOS angles, ³⁸ we assume that the measured LOS angular rates contain the non-Gaussian noise.

Models of the bounded measurement noises can be described as

$$\{\begin{array}{c}{\widehat{\stackrel{·}{\theta }}}_{{\mathrm{L}}_i}={\stackrel{·}{\theta }}_{{\mathrm{L}}_i}+{\varpi }_{{\theta }_i}\hfill \\ \multicolumn{1}{c}{{\widehat{\stackrel{·}{\psi }}}_{{\mathrm{L}}_i}={\stackrel{·}{\psi }}_{{\mathrm{L}}_i}+{\varpi }_{{\psi }_i}}\end{array}$$

(101)

where ${\widehat{\stackrel{·}{\theta }}}_{{\mathrm{L}}_i}$ and ${\widehat{\stackrel{·}{\psi }}}_{{\mathrm{L}}_i}$ are the measured values of LOS angular rates that will be used by the controller and the ESO, ${\varpi }_{{\theta }_i}$ and ${\varpi }_{{\psi }_i}$ are the random measurement noises which satisfy $\left|{\varpi }_{{\theta }_i}\right|\le {\varpi }_{max}$ and $\left|{\varpi }_{{\psi }_i}\right|\le {\varpi }_{max}$ , and ${\varpi }_{max}$ is set to 0.05 (°)/s. Other settings are the same as case 1.

• Case 4: In case 4, the wind disturbances in the real environment are taken into consideration, which will randomly act on the interceptors during their flight. Models of the bounded wind disturbances can be described as ²⁴

$$\{\begin{array}{c}{\stackrel{·}{x}}_{{\mathrm{M}}_i}=v_{{\mathrm{M}}_i}\cos {\theta }_{{\mathrm{M}}_i}\cos {\psi }_{{\mathrm{M}}_i}+{\varpi }_{x_i}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{y}}_{{\mathrm{M}}_i}=v_{{\mathrm{M}}_i}\sin {\theta }_{{\mathrm{M}}_i}+{\varpi }_{y_i}}\\ \multicolumn{1}{c}{{\stackrel{·}{z}}_{{\mathrm{M}}_i}=-v_{{\mathrm{M}}_i}\cos {\theta }_{{\mathrm{M}}_i}\sin {\psi }_{{\mathrm{M}}_i}+{\varpi }_{z_i}}\end{array}$$

(102)

where ${\varpi }_{x_i}$ , ${\varpi }_{y_i}$ , and ${\varpi }_{z_i}$ are the random wind disturbances which satisfy $\left|{\varpi }_{x_i}\right|\le {\overline{\varpi }}_{max}$ , $\left|{\varpi }_{y_i}\right|\le {\overline{\varpi }}_{max}$ , and $\left|{\varpi }_{z_i}\right|\le {\overline{\varpi }}_{max}$ . ${\overline{\varpi }}_{max}$ is set to 5 m/s, other settings are the same as case 2.

Comparative experiments will be conducted to test the performance of the designed FPESO (DMPC-FPESO), fixed-time ESO, ¹⁹ and linear ESO (LESO) ³⁹ in terms of disturbance and noise tolerance. The fixed-time ESO is tested together with the FxTCG, and the LESO will be integrated into the DMPC (DMPC-LESO) controller for testing.

Simulation results can be seen in Figures 7 and 8 and Tables 4 and 5. The consensus error of ${\mathrm{M}}_i$ in Figures 7(c) and 8(c) is defined as $x_{i2}-(\sum _{j=1}^{N_{\mathrm{M}}}{x}_{j2}/N_{\mathrm{M}})$ . From Figures 7(a), (b), 8(a) and (b), it can be seen that the flight trajectories and the variation curves of LOS angles do not exhibit significant changes compared to the ideal cases in subsection 5.1. However, due to external noises and disturbances, strong fluctuations occur during the convergence of consensus errors, as shown in Figures 7(c) and 8(c). Despite the presence of noises, DMPC-FPESO and DMPC-LESO still achieve smaller consensus errors than the FxTCG, which benefits from the reasonable design of the DMPC controller. Figures 7(d), (e), 8(d) and (e) depict the response curves and the estimation errors of different ESO, and it can be observed that FPESO benefits better noise tolerance performance than the LESO and fixed-time ESO, which is further demonstrated by the mean square error (MSE) of disturbance estimation in Tables 4 and 5. In addition, compared to the case with wind disturbances, all of the three ESO has better performance in the presence of measurement noise, which is because the impact of measurement noises is getting smaller as the LOS angular rates converge. Therefore, in the final stage of guidance, the estimated value gradually approximates the actual disturbances. By contrast, due to the lack of consideration of wind disturbances in the kinetic model of guidance, the estimation error will always float within a certain range, rather than converging to the origin. (according to Figure 7(d) and (e)). Tables 4 and 5 provide the detailed comparison among the three methods, and it can be seen that the overall guidance performance of DMPC-FPESO is superior to FxTCG and DMPC-LESO, due to the better performance of FPESO in noise suppression.

Figure 7.

Simulation results of case 3. (a) 3D flight trajectories of interceptors and the target (b) variation curves of LOS angles (c) variation curves of consensus errors (d) Response of extended state observers (e) disturbance estimation errors.

Figure 8.

Simulation results of case 4. (a) 3D flight trajectories of interceptors and the target (b) variation curves of LOS angles (c) variation curves of consensus errors (d) Response of extended state observers (e) disturbance estimation errors.

Table 4.

Simulation results of case 3.

Method	Index of interceptor	Miss distance (m)	Error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Error of LOS angle ${\psi }_{\mathrm{L}}$ (°)	Impact time (s)	MSE of disturbance estimation
DMPC-LESO	1	0.0305	4.21 × 10−3	0.0154	38.4243	667.1361
	2	0.0257	0.0429	0.0182	38.4243	639.0058
	3	5.98 × 10−3	0.0291	9.19 × 10−3	38.4243	625.8210
	4	0.0360	0.0383	0.0653	38.4243	603.8300
	Average	0.0246	0.0286	0.0270	38.4243	633.9482
FxTCG	1	0.0952	0.0288	0.0165	39.3434	329.6029
	2	0.0958	8.13 × 10−3	0.0298	39.3438	298.7480
	3	0.0970	0.0196	0.1335	39.3434	285.1601
	4	0.0950	2.41 × 10−3	7.01 × 10−3	39.3438	302.4189
	Average	0.0957	0.0147	0.0467	39.3436	303.9825
DMPC-FPESO	1	0.0120	0.0286	4.28 × 10−3	38.4664	76.4243
	2	9.04 × 10−3	7.89 × 10−3	0.0194	38.4664	71.9725
	3	0.0445	8.58 × 10−3	2.37 × 10−3	38.4664	67.3472
	4	0.0233	0.0233	4.07 × 10−3	38.4664	66.9844
	Average	0.0222	0.0171	7.53 × 10−3	38.4664	70.6821

Table 5.

Simulation results of case 4.

Method	Index of interceptor	Miss distance (m)	Error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Error of LOS angle ${\psi }_{\mathrm{L}}$ (°)	Impact time (s)	MSE of disturbance estimation
DMPC-LESO	1	0.0333	4.17 × 10−4	0.0342	27.6356	2306.5931
	2	0.0492	0.0186	0.1025	27.6355	2506.5638
	3	0.1103	0.0178	0.0441	27.6357	2372.3005
	4	0.1353	4.16 × 10−3	0.2554	27.6355	2367.3695
	Average	0.0820	0.0102	0.1091	27.6356	2388.2067
FxTCG	1	0.0844	0.0470	0.0474	28.2804	769.7578
	2	0.0769	0.0426	0.0593	28.2803	803.0019
	3	0.0705	0.0320	0.0168	28.2804	797.7900
	4	0.0768	0.0853	0.0222	28.2803	782.8878
	Average	0.0772	0.0517	0.0364	28.2803	788.3594
DMPC-FPESO	1	0.0458	0.0308	0.0364	27.3942	668.4990
	2	0.1095	0.0238	0.0143	27.3941	677.0499
	3	0.0701	1.59 × 10−3	2.14 × 10−3	27.3942	652.1748
	4	0.0939	0.0135	0.0315	27.3943	685.8194
	Average	0.0798	0.0174	0.0211	27.3942	670.8858

In the next subsection, a comparative experiment will be conducted to demonstrate the advantages of DMPC-CG over other DMPC-based cooperative guidance schemes.

Performance comparison with constraint-free models

In order to demonstrate the advantages of the proposed DMPC-CG, the comparison is made with the synchronous DMPC models which ignored the design of compatibility constraints and terminal ingredients.^24,27 The constraint-free DMPC model can be described as

• Problem 3: At the time instant $t_k$ , given the system state $x_i\left(t_k\right)$ , the assumed state trajectories of neighbors ${\tilde{\mathit{x}}}_{i-}^a\left(\tau |t_k\right)$ over the time interval $\tau \in [t_k,t_k+T]$ , find the optimal control inputs

$${\tilde{\mathit{u}}}_i^{*}\left(\tau |t_k\right)=\underset{{\tilde{\mathit{u}}}_i\left(\tau |t_k\right)}{\arg min}{\Gamma }_i(t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)$$

(103)

subject to: system dynamics (30), (31) and control input admissible constraints (32).

In (103),

$$\begin{array}{c}\multicolumn{1}{c}{{\Gamma }_i(t_k,{\tilde{\mathit{x}}}_i,{\tilde{\mathit{x}}}_{i-}^a,{\tilde{\mathit{u}}}_i)={\int }_{t_k}^{t_k+T}{\lambda }_{1i}\left(\sum _{j\in {\mathbb{N}}_i}{({\tilde{x}}_{i2}\left(\tau |t_k\right)-{\tilde{x}}_{j2}^a\left(\tau |t_k\right))}^2\right)\mathrm{d}\tau }\\ \multicolumn{1}{c}{+{\rho }_{1i}{\tilde{x}}_{i1}^2(t_k+T|t_k)+{\rho }_{2i}(\left|{\tilde{x}}_{i3}(t_k+T|t_k)\right|+\left|{\tilde{x}}_{i5}(t_k+T|t_k)\right|)}\end{array}$$

(104)

The simulation will be carried out under the same settings as case 1. Simulation results are shown in Figure 9 and Table 6.

Figure 9.

Comparative experimental results of DMPC-CG and constraint-free DMPC. (a) variation curves of LOS angles (b) variation curves of times-to-go (c) variation curves of the objective function value of constraint-free DMPC (d) variation curves of the objective function value of DMPC-CG.

Table 6.

Comparative experimental results of DMPC-CG and constraint-free DMPC.

Method	Index of interceptor	Miss distance (m)	Error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Error of LOS angle ${\psi }_{\mathrm{L}}$ (°)	Impact time (s)
Constraint-free DMPC	1	0.0372	0.0121	0.0275	36.7411
	2	0.0244	0.0150	0.0371	36.7414
	3	0.0926	0.1743	0.0267	36.7417
	4	0.3451	0.0711	0.3656	36.7428
	Average	0.1248	0.0681	0.1142	36.7418
DMPC-CG	1	4.74 × 10−3	2.56 × 10−3	9.56 × 10−4	38.4360
	2	6.56 × 10−4	1.05 × 10−3	2.52 × 10−3	38.4360
	3	1.24 × 10−3	8.33 × 10−5	4.07 × 10−3	38.4360
	4	2.40 × 10−3	3.20 × 10−3	2.08 × 10−3	38.4360
	Average	2.26 × 10−3	1.72 × 10−3	2.41 × 10−3	38.4360

Figure 9(a) shows that, the variation curves of LOS angles under constraint-free DMPC exhibit continuous fluctuations throughout the guidance process, which is mainly due to the lack of terminal auxiliary controllers (53) and (54), and an inappropriate objective function design. From Figure 9(b), the times-to-go under constraint-free DMPC seem to reach a consensus in the initial stage of guidance, however, the consensus errors among interceptors have not been fully eliminated before the end of guidance. By comparison, benefiting from the compatibility constraint (37) and the terminal set constraint (44), the DMPC-CG achieves better consensus effect throughout the entire guidance process. Figure 9(c) and (d) respectively depict the variation curves of the objective function values under the two methods. Using either of them, the objective function value shows a decreasing trend for most of the time. However, without the guarantee provided by the compatibility constraint and the terminal ingredient, the closed-loop stability under constraint-free DMPC is broken in the final stage of guidance. The detailed performance comparison between the two methods is shown in Table 6, and the DMPC-CG has huge advantages over the constraint-free DMPC model, which further demonstrates the significance of the compatibility constraint and the terminal ingredients designed in this paper.

Test for different settings of parameters

Test for parameters of the DMPC controller

In this subsection, we will provide a tuning guideline of the parameters by testing the cooperative guidance effect under different parameter settings. As shown in the stage cost (35), the weighting scalars ${\lambda }_{1i}$ , ${\lambda }_{2i}$ , ${\lambda }_{3i}$ reflect the importance of each control objective (consensus on time-to-go, convergence of LOS angles and convergence of LOS angular rates). The larger the scalar is given, the more the corresponding objective is emphasized. Given that these objectives are considered equally important in conventional cooperative guidance problems, the impact of the weighting scalars ${\lambda }_{1i}$ , ${\lambda }_{2i}$ , ${\lambda }_{3i}$ on the guidance effect will not be further explored. The main focus is on the tune of ${\rho }_{1i}$ , ${\rho }_{2i}$ , ${\rho }_{3i}$ , and $k_{1i}$ , $k_{2i}$ , as they are mutually related and have great impact on the convergence performance.

Test for ${\rho }_{1i}$ and $k_{1i}$

${\rho }_{1i}$ and $k_{1i}$ mainly affect the consensus on times-to-go. The relationship between the value range of $k_{1i}$ and ${\rho }_{1i}$ is shown in (48), and it can be seen that as the ${\rho }_{1i}$ increases, the lower bound of $k_{1i}$ will decrease. Furthermore, according to (49), The terminal region is simultaneously determined by ${\rho }_{1i}$ and $k_{1i}$ . Provided that the guidance will last for no >40 s (i.e., $\underset{\tau ,t_k}{max}{x}_{\mathrm{ir}}\left(\tau |t_k\right)<40$ ), the terminal region related scalar ${\Lambda }_{i1}$ will vary with $k_{1i}$ and the ratio of ${\rho }_{1i}/{\overline{\lambda }}_{1i}$ as shown in Figure 10. It can be observed that the terminal region becomes larger as the ratio ${\rho }_{1i}/{\overline{\lambda }}_{1i}$ increases and the control gain $k_{1i}$ decreases.

Figure 10.

Variation of ${\Lambda }_{i1}$ with respect to ${\rho }_{1i}/{\overline{\lambda }}_{1i}$ and $k_{1i}$ .

For further testing the impact of different parameter selections on the convergence performance of consensus errors, we choose the following three sets of parameters for simulation.

• Parameters 1: ${\rho }_{1i}=40{\overline{\lambda }}_{1i}$ , $k_{1i}=2$ , ${\Lambda }_{i1}=-0.2515$

• Parameters 2: ${\rho }_{1i}=40{\overline{\lambda }}_{1i}$ , $k_{1i}=4$ , ${\Lambda }_{i1}=-0.1253$

• Parameters 3: ${\rho }_{1i}=30{\overline{\lambda }}_{1i}$ , $k_{1i}=2$ , ${\Lambda }_{i1}=-0.1680$

The simulation will be carried out in the scenario of case 2, with all other parameters being the same as those given before.

Simulation results can be seen in figure 11, and it shows that the consensus errors under all three sets of parameters are the same at the end of guidance, while different parameters exhibit diverse variation trends in the convergence process.

Figure 11.

Variation of consensus errors under different parameter settings of DMPC controller.

The variation under the second set of parameters indicates that, the larger the gain $k_{1i}$ is chosen, the faster the consensus errors converge. However, an excessively large gain can also lead to the overshoot. According to the variation under the third set of parameters, a smaller ${\rho }_{1i}$ results in a smaller terminal region, and causes fluctuations in the convergence process. Thus, we should integrate the choice of ${\rho }_{1i}$ and $k_{1i}$ to obtain a proper convergence speed and a smoother variation curve. From the experiment results, the first set of parameters is superior to the other two, and can be considered as a reference for the selection of ${\rho }_{1i}$ and $k_{1i}$ .

Test for ${\rho }_{2i}$ , ${\rho }_{3i}$ , and $k_{2i}$

${\rho }_{2i}$ , ${\rho }_{3i}$ , and $k_{2i}$ determine the convergence performance of LOS angles and angular rates. For the pre given ${\lambda }_{2i}=1$ , ${\lambda }_{3i}=1$ , ${\mu }_i=1$ , the lower bound of ${\rho }_{3i}$ can be calculated with (55). And the relationship between ${\rho }_{3i}$ and the value range of ${\rho }_{2i}$ is presented in Figure 12.

Figure 12.

Variation of the value range of ${\rho }_{2i}$ with respect to ${\rho }_{3i}$ .

Based on (57), the variation curves of the upper and lower bounds of $k_{2i}$ with respect to ${\rho }_{2i}$ and ${\rho }_{3i}$ is depicted in Figure 13. It can be observed that as ${\rho }_{2i}$ and ${\rho }_{3i}$ increase simultaneously, the value range of the gain $k_{2i}$ becomes larger.

Figure 13.

Variation of the upper and lower bounds of $k_{2i}$ with respect to ${\rho }_{2i}$ and ${\rho }_{3i}$ .

In order to further test the guidance effect under different parameter settings, we select the following three sets of parameters for simulation.

• Parameters 4: ${\rho }_{2i}=4$ , ${\rho }_{3i}=1.5$ , $k_{2i}=0.125$

• Parameters 5: ${\rho }_{2i}=3$ , ${\rho }_{3i}=2$ , $k_{2i}=0.18$

• Parameters 6: ${\rho }_{2i}=5$ , ${\rho }_{3i}=1.3$ , $k_{2i}=0.1$

The simulation will be conducted in the scenario of case 2, with all other parameters being the same as those given in subsection 5.1.

Simulation results can be seen in Figure 14 and Table 7. It shows that as the control gain $k_{2i}$ increases, the convergence speed of the LOS angles becomes faster, but the overshoot also occurs. Furthermore, a larger ${\rho }_{3i}$ typically leads to a smaller terminal miss distance, while an excessively large ${\rho }_{2i}$ and a small ${\rho }_{3i}$ may result in an inappropriate gain $k_{2i}$ and poor control performance of LOS angles. Therefore, selecting ${\rho }_{2i}$ and ${\rho }_{3i}$ with proper size and proportion has great impact on the control performance of LOS angles and angular rates. Based on the above results, the fourth set of parameters achieves better guidance performance compared to the other two, and can be provided as a reference for the selection of ${\rho }_{2i}$ , ${\rho }_{3i}$ , and $k_{2i}$ .

Figure 14.

Variation curves of LOS angles under different parameter settings of the DMPC controller.

Table 7.

Simulation results under different parameter settings of the DMPC controller.

Set of parameters	Average miss distance (m)	Average error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Average error of LOS angle ${\psi }_{\mathrm{L}}$ (°)
4	0.0169	3.35 × 10−3	4.46 × 10−3
5	2.02 × 10−3	6.17 × 10−3	0.0144
6	0.0101	0.0604	0.0328

Test for parameters of ZOA

In this subsection, we will analyze the computational efficiency under different parameter settings of the Zebra optimization algorithm (including the number of zebras in population $N_{\mathrm{Z}}$ and the maximum number of iterations $K_1$ ). Generally speaking, the larger $N_{\mathrm{Z}}$ and $K_1$ are chosen, the smaller objective function value will be obtained. However, as $N_{\mathrm{Z}}$ and $K_1$ exceed certain values, continuing to increase them will not lead to the significant reduction of the objective function value, but will increase the computation time. Therefore, choosing appropriate $N_{\mathrm{Z}}$ and $K_1$ is crucial for ensuring both computational efficiency and guidance performance.

The simulation will be carried out under the settings of case 2, with the hardware of Intel (R) Core (™) i7-13620H. The following three sets of parameters are used for comparative experiment.

• Parameters 7: $N_{\mathrm{Z}}=10$ , $K_1=20$

• Parameters 8: $N_{\mathrm{Z}}=20$ , $K_1=40$

• Parameters 9: $N_{\mathrm{Z}}=5$ , $K_1=15$

The other parameters are chosen as in subsection 5.1.

As shown in Figure 15 and Table 8, although larger $N_{\mathrm{Z}}$ and $K_1$ lead to a smaller objective function value and a better guidance performance, the time required for computation will also greatly increase, which is up to almost half of the sampling period. On the other hand, when $N_{\mathrm{Z}}$ and $K_1$ are both chosen to be small, the time complexity of computation will be greatly reduced, but at the cost of the reduction of the guidance performance. Therefore, the selection of $N_{\mathrm{Z}}$ and $K_1$ should achieve the balance between the computational efficiency and the guidance performance.

Figure 15.

Simulation results under different parameter settings of ZOA.

Table 8.

Simulation results under different parameter settings of ZOA.

Set of parameters	Average/max computation time (s)	Average miss distance (m)	Average error of LOS angle ${\theta }_{\mathrm{L}}$ (°)	Average error of LOS angle ${\psi }_{\mathrm{L}}$ (°)	Average consensus error (s)
7	4.55 × 10−3/0.0156	0.0169	3.35 × 10−3	4.46 × 10−3	<10−5
8	8.72 × 10−3/0.0391	4.17 × 10−3	5.43 × 10−3	6.41 × 10−3	<10−5
9	1.50 × 10−3/0.0117	8.73 × 10−3	7.55 × 10−3	0.0152	<10−5

As a general rule, when the computation time does not exceed the one-tenth of the sampling period of controller, the control delay induced by computation can be neglected, so the parameters 7 and 9 are preferred from this perspective.

Conclusion

This paper presents a DMPC-based 3D cooperative guidance scheme for multiple interceptors simultaneously attacking a maneuvering target. The proposed cooperative guidance scheme consists of a FPESO and a DMPC controller, which can achieve the convergence of LOS angles and the consensus on impact time among all interceptors. The details can be summarized as follows.

• (1) The FPESO is constructed and used to estimate the disturbances caused by target maneuver, which shows good performance in convergence speed and noise suppression.

• (2) The ZOA is employed to solve the FHOCP of each interceptor, the optimization procedure can be carried out in a completely distributed manner owing to the use of the assumed state trajectory, whose deviation from the actual state trajectory is limited by the designed compatibility constraint.

• (3) The terminal ingredients corresponding to the coupling and decoupling stage costs are developed, respectively. Applying them to the DMPC controller, the resulting FHOCP is feasible at each sampling instant, and the stability of the multi-interceptor system is also guaranteed.

• (4) The effectiveness and the superiority of the proposed DMPC-based cooperative guidance scheme are verified by the theoretical proofs and simulation results.

Future study will aim to extend the proposed methods to tackle the cooperative guidance problem for multiple interceptors in the presence of enemy radars or anti-interception weapons.