Cooperative Integrated Guidance and Control for Active Target Protection in Three-Player Conflict

Abstract

This paper addresses the active target protection problem in a three-player (Target–Attacker–Defender, TAD) conflict by proposing a cooperative integrated guidance and control (IGC) strategy. Unlike previous studies that have designed guidance and control loops separately, this work establishes an IGC model by linearizing both the translational motion and the rotational motion of the vehicles, thereby generating actuator commands directly. This model integrates the kinematics and short-period dynamics, providing a more comprehensive and accurate representation of the vehicles’ characteristics. Based on the linearization and order reduction, differential game theory and the sweep method are employed to derive and analytically solve the Riccati differential equation, yielding an optimal control strategy with an explicit expression. The theoretical rigor of the proposed approach is ensured through a proof of optimality sufficiency. Furthermore, factors influencing the computational accuracy of the Riccati equation solution, including the singular values of the control matrix and condition numbers of the solution matrix, are analyzed. Taking into account the dynamic response and limitations of the actuators, numerical simulations demonstrate the effectiveness and superiority of the proposed IGC strategy in intercepting the attacker and protecting the target compared to traditional separated guidance and control designs.

1. Introduction

Differential game theory, a branch of mathematical control theory, explores the strategic interactions in dynamic settings where players’ decisions influence the outcome over time [1]. One of the pivotal applications of differential games is in pursuit–evasion scenarios, which involve multiple agents with conflicting objectives. These scenarios are not only academically intriguing but have substantial practical implications in areas such as military tactics, security systems, and autonomous vehicle navigation [2,3,4].

Among the various configurations of pursuit–evasion games, three-player conflicts involving a Target, Attacker, and Defender (TAD) present a particularly complex and rich problem space [5,6,7]. In the classic two-player pursuit–evasion game, the scenario typically involves an attacker pursuing a target whose primary strategy is passive evasion. Instead, in TAD configurations the target benefits from the active protection provided by the defender, whose role is to intercept or obstruct the attacker’s efforts.

As a multi-agent system, with the inherent dynamics and interdependencies, the TAD conflict has garnered increasing attention due to its relevance to real-world scenarios in which valuable assets are actively defended [8]. Garcia et al. analyzed the game’s dynamics, focusing on the regions where the attacker prevails and the defender’s optimal interception strategies [9]. Their subsequent work in [10] developed optimal strategies for the attacker while ensuring that solutions were aligned with the Hamilton–Jacobi–Isaacs equation. A barrier analysis was introduced in [11] to delineate winning regions for both the attacker and the target–defender team, facilitating strategic decision-making through explicit policy and geometric analysis. Subsequent research explored stage-based strategies for a fast nonsuicidal attacker while adjusting tactics based on capture radii [12]. In [13], the authors presented a geometric interception strategy using the Apollonius circle, offering a practical approach for a slower defender to intercept a faster attacker.

The insights gained from TAD scenarios significantly enhance active protection systems for aerial vehicles, improving their effectiveness and survivability in hostile environments. The work in [14] optimized control inputs for the target and the defender, allowing effective interception even with slower defenders. The guidance law for defense missiles of nonmaneuverable aircraft was developed in [15], demonstrating superiority in non-coplanar engagements through high-fidelity simulations. Similar work was presented in [16], focusing on an air combat scenarios in which a defending missile aims to maximize the separation between a target aircraft and an attacking missile at the point of interception.

However, the aforementioned modeling approach in aerial pursuit–evasion games often simplifies the dynamics by directly using heading angles as control inputs, overlooking critical factors such as turning radii and the dynamic response of the system to commands, such as oscillations or delays. Additionally, these models tend to ignore the physical limitations of the maneuverability and actuators, which play a crucial role in the realistic implementation of control strategies. The oversimplified approach does not adequately leverage the general principles and knowledge of dynamics, and neglects the characteristics of aerial vehicles. As a result, the obtained strategies usually lack realistic feasibility. In fact, although there are nonlinear factors such as actuator saturation, aerial vehicles typically still exhibit pronounced linear characteristics, which form the classical foundation for the design of guidance and control systems [17,18]. It is essential for the design of guidance and control systems to include considerations of these linear characteristics along with the vehicle’s dynamic response processes.

This need for realistic modeling and design has led to the development of research that focuses on model linearization and order reduction. The prevailing research trend gravitates towards one-sided optimal control or differential games. The distinction lies in that optimal control theory-based guidance laws necessitate prior knowledge of the adversary’s control tactics. Although the reliance on prior information for one-sided optimization can be reduced through information sharing between the target and the defender [19], there are still challenges, such as the impracticality of deploying numerical optimization algorithms onboard [20,21]. Conversely, the differential game has attracted broader interest for its independence from presumptions regarding the opponent’s strategy [22]. The differential game can obtain the game strategy of the two opponents by finding the saddle point solution; under the condition of accurate modeling, this can guarantee the optimality of the strategy against the opponent’s arbitrary maneuver [23,24]. Considering the drawback that the control of linear quadratic differential game guidance law may exceed operational bounds, the bounded differential game was proposed and verified in a two-dimensional plane and three-dimensional space [25,26]. In addition, in order to be closer to the real combat scenarios, recent studies have accounted for constraints such as limited operational capabilities [27], state estimation under imperfect information through Kalman filtering [28], relative intercept angle limitations on attack protocols [24,29,30], cooperative multi-vehicle strategy against an active defense target [24,31], and weapon–target allocation strategies [32].

However, traditional active protection research has focused on designing guidance laws that generate acceleration commands without considering how these commands are effectively tracked by the autopilot. Such approaches overlook the control loop required to generate actuator commands. With the introduction of advanced guidance laws and increased confrontation, oversimplification can lead to significant distortion [33]. Moreover, in the endgame phase of engagements, where the relative motion between vehicles changes rapidly, the assumption of spectral separation between the guidance and control loops becomes untenable [34,35]. Therefore, integrated guidance and control (IGC) was proposed to address these issues by merging the guidance and control loops, thereby mitigating the negative impacts of their coupling [36,37]. IGC ensures that the commands sent to actuators are both feasible and effective within the physical constraints of the system. In the view of engineering designers, IGC is expected to significantly improve missile performance, resulting in lower weight and enhanced lethality [18]. In recent years, the development of IGC has progressed rapidly, with research increasingly focusing on a variety of advanced methods, including sliding mode control [38,39], back-stepping [40], robust control [41], and small-gain theorem [42], among others.

Building on the previously discussed challenges and advancements in three-player conflicts and IGC, this article aims to delve deeper into an IGC design specifically tailored for aerial vehicle active protection scenarios. We present an IGC model for active protection distinguished by the linearization of both the center of mass motion and rotation around the center of mass. Our approach leverages differential game theory and the sweep method to derive and solve the Riccati differential equation analytically. The main contributions of this paper are summarized as follows:

• Advanced IGC Modeling for Active Protection: An advanced modeling approach for active protection is introduced that uniquely combines linearization of both translational and rotational dynamics around the center of mass. Although the linearization technique is a well-established practice in aerial vehicle design, our contribution lies in applying this approach to the IGC design in a three-player conflict scenario. This modeling framework allows for a more comprehensive and accurate representation of the characteristics of aerial vehicles compared to previous studies.
• Analytical Derivation and Solution of Riccati Equation: Through the application of differential game theory and the sweep method, this study derives and solves the Riccati differential equation, providing an analytical expression for the optimal control strategy and in turn presenting possibilities for real-time onboard calculation.
• Theoretical Rigor in Solution Analysis: The theoretical rigor of the proposed approach is provided through a proof of optimality sufficiency. Additionally, we examine factors that influence the computational accuracy of the Riccati equation solution using singular values of the control matrix and condition numbers of the solution matrix.

The rest of the paper is organized as follows. In Section 2, we describe the three-player conflict in the IGC framework, focusing on model linearization and order reduction. Section 3 is dedicated to the derivation of the IGC active protection strategy. In Section 4, we conduct a theoretical analysis of the IGC strategy. Following this, Section 5 presents simulations to demonstrate the practical applicability and effectiveness of the strategy. Finally, Section 6 concludes the paper by summarizing the findings.

2. Problem Formulation

2.1. Nonlinear Engagement Model

In the TAD three-player conflict, the adversarial attacker aims at the target, which is protected by a defender aiming to intercept the incoming threat mid-course. In general, the target is weak in maneuvering and has difficulty avoiding being hit by the attacker through simple maneuvering strategies; thus, the target adopts an active defense strategy in order to survive.

As shown in Figure 1, the engagement is refined into a geometric representation within an inertial coordinate system $X_{I}O{Y}_I$. In the aerial conflict, each entity’s flight state is characterized by its flight path angle ($\gamma$), flight velocity (V), and actual acceleration (a), denoted by the respective subscripts A, T, and D for attacker, target, and defender. The engagement is divided into two distinct but interrelated processes: the attacker’s pursuit of the target ($AT$) and the defender’s interception of the attacker to protect the target ($AD$). In these engagements, r stands for the distance between the respective vehicles, $LOS$ denotes the line-of-sight, and $\lambda$ represents the line-of-sight angle. Specifically, $\dot{\gamma }$, $\dot{\lambda }$, and $\dot{r}$ denote the respective derivatives of $\gamma$, $\lambda$, and r. Initial conditions at the start of the engagement are indicated by the subscript 0.

As shown in Figure 1, the nonlinear engagement model of attacker–target and attacker–defender can be represented as

$$\begin{array}{cc}\hfill & {\dot{r}}_{Ai}=-V_{i}cos\left({\gamma }_i-{\lambda }_{Ai}\right)-V_{A}cos\left({\gamma }_A+{\lambda }_{Ai}\right)\hfill \\ \hfill & {\dot{\lambda }}_{Ai}=\frac{-V_{i}sin\left({\gamma }_i-{\lambda }_{Ai}\right)+V_{A}sin\left({\gamma }_A+{\lambda }_{Ai}\right)}{r_{Ai}},\hfill \end{array}$$

(1)

where A stands for the attacker and i represents the target or defender, i.e., $i\in \left\{T,D\right\}$.

The rate of change of the flight path angle for each vehicle can be expressed as

$${\dot{\gamma }}_j=\frac{a_j}{V_j}=\frac{F_j({\alpha }_j,{\delta }_j)}{m_j{V}_j},j\in \{A,T,D\},$$

(2)

where $m_j$ represents the mass of the vehicle j and $F_j$ is the normal aerodynamic force acting on the airframe, which is a function of the angle of attack ${\alpha }_j$ and the fin deflection ${\delta }_j$ from the actuator.

Unlike previous studies focusing solely on translational motion in active protection, the IGC design incorporates rotational kinematics and dynamics around the center of mass involving the pitch angle $\theta$ and pitch rate q. The angular acceleration ${\dot{q}}_j$ is generated by the torque applied to the airframe divided by the moment of inertia $I_j$. Thus, it is possible to obtain

$$\begin{array}{cc}\hfill \dot{\alpha }& =q-\frac{F_j({\alpha }_j,{\delta }_j)}{m_j{V}_j}\hfill \\ \hfill \ddot{\theta }& =\dot{q}=\frac{M_j({\alpha }_j,{\delta }_j)}{I_j}.\hfill \end{array}$$

(3)

Although the nonlinear differential Equations (1)–(3) can be solved numerically, adopting an analytical approach provides clearer insights and understanding of the engagement dynamics. To facilitate this, we linearize the equations around an operating condition, which allows for analysis and design based on linear systems theory.

2.2. Linearization and Order Reduction

To facilitate analysis of the engagement dynamics in the terminal phase, the model is simplified based on the following set of assumptions:

• Constant Velocity: The velocities of the entities are assumed to be constant.
• Co-planar Movement: The three entities are considered to move within the same plane; the defender is launched by the target, meaning that $LO{S}_{MT0}$ and $LO{S}_{MD0}$ coincide.
• Linear Plant Models: The dynamics of each entity are approximated as linear systems.
• Linearization of Collision Triangles: It is assumed that the two collision triangles can be linearized along their respective initial lines of sight (LOS).

The dynamics model of each vehicle can be represented by a linear equation of arbitrary order:

$$\begin{array}{c}{\mathit{x}}_j={\mathit{A}}_j{\mathit{x}}_j+{\mathit{B}}_j{u}_j,j\in \left\{A,T,D\right\}\hfill \\ a_j={\mathit{C}}_j{\mathit{x}}_j+d_j{u}_j,j\in \left\{A,T,D\right\}\hfill \end{array}$$

(4)

where ${\mathit{x}}_j$ represents the internal state variables of each vehicle, $a_j$ is the acceleration perpendicular to the initial line-of-sight direction, and $u_j$ is the corresponding control input.

In previous studies, vehicles described by Equation (4) have generally been simplified to ideal dynamics or first-order delay dynamics. Ideal dynamics implies that ${\mathit{A}}_j={\mathit{B}}_j={\mathit{C}}_j=0$ and $d_j=1$. First-order delay dynamics implies that ${\mathit{A}}_j=1/{\tau }_j$, ${\mathit{B}}_j=1/{\tau }_j$, ${\mathit{C}}_j=1/{\tau }_j$, and $d_j=0$ with the time constant of ${\tau }_j$. In both cases, $u_j$ represents the command acceleration. We further extend the scope of Equation (4) to cover the short-period dynamics described by Equation (3) for the target and the defender [43], which means that

$$\begin{array}{cccc}\hfill & {\mathit{A}}_i=\left[\begin{array}{cc}{Z}_{\alpha }^{\left(i\right)}& 1\\ M_{\alpha }^{\left(i\right)}& 0\end{array}\right]\hfill & \hfill & {\mathit{B}}_i=\left[\begin{array}{c}{Z}_{\delta }^{\left(i\right)}\\ M_{\delta }^{\left(i\right)}\end{array}\right]\hfill \\ \hfill & {\mathit{C}}_i=\left[\begin{array}{cc}-Z_{\alpha }^{\left(i\right)}{V}^{\left(i\right)}& 0\end{array}\right]\hfill & \hfill & d_i=-Z_{\delta }^{\left(i\right)}{V}^{\left(i\right)}i\in \left\{T,D\right\},\hfill \end{array}$$

(5)

where the numerical coefficients are defined by

$$\begin{array}{cccc}\hfill & Z_{\alpha }^{\left(i\right)}=-\frac{1}{m_i{V}_i}\frac{\partial F_i\left({\alpha }_i,{\delta }_i\right)}{\partial {\alpha }_i}\hfill & \hfill & Z_{\delta }^{\left(i\right)}=-\frac{1}{m_i{V}_i}\frac{\partial F_i\left({\alpha }_i,{\delta }_i\right)}{\partial {\delta }_i},\hfill \\ \hfill & M_{\alpha }^{\left(i\right)}=\frac{1}{I_i}\frac{\partial M_i({\alpha }_i,{\delta }_i)}{\partial {\alpha }_i}\hfill & \hfill & M_{\delta }^{\left(i\right)}=\frac{1}{I_i}\frac{\partial M_i({\alpha }_i,{\delta }_i)}{\partial {\delta }_i}.\hfill \end{array}$$

(6)

The target and defender are modeled as IGC configurations and the attacker is modeled with first-order lag dynamics. Choosing the state variables as ${\mathit{x}}_{Ai}={\left[y_{Ai},{\dot{y}}_{Ai},{\alpha }_i,q_i,a_A\right]}^{\top }$, the equations of motion for the attacker to engage the target or defender can be expressed as

$${\dot{\mathit{x}}}_{Ai}={\mathit{A}}_{Ai}\mathit{x}+{\tilde{\mathit{B}}}_i{u}_i+{\tilde{\mathit{B}}}_A{u}_A,i\in \left\{T,D\right\}$$

(7)

where

$$\begin{array}{cc}\hfill & {\mathit{A}}_{Ai}=\left[\begin{array}{ccccc}0& 1& 0& 0& 0\\ 0& 0& \mp Z_{\alpha }^{\left(i\right)}{V}_i& 0& \mp 1\\ 0& 0& Z_{\alpha }^{\left(i\right)}& 1& 0\\ 0& 0& M_{\alpha }^{\left(i\right)}& 0& 0\\ 0& 0& 0& 0& -1/{\tau }_A\end{array}\right]\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{B}}}_T={\left[\begin{array}{ccccc}0& -Z_{\delta }^{\left(T\right)}{V}_T& Z_{\delta }^{\left(T\right)}& M_{\delta }^{\left(T\right)}& 0\end{array}\right]}^{\top }\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{B}}}_D={\left[\begin{array}{ccccc}0& Z_{\delta }^{\left(D\right)}{V}_D& Z_{\delta }^{\left(D\right)}& M_{\delta }^{\left(D\right)}& 0\end{array}\right]}^{\top }\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{B}}}_A={\left[\begin{array}{ccccc}0& 0& 0& 0& 1/{\tau }_A\end{array}\right]}^{\top }.\hfill \end{array}$$

(11)

The symbols ∓ in Equation (8) correspond to $i=T$ and $i=D$, respectively. For the target and the defender, the control input $u_i$ is defined as the fin deflection of the actuator ${\delta }_i$. For the attacker, the control input $u_A$ is the commanded acceleration $a_A^c$.

The time-to-go between the attacker–target pair and the defender–attacker pair is denoted by $t_{go1}$ and $t_{go2}$, respectively. The time-to-go can be calculated by $t_{go1}=t_{f1}-t$ and $t_{go2}=t_{f2}-t$, where the interception time is defined as

$$\begin{array}{cc}\hfill t_{f1}& =\frac{r_{AT0}}{\left[V_{A}cos\left({\gamma }_{A0}+{\lambda }_{AT0}\right)+V_{T}cos\left({\gamma }_{T0}-{\lambda }_{AT0}\right)\right]},\hfill \\ \hfill t_{f2}& =\frac{r_{AD0}}{\left[V_{A}cos\left({\gamma }_{A0}+{\lambda }_{AD0}\right)+V_{D}cos\left({\gamma }_{D0}-{\lambda }_{AD0}\right)\right]}.\hfill \end{array}$$

(12)

We assume that the engagement of the attacker with the defender precedes the engagement of the attacker with the target, i.e., $t_{f1}$ and $t_{f2}$ satisfy $t_{f1}-t_{f2}>0$ in the timeline. This is because once the attacker hits or misses the target, it means that the game is over and the defender is no longer needs to continue the engagement.

Therefore, unlike previous studies that only considered guidance loops, the linearization encompasses two aspects: linearization of the airframe and linearization of engagement geometry.

• Linearization of the airframe: Linearization of the airframe dynamics is typically performed around a trim condition, which represents a steady-state flight condition; in our case, this corresponds to flight at a constant velocity and a small angle of attack. This approach is widely accepted in aerospace control systems design, especially for short-period dynamics, as deviations from this trim condition during maneuvers are usually small enough to maintain the validity of the linear model.
• Linearization of engagement geometry: For the engagement geometry, we employ linearization around a nominal collision course. This approach is commonly used and justified in terminal guidance problems, particularly in the endgame phase, where the relative geometry changes are relatively small and can be approximated linearly.

This situation can be formulated as a three-player differential game between two opposing sides with different objectives. The cost function that we utilize incorporates both miss distance and energy consumption through a linear–quadratic (LQ) form, striking a balance between the two factors. This also means that it is a soft constraint, ensuring that the control effort does not exceed the maneuverability. The specific formulation of the cost function is as follows:

$$\underset{\left[u_T,u_D\right]}{min}\underset{u_A}{max}J=-\frac{1}{2}{\alpha }_{AT}{y}_{AT}^2\left(t_{f1}\right)+\frac{1}{2}{\alpha }_{AD}{y}_{AD}^2\left(t_{f2}\right)+\frac{1}{2}{\int }_0^{t_{f1}}{\beta }_T{u}_T^2\mathrm{d}t+\frac{1}{2}{\int }_0^{t_{f2}}{\beta }_D{u}_D^2\mathrm{d}t-\frac{1}{2}{\int }_0^{t_{f1}}{u}_A^2\mathrm{d}t$$

(13)

where ${\alpha }_{AT}$ and ${\alpha }_{AD}$ are non-negative weight coefficients related to miss distance and the coefficients ${\beta }_U$ and ${\beta }_D$ reflect the maneuverability of the target and defender relative to the attacker. When the target has weaker maneuverability compared to the attacker, ${\beta }_T$ takes a larger value; ${\beta }_T\to \infty$ indicates that the target has no maneuverability. When the defender has stronger maneuverability compared to the attacker, ${\beta }_D$ takes a smaller value, with ${\beta }_D\to 0$ signifying strong maneuverability of the defender.

It is still difficult to obtain the analytical solution of the above linear optimal control problem, so reducing the model’s order is essential. In the context of guidance analysis, the zero-effort miss (ZEM) is a commonly used concept, which physically represents the miss distance of the interceptor if no maneuver is performed, i.e., the control command is zero from the current moment. From the perspective of signal processing, this is equivalent to the zero-input response of the linear system. Using the terminal projection transformation of the linear system, the ZEM $z_1$ and $z_2$ can be expressed as

$$\begin{array}{cc}\hfill & z_1\left(t\right)=\mathit{L}{\mathbf{\Phi }}_{AT}(t_{f1},t){\mathit{x}}_{AT},\hfill \\ \hfill & z_2\left(t\right)=\mathit{L}{\mathbf{\Phi }}_{AD}(t_{f2},t){\mathit{x}}_{AD},\hfill \end{array}$$

(14)

where $\mathit{L}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\end{array}\right]$ is a coefficient matrix and ${\mathbf{\Phi }}_{Ai}(t_f,t)$ is the state transition matrix.

The derivatives of $z_1$ and $z_2$ with respect to time in compact matrix form are

$$\begin{array}{cc}\hfill \dot{\mathit{z}}=\left[\begin{array}{c}{\dot{z}}_1\left(t\right)\hfill \\ {\dot{z}}_2\left(t\right)\hfill \end{array}\right]& =G_1\left[\begin{array}{c}{u}_T\hfill \\ u_D\hfill \end{array}\right]+G_2{u}_A\hfill \\ \hfill & =G_1{u}_{T/D}+G_2{u}_A,\hfill \end{array}$$

(15)

where

$$\begin{array}{cc}\hfill & G_1=\left[\begin{array}{c}\hfill {\Phi }_{AT,y}\left(t_{go1}\right)\\ \hfill {\Phi }_{AD,y}\left(t_{go2}\right)\end{array}\right]·\left[\begin{array}{cc}{\tilde{\mathit{B}}}_T& {\tilde{\mathit{B}}}_D\end{array}\right]\hfill \\ \hfill & G_2=\left[\begin{array}{c}{\Phi }_{AT,y}\left(t_{go1}\right)\hfill \\ {\Phi }_{AD,y}\left(t_{go2}\right)\hfill \end{array}\right]·{\tilde{\mathit{B}}}_A\hfill \\ \hfill & {\Phi }_{AT,y}\left(t_{go1}\right)=\mathit{L}{\mathbf{\Phi }}_{AT}(t_{f1},t)\hfill \\ \hfill & {\Phi }_{AD,y}\left(t_{go2}\right)=\mathit{L}{\mathbf{\Phi }}_{AD}(t_{f2},t).\hfill \end{array}$$

(16)

Correspondingly, the cost function of the reduced-order differential game problem is provided by

$$\underset{u_{T/D}}{min}\underset{u_A}{max}{J}_z=\frac{1}{2}{\mathit{z}}^{\top }\left(t_f\right)Q\mathit{z}\left(t_f\right)+\frac{1}{2}{\int }_{t_0}^{t_f}\left(u_{T/D}^{\top }{R}_1{u}_{T/D}+u_A^{\top }{R}_2{u}_A\right)dt,$$

(17)

where

$$Q=\left[\begin{array}{cc}-{\alpha }_{AT}& 0\\ 0& {\alpha }_{AD}\end{array}\right]R_1=\left[\begin{array}{cc}{\beta }_T& 0\\ 0& {\beta }_D\end{array}\right]R_2=-1$$

3. Derivation of IGC Active Protection Strategy

The differential game approach is applied to solve for the optimal control of each vehicle. The Hamiltonian of the linearized system is

$$H=\frac{1}{2}\left(u_{T/D}^{\top }{R}_1{u}_{T/D}+u_A^{\top }{R}_2{u}_A\right)+{\mathit{\lambda }}^{\top }\left(G_1{u}_{T/D}+G_2{u}_A\right),$$

(18)

where $\mathit{\lambda }$ stands for the introduced costates. According to the calculus of variations, by equating the derivatives of the Hamiltonian with respect to the players’ controllers to zero we obtain

$$\begin{array}{cc}\hfill & \frac{\partial H}{\partial u_{T/D}}=R_1{u}_{T/D}+G_1^{\top }\mathit{\lambda }=0\Rightarrow u_{T/D}=-R_1^{-1}{G}_1^{\top }\mathit{\lambda },\hfill \\ \hfill & \frac{\partial H}{\partial u_A}=R_2{u}_A+G_2^{\top }\mathit{\lambda }=0\Rightarrow u_A=-R_2^{-1}{G}_2^{\top }\mathit{\lambda }.\hfill \end{array}$$

(19)

It follows that in order to solve for the expression of the optimal control we need to obtain an expression for the costate $\mathit{\lambda }$. From the adjoint equation and the transversality condition, we have

$$\dot{\mathit{\lambda }}=-\frac{\partial H}{\partial \mathit{z}}=0\mathit{\lambda }\left(t_f\right)=\left[\begin{array}{c}-{\alpha }_{AT}{z}_1\left(t_{f1}\right)\\ {\alpha }_{AD}{z}_2\left(t_{f2}\right)\end{array}\right]=\left[\begin{array}{c}{\lambda }_1\left(t\right)\hfill \\ {\lambda }_2\left(t\right)\hfill \end{array}\right].$$

(20)

As shown in Equation (20), the expression for $\mathit{\lambda }$ includes $z_1\left(t_{f1}\right)$ and $z_2\left(t_{f2}\right)$. However, unlike previous studies such as [22] which tackled more idealized dynamic models, the complexity inherent to the IGC problem precludes a straightforward analytical solution. Therefore, to address these challenges we adopt the sweep method to solve for the costates. We assume that $\mathit{\lambda }\left(t\right)$ can be written in the form of linear feedback of reduced-order states as

$$\mathit{\lambda }\left(t\right)=P\left(t\right)\mathit{z}\left(t\right),$$

(21)

i.e.,

$$\left[\begin{array}{c}{\lambda }_1\left(t\right)\\ {\lambda }_2\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{P}_{11}\left(t\right)& P_{12}\left(t\right)\\ P_{21}\left(t\right)& P_{22}\left(t\right)\end{array}\right]\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right],$$

(22)

with the transversality condition

$$\mathit{\lambda }\left(t_f\right)=Q\mathit{z}\left(t_f\right).$$

(23)

Taking the derivative of $\mathit{\lambda }$, as shown in Equation (21), yields

$$\dot{\mathit{\lambda }}=\dot{P}\mathit{z}+P\dot{\mathit{z}}=0.$$

(24)

Next, incorporating the reduced-order dynamics of the system, represented by Equation (15), results in

$$\dot{P}\mathit{z}+P\dot{\mathit{z}}=\dot{P}\mathit{z}+P\left[-G_1\left(t\right)R_1^{-1}{G}_1^{\top }\left(t\right)P\left(t\right)-G_2\left(t\right)R_2^{-1}{G}_2^{\top }\left(t\right)P\left(t\right)\right]\mathit{z}\left(t\right)=0.$$

(25)

Upon rearranging and simplifying the terms, we obtain the Riccati differential equation

$$\dot{P}\left(t\right)=P\left(t\right)\left[G_1\left(t\right)R_1^{-1}{G}_1^{\top }\left(t\right)P\left(t\right)+G_2\left(t\right)R_2^{-1}{G}_2^{\top }\left(t\right)P\left(t\right)\right].$$

(26)

Equation (26) characterizes the evolution of the matrix $P\left(t\right)$ over time. The terms within the brackets represent the combined influence of the system dynamics and the control effort costs on the evolution of $P\left(t\right)$. However, the $P\left(t\right)$ matrix in Equation (26) has four elements and the variables are coupled, which makes it impossible to solve it directly. In order to find the analytical expression for $P\left(t\right)$, the matrix $K\left(t\right)$ is introduced to satisfy

$$K\left(t\right)=P^{-1}\left(t\right)K\left(t\right)P\left(t\right)=I,$$

(27)

where I is the identity matrix. Therefore, the derivative of $K\left(t\right)$ can be expressed in terms of $P\left(t\right)$:

$$\dot{K}\left(t\right)P\left(t\right)+K\left(t\right)\dot{P}\left(t\right)=0.$$

(28)

This leads to the fully decoupled matrix differential equation for $K\left(t\right)$:

$$\begin{array}{cc}\hfill \dot{K}\left(t\right)& =-K\left(t\right)\dot{P}\left(t\right)P^{-1}\left(t\right)\hfill \\ \hfill & =-\underset{I}{\underbrace{K\left(t\right)P\left(t\right)}}{G}_1\left(t\right)R_1^{-1}{G}_1^{\top }\left(t\right)-\underset{I}{\underbrace{K\left(t\right)P\left(t\right)}}{G}_2\left(t\right)R_2^{-1}{G}_2^{\top }\left(t\right)\hfill \\ \hfill & =-G_1\left(t\right)R_1^{-1}{G}_1^{\top }\left(t\right)-G_2\left(t\right)R_2^{-1}{G}_2^{\top }\left(t\right).\hfill \end{array}$$

(29)

This decoupled differential equation for $K\left(t\right)$ is more tractable and its terms purely depend on the matrices $G_1\left(t\right)$, $G_2\left(t\right)$, $R_1$, and $R_2$, which are typically known from the system dynamics and cost function. The terminal condition for the matrix $P\left(t\right)$ can be determined from the costate’s terminal condition stated in (23):

$$P\left(t_f\right)=Q=\left[\begin{array}{cc}-{\alpha }_{AT}& 0\\ 0& {\alpha }_{AD}\end{array}\right].$$

(30)

Correspondingly, the terminal condition for $K\left(t\right)$ is provided by

$$K\left(t_f\right)=P^{-1}\left(t_f\right)=\left[\begin{array}{cc}-1/{\alpha }_{AT}& 0\\ 0& 1/{\alpha }_{AD}\end{array}\right].$$

(31)

This setup allows for backward integration from $t_f$ to the initial time while using $K\left(t_f\right)$ to calculate $K\left(t\right)$, providing a method for determining $P\left(t\right)$ across the engagement timeline.

In light of the computational challenges posed by backward integration and the specific representation of the state transition matrix in terms of time-to-go, it is useful to reverse the time in the differential equation for the matrix $K\left(t\right)$. Moreover, reversing the time allows for dynamic adjustment of control inputs based on changes in the engagement through state feedback, rather than relying on a precomputed offline $P\left(t\right)$. This adjustment leads to

$$\frac{\mathrm{d}K}{\mathrm{d}{t}_{go}}=-\frac{\mathrm{d}K}{\mathrm{d}t}=G_1\left(t\right)R_1^{-1}{G}_1^{\top }\left(t\right)+G_2\left(t\right)R_2^{-1}{G}_2^{\top }\left(t\right).$$

(32)

Simultaneously, the terminal conditions in (31) are converted into initial conditions for the integration process in terms of $t_{go}$. Thus, the problem is transformed into an initial value problem for system ordinary differential equations in a decoupled form. The decoupling facilitates solving for $K\left(t\right)$ directly through indefinite integration, which in turn allows for determining the solution $P\left(t\right)$ of the Riccati differential Equation (26).

Furthermore, in order to obtain the analytical expression of the optimal control it is necessary to solve for the relevant elements in Equation (16). Because ${\mathit{A}}_{Ai}$ in (7) is time-invariant, the state transition matrix can be solved by ${\Phi }_{A\mathit{i}}={\mathcal{L}}^{-1}\left[{\left(s\mathit{I}-{\mathit{A}}_{Ai}\right)}^{-1}\right]$. To find the inverse of $\left(s\mathit{I}-{\mathit{A}}_{Ai}\right)$, we can augment the given matrix with the identity matrix and perform Gaussian elimination. By performing elementary row and column operations on the augmented matrix, we obtain

$$\mathit{L}·{\left(s\mathit{I}-{\mathit{A}}_{Ai}\right)}^{-1}=\left[\begin{array}{ccccc}{\displaystyle \frac{1}{s}}& {\displaystyle \frac{1}{s^2}}& \pm {\displaystyle \frac{V_i{Z}_{\alpha }^{\left(i\right)}}{s(M_{\alpha }^{\left(i\right)}+s(-s+Z_{\alpha }^{\left(i\right)}))}}& \pm {\displaystyle \frac{V_i{Z}_{\alpha }^{\left(i\right)}}{s^2(M_{\alpha }^{\left(i\right)}+s(-s+Z_{\alpha }^{\left(i\right)}))}}& \mp {\displaystyle \frac{{\tau }_A}{s^2+s^3{\tau }_A}}\end{array}\right].$$

(33)

Thus, ${\Phi }_{AT,y}\left(t_{go1}\right)$ and ${\Phi }_{AD,y}\left(t_{go2}\right)$ in $G_1$ and $G_2$ can be given by

$$\begin{array}{cc}\hfill {\Phi }_{AT,y}\left(t_{go1}\right)& ={\mathcal{L}}^{-1}\left[\mathit{L}·{\left(s\mathit{I}-{\mathit{A}}_{AT}\right)}^{-1}\right],\hfill \\ \hfill {\Phi }_{AD,y}\left(t_{go2}\right)& ={\mathcal{L}}^{-1}\left[\mathit{L}·{\left(s\mathit{I}-{\mathit{A}}_{AD}\right)}^{-1}\right].\hfill \end{array}$$

(34)

In order to solve the inverse Laplace transform of the above equation, the following transformation is provided:

$$\frac{Z_{\alpha }}{M_{\alpha }+Z_{\alpha }s-s^2}=-\frac{Z_{\alpha }}{{\kappa }_1}·\frac{{\kappa }_1}{{\left(s+{\kappa }_2\right)}^2+{\kappa }_1^2}$$

(35)

$$\frac{Z_{\alpha }}{s\left(M_{\alpha }+Z_{\alpha }s-s^2\right)}=\frac{Z_{\alpha }}{M_{\alpha }}\left[\frac{1}{s}-\frac{s+{\kappa }_2}{{\left(s+{\kappa }_2\right)}^2+{\kappa }_1^2}+\frac{Z_{\alpha }}{2a}·\frac{{\kappa }_1}{{\left(s+{\kappa }_2\right)}^2+{\kappa }_1^2}\right]$$

(36)

$$\frac{Z_{\alpha }}{s^2\left(M_{\alpha }+Z_{\alpha }s-s^2\right)}=\frac{Z_{\alpha }}{M_{\alpha }}\left[\frac{1}{s^2}-\frac{Z_{\alpha }}{M_{\alpha }s}+\frac{Z_{\alpha }}{M_{\alpha }}\frac{s+{\kappa }_2}{{\left(s+{\kappa }_2\right)}^2+{\kappa }_1^2}-\frac{1+\frac{Z_{\alpha }^2}{2M_{\alpha }}}{{\kappa }_1}·\frac{{\kappa }_1}{{\left(s+{\kappa }_2\right)}^2+{\kappa }_1^2}\right]$$

(37)

where ${\kappa }_1=\sqrt{-M_{\alpha }-\frac{Z_{\alpha }^2}{4}}$, ${\kappa }_2=-\frac{Z_{\alpha }}{2}$. It should be noted that Equations (35)–(37) are not a partial fraction decomposition of the original fraction, as we have not found the roots of the denominator. The rational fractions in Equations (35)–(37) allow us to utilize the following inverse Laplace transform:

$${\mathcal{L}}^{-1}\left[\frac{a}{{\left(s+b\right)}^2+a^2}\right]=e^{-bt}sinat{\mathcal{L}}^{-1}\left[\frac{s+b}{{\left(s+b\right)}^2+a^2}\right]=e^{-bt}cosat.$$

(38)

Let ${\kappa }_1^{\left(T\right)}=\sqrt{-M_{\alpha }^{\left(T\right)}-\frac{{\left(Z_{\alpha }^{\left(T\right)}\right)}^2}{4}}$, ${\kappa }_2^{\left(T\right)}=-\frac{Z_{\alpha }^{\left(T\right)}}{2}$ and ${\kappa }_1^{\left(D\right)}=\sqrt{-M_{\alpha }^{\left(D\right)}-\frac{{\left(Z_{\alpha }^{\left(D\right)}\right)}^2}{4}}$, ${\kappa }_2^{\left(D\right)}=-\frac{Z_{\alpha }^{\left(D\right)}}{2}$; then, the first row of the state transition matrix used to compute the ZEM ${\Phi }_{AT,y}\left(t,t_{f1}\right)$ and ${\Phi }_{AD,y}\left(t,t_{f2}\right)$ can be derived from the inverse Laplace transform lookup table, which can be represented as

$$\begin{array}{cc}\hfill & {\Phi }_{AT,y}\left(t,t_{f1}\right)={\Phi }_{AT,y}\left(t_{go1}\right)\hfill \\ \hfill =& {\left[\begin{array}{c}1\\ t_{go1}\\ \frac{Z_{\alpha }^{\left(T\right)}{V}_T}{M_{\alpha }^{\left(T\right)}}\left(1-e^{-{\kappa }_2^{\left(T\right)}{t}_{go1}}cos\left({\kappa }_1^{\left(T\right)}{t}_{go1}\right)+\frac{e^{-{\kappa }_2^{\left(T\right)}{t}_{go1}}{Z}_{\alpha }^{\left(T\right)}}{2{\kappa }_1^{\left(T\right)}}sin\left({\kappa }_1^{\left(T\right)}{t}_{go1}\right)\right)\\ \frac{Z_{\alpha }^{\left(T\right)}{V}_T}{M_{\alpha }^{\left(T\right)}}\left(t_{go1}-\frac{Z_{\alpha }^{\left(T\right)}}{M_{\alpha }^{\left(T\right)}}+\frac{Z_{\alpha }^{\left(T\right)}}{M_{\alpha }^{\left(T\right)}}{e}^{-{\kappa }_2^{\left(T\right)}{t}_{go1}}cos\left({\kappa }_1^{\left(T\right)}{t}_{go1}\right)-\frac{1+\frac{{\left(Z_{\alpha }^{\left(T\right)}\right)}^2}{2M_{\alpha }^{\left(T\right)}}}{{\kappa }_1^{\left(T\right)}}{e}^{-{\kappa }_2^{\left(T\right)}{t}_{go1}}sin\left({\kappa }_1^{\left(T\right)}{t}_{go1}\right)\right)\\ -{\tau }_A\left(t_{go1}-{\tau }_A+e^{-\frac{t_{go1}}{{\tau }_A}}{\tau }_A\right)\end{array}\right]}^{\top }\hfill \end{array}$$

(39)

and

$$\begin{array}{cc}\hfill & {\Phi }_{AD,y}\left(t,t_{f2}\right)={\Phi }_{AD,y}\left(t_{go2}\right)\hfill \\ \hfill =& {\left[\begin{array}{c}1\\ t_{go2}\\ -\frac{Z_{\alpha }^{\left(D\right)}{V}_D}{M_{\alpha }^{\left(D\right)}}\left(1-e^{-{\kappa }_2^{\left(D\right)}{t}_{go2}}cos\left({\kappa }_1^{\left(D\right)}{t}_{go2}\right)+\frac{e^{-{\kappa }_2^{\left(D\right)}{t}_{go2}}{Z}_{\alpha }^{\left(D\right)}}{2{\kappa }_1^{\left(D\right)}}sin\left({\kappa }_1^{\left(D\right)}{t}_{go2}\right)\right)\\ -\frac{Z_{\alpha }^{\left(D\right)}{V}_D}{M_{\alpha }^{\left(D\right)}}\left(t_{go2}-\frac{Z_{\alpha }^{\left(D\right)}}{M_{\alpha }^{\left(D\right)}}+\frac{Z_{\alpha }^{\left(D\right)}}{M_{\alpha }^{\left(D\right)}}{e}^{-{\kappa }_2^{\left(D\right)}{t}_{go2}}cos\left({\kappa }_1^{\left(D\right)}{t}_{go2}\right)-\frac{1+\frac{{\left(Z_{\alpha }^{\left(D\right)}\right)}^2}{2M_{\alpha }^{\left(D\right)}}}{{\kappa }_1^{\left(D\right)}}{e}^{-{\kappa }_2^{\left(D\right)}{t}_{go2}}sin\left({\kappa }_1^{\left(D\right)}{t}_{go2}\right)\right)\\ {\tau }_A\left(t_{go2}-{\tau }_A+e^{-\frac{t_{go2}}{{\tau }_A}}{\tau }_A\right)\end{array}\right]}^{\top }.\hfill \end{array}$$

(40)

Substituting ${\Phi }_{AT,y}\left(t_{go1}\right)$ and ${\Phi }_{AD,y}\left(t_{go2}\right)$ into $G_1$ and $G_2$ in (32), it can be found that the elements of $\frac{dK}{d{t}_{go}}$ share the following form and that the indefinite integral can be expressed in terms of elementary functions:

$$\begin{array}{cc}\hfill & \int {\left(n+tp-r{e}^{tq}cos\left(mt\right)+v{e}^{tq}sin\left(mt\right)\right)}^{2}dt\hfill \\ \hfill & =\frac{1}{12}\left(\begin{array}{c}4t\left(3n^2+3npt+p^2{t}^2\right)+\frac{3e^{2qt}(r^2+v^2)}{q}\\ +\frac{3e^{2qt}\left(2mrv+q(r^2-v^2)\right)cos\left(2mt\right)}{m^2+q^2}-\frac{3e^{2qt}\left(2qrv+m(-r^2+v^2)\right)sin\left(2mt\right)}{m^2+q^2}\\ -\frac{24e^{qt}\left(m^{2}r\left(p+nq+pqt\right)+q^{2}r\left(nq+p(-1+qt)\right)+m^3\left(n+pt\right)v+mq\left(nq+p(-2+qt)\right)v\right)cos\left(mt\right)}{{(m^2+q^2)}^2}\\ +\frac{24e^{qt}\left(-m^{3}r\left(n+pt\right)-mqr\left(nq+p(-2+qt)\right)+m^2\left(p+nq+pqt\right)v+q^2\left(nq+p(-1+qt)\right)v\right)sin\left(mt\right)}{{(m^2+q^2)}^2}\end{array}\right).\hfill \end{array}$$

(41)

This indefinite integral can be obtained by expanding the square and employing general integration techniques.

At this point, we can obtain the expression of $K\left(t\right)$ directly through an indefinite integral, thereby obtaining the solution of the Riccati equation.

4. Theoretical Analysis of IGC Active Protection Strategy

4.1. Saddle Point Sufficiency Proof

In the development of IGC strategies, establishing the optimality of a control solution is crucial. This not only ensures the effectiveness of the control in achieving its goals, thereby minimizing miss distance and control effort, it provides theoretical validation that the proposed solution indeed performs better than any other strategy under the given model assumptions. To this end, proving that a solution to the Riccati equation serves as a sufficient condition for the existence of an optimal saddle point solution is essential. Such a proof confirms that the derived control law is not just a feasible solution but the best possible strategy within the framework of the IGC problem [44].

Theorem 1.

The existence of a solution to the Riccati equation (Equation (26)) is a sufficient condition for the existence of a saddle point solution, and the control derived from the Riccati equation is optimal.

Proof. 

Consider the reduced-order active protection problem defined by

$$\underset{u_{T/D}}{min}\underset{u_A}{max}{J}_z=\frac{1}{2}{\mathit{z}}^{\top }\left(t_f\right)Q\mathit{z}\left(t_f\right)+\frac{1}{2}{\int }_{t_0}^{t_f}\left(u_{T/D}^{\top }{R}_1{u}_{T/D}+u_M^{\top }{R}_2{u}_M\right)dt$$

(42)

subject to the dynamics

$$\dot{\mathit{z}}\left(t\right)=G_1{u}_{T/D}+G_2{u}_M.$$

(43)

To analyze this, we construct the identity

$$S={\mathit{z}}_0^{\top }P\left(t_0\right){\mathit{z}}_0-{\mathit{z}}^{\top }\left(t_f\right)Q_f\mathit{z}\left(t_f\right)+{\int }_0^{t_f}\frac{d}{dt}\left({\mathit{z}}^{\top }P\mathit{z}\right)dt\equiv 0.$$

(44)

Using the chain rule for differentiation, we obtain

$$\frac{d}{dt}\left({\mathit{z}}^{\top }P\mathit{z}\right)=\left({\dot{\mathit{z}}}^{\top }P\mathit{z}\right)+\left({\mathit{z}}^{\top }\dot{P}\mathit{z}\right)+\left({\mathit{z}}^{\top }P\dot{\mathit{z}}\right).$$

(45)

Substituting and expanding within S, we obtain

$$\begin{array}{cc}\hfill S& ={\mathit{z}}_0^{\top }P\left(0\right){\mathit{z}}_0-{\mathit{z}}^{\top }\left(t_f\right)Q_f\mathit{z}\left(t_f\right)-{\int }_0^{t_f}\left({\mathit{z}}^{\top }Q\mathit{z}+u_{T/D}^{\top }{R}_1{u}_{T/D}+u_M^{\top }{R}_2{u}_M\right)dt\hfill \\ \hfill & +{\int }_0^{t_f}{\left(u_{T/D}+R_1^{-1}{G}_1^{\top }P\mathit{z}\right)}^{\top }{R}_1\left(u_{T/D}+R_1^{-1}{G}_1^{\top }P\mathit{z}\right)dt\hfill \\ \hfill & -{\int }_0^{t_f}{\left(u_M+R_2^{-1}{G}_2^{\top }P\mathit{z}\right)}^{\top }{R}_2\left(u_M+R_2^{-1}{G}_2^{\top }P\mathit{z}\right)dt.\hfill \end{array}$$

(46)

From this, we can deduce that the performance index $J_z$ satisfies

$$\begin{array}{cc}\hfill 2J_Z& ={\mathit{z}}_0^{\top }P\left(0\right){\mathit{z}}_0+{\int }_0^{t_f}{\left(u_{T/D}+R_1^{-1}{G}_1^{\top }P\mathit{z}\right)}^{\top }{R}_1\left(u_{T/D}+R_1^{-1}{G}_1^{\top }P\mathit{z}\right)dt\hfill \\ \hfill & -{\int }_0^{t_f}{\left(u_M+R_2^{-1}{G}_2^{\top }P\mathit{z}\right)}^{\top }{R}_2\left(u_M+R_2^{-1}{G}_2^{\top }P\mathit{z}\right)dt\hfill \\ \hfill & ={\mathit{z}}_0^{\top }P\left(0\right){\mathit{z}}_0+{\int }_0^{t_f}{\left(u_{T/D}-u_{T/D}^{*}\right)}^{\top }{R}_1\left(u_{T/D}-u_{T/D}^{*}\right)dt-{\int }_0^{t_f}{\left(u_M-u_M^{*}\right)}^{\top }{R}_2\left(u_M-u_M^{*}\right)dt.\hfill \end{array}$$

(47)

The expression clearly shows that $J_Z(u_{T/D}^{*},u_M)\le J_Z(u_{T/D}^{*},u_M^{*})\le J_Z(u_{T/D},u_M^{*})$, indicating that the performance index is minimized for the optimal controls $u_{T/D}^{*}$ and $u_M^{*}$, confirming the saddle point nature of the solution. This completes the proof. □

4.2. Analysis of the Invertibility of Matrix $K\left(t\right)$

The invertibility of the matrix $K\left(t\right)$ is critical in solving the Riccati equation and ensuring the computation accuracy of the proposed strategy. The ill-conditioning of $K\left(t\right)$ significantly impacts the numerical stability, especially during matrix inversion operations. Ill-conditioning refers to scenarios where small variations in inputs cause large changes in outputs, leading to significant numerical errors in calculations such as matrix inversion. An ill-conditioned $K\left(t\right)$ can undermine the reliability of computing its inverse.

It can be seen through Equation (32) that $K\left(t\right)$ are influenced by the matrices $G_1$ and $G_2$, which incorporate control inputs from ${\tilde{\mathit{B}}}_T$, ${\tilde{\mathit{B}}}_D$, and ${\tilde{\mathit{B}}}_M$. Significant differences in the control matrices can be amplified through quadratic operations in Equation (32), leading to wide disparities in the magnitudes of the elements in $K\left(t\right)$ and exacerbating its ill-conditioning.

In practical terms, control inputs such as fin deflection $u_{T/D}$ and commanded accelerations $u_A$ differ fundamentally in scale and system impact. For example, a single-degree change in fin deflection might equate to a 20 $\mathrm{m}/{\mathrm{s}}^2$ change in acceleration. The difference in control efficiency affects how sensitive the system state is to increments in $u_{T/D}$ and $u_M$.

In multiple-input–multiple-output (MIMO) systems, the disparity in control effectiveness across different actuators can be quantitatively expressed through the singular values of their control matrices. The largest singular value of control matrices offers a measure of the “control efficiency”, reflecting the sensitivity of the system’s state to changes in respective control inputs. This is because the largest singular value in a control matrix provides critical insight into the maximum stretch factor that the matrix imparts to any vector it multiplies.

To illustrate the relationship between matrix scaling and the largest singular values, we consider a scenario where the short-period dynamics parameters for the target and defender are taken from [45] and [17], respectively. We define the scaled control matrices as follows:

$$\begin{array}{cc}\hfill {\tilde{\mathit{B}}}_T^{*}& ={\mathit{B}}^{*}·{\tilde{\mathit{B}}}_T\hfill \\ \hfill {\tilde{\mathit{B}}}_D^{*}& ={\mathit{B}}^{*}·{\tilde{\mathit{B}}}_D\hfill \\ \hfill {\tilde{\mathit{B}}}_A^{*}& ={\mathit{B}}^{*}·{\tilde{\mathit{B}}}_A\hfill \end{array}$$

(48)

where ${\mathit{B}}^{*}$ is a scaling factor. The sensitivity of the system’s state to different control inputs is demonstrated by comparing the largest singular values of the control matrices under various scaling factors.

Figure 2 presents the variation of the largest singular value $\overline{\sigma }$ for each control matrix under different scaling factors. Given that the defender’s control input typically exhibits the highest control efficiency, we use the largest singular value of ${\tilde{\mathit{B}}}_D$ as the benchmark. The graph shows that the largest singular values approximate a linear relationship on a logarithmic scale. For the flight characteristics considered here, when the scaling factor for ${\tilde{\mathit{B}}}_T$ is about 10 and for ${\tilde{\mathit{B}}}_A$ about 250 the relationships satisfy

$$\overline{\sigma }\left({\tilde{\mathit{B}}}_T^{*}\right)\approx \overline{\sigma }\left({\tilde{\mathit{B}}}_D\right)\approx \overline{\sigma }\left({\tilde{\mathit{B}}}_A^{*}\right).$$

(49)

This finding suggests that the magnitude differences in control matrices must be considered when applying the IGC active protection strategy, as they directly affect the ill-conditioning of $K\left(t\right)$ and potentially the robustness of the control solution.

Furthermore, in order to check the accuracy of the calculation for $K^{-1}\left(t\right)$ during the simulation, we introduce the condition number as a real-time monitor. The condition number of a matrix K is defined as

$$\mathrm{cond}\left(K\right)=\parallel K\parallel \parallel K^{-1}\parallel ,$$

(50)

where $\parallel K\parallel$ is a norm of K. The condition number quantifies the sensitivity of the solution of a system of linear equations to changes in the input or errors in the data. A low condition number indicates that the matrix is well-conditioned, meaning that it is stable and reliable for numerical calculations. Conversely, a high condition number suggests that the matrix is ill-conditioned, which can lead to significant numerical errors in computations, especially inverses. Matrices with very high condition numbers may result in large errors in computed solutions, indicating that even small errors in data can lead to unreliable or incorrect results.

In Section 5, we demonstrate through specific examples how the methodologies proposed in this paper effectively reduce the condition number of matrix $K\left(t\right)$. By improving the conditioning of $K\left(t\right)$, we enhance the reliability and accuracy of the matrix inversion operations critical to the IGC active protection strategies.

4.3. Simplified IGC Active Protection Strategy

In the development of the IGC active protection strategy, while modeling both the target and the defender within the IGC framework is feasible, it tends to increase computational complexity and numerical instability. A more streamlined approach can be adopted by focusing the IGC model primarily on the defender.

In general, the defender exhibits more significant agility and dynamic response than the target in active protection scenarios. By concentrating the IGC model on the defender, we can effectively capture and manage these dynamics while reducing the model’s overall complexity. This simplified strategy simplifies calculations while addressing the more pronounced coupling effects between guidance and control loops in high-speed and highly maneuverable defenders.

The state vector for this simplified model is chosen as

$$\mathit{x}={\left[y_{AT},{\dot{y}}_{AT},a_T,y_{AD},{\dot{y}}_{AD},{\alpha }_D,q_D,a_A\right]}^{\top }.$$

(51)

This state vector includes the relative positions, velocities, and accelerations of the attacker and target as well as the angle of attack and pitch rate relevant to the defender.

The linearized dynamics model representing the three-player conflict can be expressed as

$$\dot{\overline{\mathit{x}}}=\overline{\mathit{A}}\mathit{x}+{\overline{\mathit{B}}}_T{u}_T+{\overline{\mathit{B}}}_D{u}_D+{\overline{\mathit{B}}}_A{u}_A,$$

(52)

where

$$\overline{\mathit{A}}=\left[\begin{array}{ccc}\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -1/{\tau }_T\end{array}& {\left[\mathbf{0}\right]}_{3\times 4}& \begin{array}{c}0\\ -1\\ 0\end{array}\\ {\left[\mathbf{0}\right]}_{4\times 3}& \begin{array}{ccc}0& 1& {\left[\mathbf{0}\right]}_{1\times 2}\\ 0& 0& {\mathit{C}}_D\\ {\left[\mathbf{0}\right]}_{2\times 1}& {\left[\mathbf{0}\right]}_{2\times 1}& {\mathit{A}}_D\end{array}& \begin{array}{c}0\\ 1\\ 0\end{array}\\ {\left[\mathbf{0}\right]}_{1\times 3}& {\left[\mathbf{0}\right]}_{1\times 4}& -1/{\tau }_M\end{array}\right]$$

(53)

$${\overline{\mathit{B}}}_T={\left[\begin{array}{ccc}\begin{array}{ccc}0& 0& 1/{\tau }_T\end{array}& {\left[\mathbf{0}\right]}_{1\times 4}& 0\end{array}\right]}^{\top }$$

(54)

$${\overline{\mathit{B}}}_D={\left[\begin{array}{ccc}{\left[\mathbf{0}\right]}_{1\times 3}& \begin{array}{ccc}0& -d_D& {\mathit{B}}_D^{\top }\end{array}& 0\end{array}\right]}^{\top }$$

(55)

$${\overline{\mathit{B}}}_A={\left[\begin{array}{ccc}{\left[\mathbf{0}\right]}_{1\times 3}& {\left[\mathbf{0}\right]}_{1\times 4}& 1/{\tau }_M\end{array}\right]}^{\top }.$$

(56)

This model simplifies the integration of guidance and control for the defender by isolating and directly addressing the most dynamically significant elements of the system. The derivation follows procedures similar to those previously discussed; for the sake of brevity, it is not repeated here.

5. Simulation Results and Analysis

This section describes the simulation setup and presents the results of two specific simulation examples to validate the proposed IGC approach within the context of aerial vehicle active protection scenarios. The results are compared to traditional methods where guidance and control loops are designed separately.

5.1. Simulation Setup

The simulation is designed to evaluate the effectiveness of the IGC strategy proposed in this paper. The defender model is based on the one described in [17]. The fin deflection commands generated by the autopilot (AP) or IGC are transmitted to the actuator. In aerial vehicle systems, the actuator is a critical component that converts the autopilot’s commands into physical movements, such as the deflection of tail fins, to achieve the desired motion. Considering the need for rapid response, actuators in aerial vehicles are typically high-bandwidth devices which exceed the bandwidth of the flight control loop itself. Although most modern actuators are electromechanical, hydraulic options are used in specific applications.

For the initial design and analysis phases, it is common to model the actuator dynamics using a simplified second-order transfer function. While this model does not capture the full complexity of the actuator hardware, it provides a reasonable approximation for understanding how the actuator’s response characteristics can influence the performance of the guidance and control loops. The second-order modeling is expressed by the transfer function

$$\frac{\delta \left(s\right)}{{\delta }_c\left(s\right)}=\frac{{\omega }_a^2}{s^2+2{\zeta }_a{\omega }_{a}s+{\omega }_a^2},$$

(57)

where $\delta \left(s\right)$ is the actuator output, ${\delta }_c\left(s\right)$ is the command from the autopilot, ${\omega }_a$ is the natural frequency, and ${\zeta }_a$ is the damping ratio of the actuator. Reasonably designing the parameters of actuators is also crucial for the stability of a control system. A series of curves comparing different settings of the second-order response and the characteristics of these responses are shown in Figure 3.

For subsequent simulations, the parameters chosen for the actuator model are ${\omega }_a=150$ rad/s and ${\zeta }_a=0.7$. These parameters are selected to balance responsiveness and stability, ensuring that the actuator can adequately follow the control commands without inducing undue oscillations or delays in the vehicle’s response.

5.2. Simulation Scenario I

In this simulation example, the attacker directly employs the Proportional Navigation (PN) guidance law to engage the target. The defender’s performance is evaluated under two different configurations: an IGC design, and a traditional approach in which guidance and control loops are designed separately, i.e., PN+AP. In the traditional setup, the defender uses a PN guidance law coupled with a classic three-loop autopilot to track acceleration commands. The initial conditions for this simulation are detailed in

Table 1. Initial parameters for simulation scenario I.

Parameter	Value
$(x_{T0},y_{T0})$	(0, 1000) m
$(x_{D0},y_{D0})$	(0, 1000) m
$(x_{A0},y_{A0})$	(2000, 1200) m
${\tau }_A$	0.5 s
$u_A^{max}$	40 g

:

The engagement trajectories of the vehicles are depicted in Figure 4. Under the IGC design, the defender successfully intercepts the attacker at 2.55 s with a miss distance of only 0.16 m. This high-precision interception effectively neutralizes the threat, preventing any damage to the target. In stark contrast, the defender utilizing the traditional AP+PN design results in a significantly larger miss distance of 28.67 m. This considerable deviation results in the attacker successfully hitting the target, demonstrating the limitations of the traditional approach in high-stakes scenarios requiring precise interception.

In addition, the defender employing the IGC design follows a more direct and efficient flight path. This path not only demonstrates the system’s enhanced responsiveness to dynamic threats, it indicates a lower demand for the defender’s maneuvering capabilities; conversely, the trajectory of the defender using the traditional PN+AP design exhibits more curvature and requires greater maneuvering effort, which can lead to higher fuel consumption and greater stress on the defender.

The ZEM is a critical metric in three-player aerial engagements, providing insight into the effectiveness of the active protection strategy by measuring the predicted miss distance (assuming that no further maneuvers are performed). An analysis of the ZEM for both the IGC strategy and the PN+AP reveals significant differences in performance and stability.

For the defender utilizing the IGC design, the ZEM $z_2$ exhibits a stable and smooth convergence towards zero, as shown in Figure 5. This stable convergence indicates that the defender effectively locks the attacker within the interception triangle, maintaining a consistent trajectory that ensures successful interception. Conversely, the ZEM values for the defender employing the traditional PN+AP display considerable fluctuations, highlighting instability and a lack of control precision. The erratic nature of the ZEM curve when using PN+AP suggests that the defender frequently adjusts their path, reacting to the attacker’s maneuvers rather than anticipating, which leads to a failure to intercept.

In the proposed IGC design, computing the inverse of matrix $K\left(t\right)$ is crucial for solving the Riccati differential equation in order to derive the solution matrix $P\left(t\right)$. The condition number of matrix $K\left(t\right)$ serves as a critical indicator of its invertibility and the potential for numerical instability. A low condition number suggests that the matrix is not sensitive to small perturbations or numerical errors, ensuring the reliability of the matrix inversion operation.

Throughout the engagement process, the condition number of matrix $K\left(t\right)$ maintained values typically below 10, as depicted in Figure 6. This indicates that matrix $K\left(t\right)$ remained non-singular and well-conditioned, ensuring the accuracy and reliability of the inversion computations required for solving $P\left(t\right)$. This stability in the condition number affirms the method’s robustness against potential numerical issues, particularly when real-time onboard calculations are required for dynamic aerial engagements.

Our approach based on the differential game framework does not necessitate prior knowledge of the adversary’s exact strategies or their functional form; however, it does demand specific battlefield information, such as the attacker’s time constant ${\tau }_A$ and acceleration $a_A$. Although ${\tau }_A$ can be inferred from open-source information or engineering experience, and $a_A$ can be estimated in real-time using filtering algorithms, these methods introduce inevitable errors. In order to assess potential disadvantages, we conducted Monte Carlo simulations considering deviation in estimates of ${\tau }_A$ and $a_A$. In addition, we introduced varying levels of noise interference in the measurement of the Line-of-Sight (LOS) rotation rate. The defender’s kill radius was set to 5 m, meaning that an interception is successful if the miss distance is less than 5 m. The attacker’s initial position was varied between 1500–3000 m laterally and 200–1800 m longitudinally.

We compiled the defender’s interception success rates under different noise levels, as shown in

Table 2. Success rate under noise and deviation.

Noise and Deviation	0	10%	15%	25%
Success Rate	100%	100%	100%	100%

. The results indicate that when the measurement noise and model deviation are below 25% the proposed IGC strategy’s performance remains largely unaffected, demonstrating the method’s robustness to noise interference.

These findings address concerns about sensor noise while providing a more comprehensive evaluation of our IGC strategy’s performance under realistic conditions. In addition, they demonstrate the practical applicability of our approach in scenarios where sensor imperfections and countermeasures are present.

5.3. Simulation Scenario II

The vehicles’ simulation parameters are detailed in

Table 3. Initial parameters of the vehicles.

Initial Parameters	Value
$(x_{T0},y_{T0})$	(0, 1000) m
$(x_{D0},y_{D0})$	(0, 1000) m
$(x_{A0},y_{A0})$	(3500, 800) m
${\tau }_A$	0.5 s
$u_A^{max}$	20 g

. Faced with the threat of interception, the attacker employs weave maneuvers to evade the defender. If the attacker succeeds in evading the defender, it then switches to the proportional navigation (PN) guidance law to strike the target.

With the implementation of the IGC design, the trajectories of the three vehicles are illustrated in Figure 7a. The defender successfully intercepts the attacker with a miss distance of 0.33 m, demonstrating the effectiveness of the IGC approach in protecting the target.

To highlight the superiority of the proposed IGC approach, we compare it against a traditional scheme where the guidance and control loops are designed separately. In this traditional configuration, the defender employs a PN guidance law coupled with a three-loop autopilot to follow the acceleration commands. The trajectory for this separation design scheme is depicted in Figure 7b. In contrast to the IGC approach, the defender using the PN+AP scheme fails to intercept the attacker, resulting in a significant miss distance of 49.64 m. Consequently, this failure allows the attacker to successfully hit the target.

The ZEM $z_2$ comparison is illustrated in Figure 8. The defender utilizing IGC shows markedly less influence from the attacker’s evasive maneuvers, successfully locking the attacker within the interception triangle. In contrast, the defender using the separate PN+AP approach experiences significant fluctuations in $z_2$, ultimately failing to intercept the attacker.

The acceleration curves of the flight vehicles under both the integrated and separated guidance and control schemes are presented in Figure 9a,b. The simulation with the IGC design ends at 4.44 s, whereas the simulation with the PN+AP approach extends to 5.94 s. Notably, the attacker’s acceleration profile before evading the defender remains consistent in both figures until it switches to a PN guidance law heading to the target.

The acceleration of the defender under these two scenarios is particularly noteworthy. With the IGC approach, the defender experiences brief initial oscillations before quickly stabilizing. Throughout the engagement, the defender maintains relatively low overload levels generally not exceeding 5 g, with a maximum peak below 10 g. In contrast, under the traditional separated design the defender’s demand for maneuverability significantly increases, with peak overloads reaching up to 30 g.

The comparison in acceleration clearly illustrates the advantages of integrating guidance and control loops in the IGC design. By considering the dynamic characteristics of the vehicle comprehensively, the IGC approach effectively reduces the demand for maneuverability in high-dynamic environments, leading to more stable and efficient flight performance. Such an integrated design enhances interception efficacy while minimizing the stress on the vehicle’s structural and operational limits.

The fin deflection of the defender is depicted in Figure 10, showing a smooth variation that accounts for actuator response delays. This controlled adjustment ensures compatibility with realistic hardware capabilities, demonstrating the practical applicability of the IGC approach.

To ensure the numerical accuracy during the inversion operation of matrix $K\left(t\right)$, we analyzed the changes in the elements of matrix $K\left(t\right)$ and its condition number. As displayed in Figure 11, the variation curves of each element confirm that the matrix does not exhibit any singularity, validating the feasibility of obtaining matrix $P\left(t\right)$ through inversion. In addition, we compared the theoretical solution obtained from the linearized model with the actual values observed in the nonlinear model simulations. The close agreement between the two sets of results shown in Figure 11 serves as empirical evidence supporting the validity and feasibility of the linearization method employed in this work.

Further, as shown in Figure 12, the condition number of matrix $K\left(t\right)$ generally remains below 10 and does not exceed 100. This indicates that matrix $K\left(t\right)$ is not ill-conditioned, ensuring that the precision of the inversion operation is safeguarded.

6. Conclusions

This paper has proposed a cooperative IGC strategy for active target protection in a three-player conflict scenario. The IGC model was established by linearizing both the translational and rotational dynamics of the target and defender vehicles. Differential game theory and the sweep method were applied to derive and solve the Riccati differential equation analytically, providing an optimal control solution. The theoretical optimality of the derived strategy was proven and the computational accuracy of the solution was analyzed in terms of the control matrix singular values and solution matrix condition numbers. Compared to traditional separated guidance and control designs, the IGC strategy demonstrates superior performance, smoother control trajectories, and reduced maneuverability demands on the defender, leading to enhanced interception efficacy and minimized stress on the vehicle’s structural and operational limits. The proposed IGC strategy provides a promising framework for active target protection in aerial vehicle engagements. Future work could explore extending the model to incorporate more complex vehicle dynamics, multi-vehicle cooperation, and robustness against uncertainties in system parameters.