Multiple Tactical Missiles Cooperative Attack With Formation-Containment Tracking Requirement Along the Planned Trajectory

This paper studies the cooperative guidance and control problems with formation-containment tracking requirement under directed topologies for multiple tactical missiles. The states of the formation-leaders can not only keep a parallel triangle formation but also make tracks for the state trajectory generated by the tracking-leader, while the states of followers converge to the convex hull formed by those of the formation-leaders. Firstly, an integrated cooperative guidance and control framework is proposed, and combat missions are distinguished in line with tactical missiles rewarded different functions in practical applications. Then, a terminal slide mode guidance law with impact angle is presented in three-dimensional space, which ensures attack missile to attack the target in finite time. Sufficient conditions for multiple tactical missiles to achieve formation-containment tracking are derived. On the basis of Lyapunov theory, it is shown that the expected formation-containment tracking can be realized by reconnaissance and decoy missiles in the presence of the tracking-leaders’ unknown control input. In addition, in the case of unknown drag force, the six degrees of freedom missile controller is designed by combining genetic algorithm and disturbance observer to dynamically adjust the rudder deflection and thrust, ensuring stable tracking of overload command during cooperative attack. Finally, numerical simulations validate the feasibility and effectiveness of the proposed results.


I. INTRODUCTION
Recently salvo attack of tactical missiles has been regarded as an effective approach to penetrate anti-air defense systems. Researches on cooperative guidance and control under the background of multiple missiles combat application have been hot and attractive issues.
To achieve salvo attack, we can perform cooperative guidance and control of tactical missiles through two different strategies, i. e., individual homing and cooperative homing. For individual homing, a suitable common impact time is preprogrammed manually beforehand, and then the open-loop or closed-loop guidance and control laws are adopted for each tactical missile individually to achieve the designated The associate editor coordinating the review of this manuscript and approving it for publication was Xiwang Dong. impact time independently. At present, the scheme used in this kind of guidance law design can be divided into two main categories. One is based on the traditional proportional guidance law. In [1], a feedback term of time-to-go error was added to the guidance law, and an impact-time-control guidance (ITCG) law was proposed. [2] designed a proportional navigation with a variable coefficient, which was based on the relationship between the remaining flight time and the effective navigation ratio. Another category is based on the modern control theory. In addition to impact time, the coordinated impact angle needs to be considered in order to make the missile warhead obtain better killing effects. [3] linearized the model under the small heading error assumption, the impact time and impact angle can be controlled simultaneously using the guidance law based on the optimal control theory. In [4], the guidance law based on the LQR theory was presented to VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ attack the uniform motion target with expected impact angle and minimum miss-distance. Unlike individual homing, a hierarchical framework free from predesigned impact time is built for cooperative homing via online data links and communication topology. In [5], missiles were classified into many groups, where each group belonged to the centralized leader-follower framework, and the leaders of different groups communicated with each other through the nearest-neighbor topology. Finite-Time Cooperative Guidance (FTCG) laws considering the saturation constraint on field-of-view (FOV) were firstly addressed to accommodate the communication topology in a single group, and then an improved sequential approach was proposed to FTCG-FOV to accommodate the communication topology between groups. [6] investigated the uncertainty of communication composed of stochastic network and additive noise, analysed the time-to-go error of each missile in the communication network, and set a new impact time through the data link using the mean square error.
Although the aforementioned cooperative homing to realize salvo attack has been widely investigated, combat missions are not distinguished according to various characteristics of tactical missiles in practical applications, and the formation-containment tracking issues are not considered. To improve the penetration effectiveness under the condition of system-of-systems combat, combat scenarios are proposed as follows: (i) Attack missile with warhead should penetrate the enemy's anti-air defense systems to destroy a target either moving at a relatively low speed or non-maneuvering, and more space on the missile can be taken up to load explosives. (ii) Reconnaissance missiles mounting detecting equipment are required to realize an expected formation pattern while tracking the designated trajectory generated by attack missile simultaneously, which can promote the detection precision in the way of detecting the target cooperatively. (iii) Decoy missiles will converge to the convex hull spanned by reconnaissance missiles synchronously, and interfere the enemy defenses by the loaded interference power amplifiers. In addition, the safety of attack missile can be guaranteed due to the fact that reconnaissance and decoy missiles will disperse the fire power of the enemy. Therefore, the proposed methods outperform the conventional multi-missile combat mode. As far as the author knows, there are few literatures on the formation-containment tracking for multiple tactical missiles. With the rapid development of consensus theory (see [7]- [18] and references therein), more studies related with consensus strategies are extended to solve the formation or containment problems of the multi-agent system. [7] presented formation control problems for the multi-agent system utilizing consensus-based approaches, where a given geometry was approached by the states or outputs of multiple agents, and provided proofs that some traditional formation methods such as virtual structure, behavior-based ones, and leader-follower can be converted into the consensus based ones. Time-varying formation protocols for the multi-agent system with general linear dynamics and switching directed interaction topologies were dealt with in [8]. Besides achieving the time-varying formation, the entire multi-agent system can also make tracks for the desired trajectory through formation tracking protocols in [9], [10]. In [11], containment control problems were investigated for general linear multi-agent systems with timevarying delays, where the states or outputs of the followers were required to converge to the convex hull spanned by those of leaders. Recently, due to the increasing researches on formation control and containment control, the formationcontainment problem has been frequently considered, which requires that the states or outputs of leaders achieve desired formation while the states or outputs of followers converge to the convex hull formed by those of leaders. [12], [13] addressed formation-containment analysis and design problems for time-delayed high-order linear time-invariant multi-agent systems with directed interaction topologies or heterogeneous architecture. Formation-containment control problems for multi-UAV systems with directed topologies were presented in [14], and a formation-containment platform with multiple UAVs was demonstrated.
To the best of our knowledge, even though each of the four aspects listed above including formation, containment, formation-tracking and formation-containment has been taken into account in the existing literatures, formationcontainment tracking problems considering all of the four aspects remain to be further investigated.
In this paper, formation-containment tracking problems for multiple tactical missiles subjected to directed topologies are proposed. Firstly, an integrated cooperative guidance and control framework is presented, in which the guidance layer, coordination layer, and control layer are assigned their own engineering problems by guidance engagement geometry, kinematic dynamics and actuator control loop, respectively. At the theoretical level of formation-containment tracking, tactical missiles rewarded different functions act as the tracking-leader, formation-leader and follower. Then a terminal slide mode variable structure guidance law with impact angle constraint is designed for the tracking-leader, which guarantees attack missile hit the target under the proposed guidance law in finite time. Moreover, formationcontainment tracking protocols are derived utilizing neighboring information of different agents. It is proved that formation-leaders can maintain a time-varying formation with target enclosing achieved, meanwhile the followers are supposed to converge to the convex hull specified by the formation-leaders. Genetic algorithm and disturbance observer (DOB) are employed to design the control systems for salvo attack of different missiles, which allow all six degrees of freedom (DOF) of the missile to be controlled using four control input through thrust and rudder deflection. Finally, a formation-containment tracking anti-ship scenario is introduced, where 6 tactical missiles are considered to perform a cooperative attack against a stationary target. The numerical simulation results are addressed to illustrate the performance of the obtained results.
Compared with the relevant results on formationcontainment tracking control, the main contributions of the current paper are threefold. Firstly, an integrated cooperative guidance and control framework is proposed. Multiple tactical missiles can not only realize cooperative attack autonomously but also constitute a predefined time-varying formation-containment geometry while tracking the trajectory of attack missile. In [2], [5], [6], [19], [20], [23], only time or space cooperative guidance problems were considered and there existed no formation or containment process. Secondly, the tracking-leader is free from certain input. Attack missile adopts the sliding mode control method to construct the homing guidance law, which acts as an unknown control input for the tracking-leader on the formation-containment tracking control level. In [10], the formation-leader performed a formation tracking task without the tracking-leader control input. The protocols to cope with cases of a tracking-leader with naught or available input cannot be extended to the non-cooperative tracking-leader case. Thirdly, the dynamic characteristics of the tactical missile in control loop are taken into account in parallel to the consideration of formation-containment tracking in guidance loop. Genetic algorithm can quickly adjust the parameters of rudder control loop, which enables the tactical missile to track the specific overload command required to follow the desired cooperative trajectory smoothly and quickly. A DOB in axial overload control loop is provided with no prior knowledge of the drag force. In [7]- [9], [11]- [13], [15]- [17], the dynamics for each agent were restricted to swarm systems without concrete objects, and the practical application of formation or containment approaches for missiles were not investigated. In [20], [21], [23], the missiles were assumed to be mass points. Because the dynamic characteristics will affect the reaction speed of the tactical missile to the cooperative guidance command dramatically, the results for aircrafts free from the dynamic characteristics cannot be directly applied to salvo attack considering missiles' dynamic characteristics and the cooperative guidance requirement simultaneously. Only 3 DOF cooperative problems were considered in [22]. Although 6 DOF cooperative problems were investigated in [24], [25], explicit knowledge of the drag coefficient was required in [24], while controller involving turbojet engine thrust was not considered in [25] with respect to axial overload control. It should be noted that [26] only focused on the mid-course guidance of cooperative attack of multiple missiles formation, and cooperative engagement process involving dynamic characteristics was ignored.
The remainder of this paper is organized as follows. Section II introduces some useful results and basic concepts and gives the problem formulation of the formation-containment tracking cooperative attack for multiple tactical missiles. In Section III, the integrated cooperative guidance and control framework is proposed, and the approach to implement the guidance law during the engagement, design the formation-containment tracking protocols and conduct the control of each missile through the rudder/axial overload control loop are presented. Numerical simulations are discussed in detail in Section IV. At last, Section V concludes this paper.

A. BASIC CONCEPTS ON GRAPH THEORY
A directed graph with M + N + 1 nodes is described by represents the adjacency matrix with nonnegative weights w ij . An edge of graph G is defined as e ij = (v i , v j ). The weight w ij > 0 if and only if e ij ∈ E, and w ii = 0(i = 0, 1, · · · , M + N ). Denote N i = {v j ∈ V : (v j , v i ) ∈ E} by the set of neighbors of node v i . A sequence of ordered edges (v i k , v i k+1 )(k = 1, 2, · · · , g − 1) is said to be a directed path The multi-agent system is made up of 1 tracking-leader, N formation-leaders and M followers. If an agent has no neighbors, it is called a tracking-leader. The formation-leader can only receive information from other formation-leaders or the tracking-leader. The neighbors of a follower cannot be anything but a cluster of other followers or formation-leaders. For notational brevity, denote by E = {1, 2, · · · , N } and F = {1 + N , 2 + N , · · · , M + N }. In this paper, consider a multiple tactical missile system with M + N + 1 missiles, and each is denoted as an agent within the multi-agent system.

B. GUIDANCE AND CONTROL GEOMETRY
The three-dimensional (3D) homing guidance geometry is depicted in Fig.1, where OXYZ , OX 4 Y 4 Z 4 are the inertial coordinate system and the line-of-sight (LOS) coordinate system, respectively. R is the range-to-go between the missile and the target along the LOS direction, q ε represents the LOS angle in the pitching direction, q β is the LOS angle in the yawing direction. One obtains The coordinate transformation from the LOS coordinate system to the inertial system can be described as follows cos q ε cos q β sin q ε − cos q ε sin q β − sin q ε cos q β cos q ε sin q ε sin q β sin q β 0 cos q β   (3) The control system design geometry is given in Fig.2, where OX 2 Y 2 Z 2 is the ballistic coordinate system. The flight-path angle and speed bank angle are described by θ and ψ v . Like the case of transformation from OX 4 Y 4 Z 4 to OXYZ , the coordinated transformation matrix from the inertial system to the ballistic coordinate system is given by Let the projection of the absolute acceleration and the relative acceleration vector between missile and target in Then the relative movement between missile and target can be written as follows The latter two formulas of Eq. (5) can be transformed into Let state variables be where q εd and q βd are the expected LOS angle within 3D space. Therefore, the equations of missile-target engagement take the form In this paper, taking the target as stationary or shifting at a relatively low speed, then let a T ε = a T β = 0, and further define Eq. (7) is rewritten as Coordinate transformations Eq. (3) and (4) are further adopted to transfer the control input from the LOS coordinate system to the ballistic coordinate system, and assuming the overload in X 2 -direction for the formation-leader is steered to be zero, one obtains that where D = [n yc , n zc ] T stands for overload command along Y 2 -direction and Z 2 -direction in OX 2 Y 2 Z 2 and (10), as shown at the bottom of the next page.

1) DYNAMIC MODELING
Based on the relationship between missile motion state and overload, one can obtain where n x , n y , n z T is overload of missile in the ballistic coordinate system. The speed, flight-path angle and speed bank angle are described by v, θ and ψ v . g is gravity acceleration. According to dynamical equations, overload can be written as where T is engine thrust. Drag force X , lift force Y , lateral force Z are the components of aerodynamic force. α, β, γ v represent attack angle, sideslip angle and velocity slope angle. We can get these three angles based on the following geometric equation where γ is the roll angle, ϕ is the pitch angle, ψ is the yaw angle. According to the transformation relationship between inertial coordinate system and body coordinate system [27], one can obtain that  φ ψ γ The angular velocities are described in the body fixed frame attached to the center of mass, and the rotational motion equations are written as where m is mass of missile, J ∈ R 3×3 is the moment of inertia, M x , M y , M z T is the force and aerodynamic torque applied to the missile. Ignoring the effect of damping torque, the approximations of the forces and moments are given by where q = 1 2 ρv 2 is dynamic pressure, and ρ stands for air density. S, L, δ x , δ y and δ z denote reference area, reference length, aileron angular deflection, rudder angular deflection and elevator angular deflection.
x , m α z , m δz z are the aerodynamic moment coefficients.

C. PROBLEM DESCRIPTION
For a multiple tactical missile system with M + N + 1 tactical missiles, the communication topology of the multi tactical missile system can be represented by a directed graph G with node i denoting the tactical missile i and edge e ij denoting the communication between tactical missile i and j. The communication strength from tactical missile j and tactical missile i is assumed to be w ij .
In the cooperative control protocol level, the second-order system is introduced to describe the kinematic dynamics of the tracking-leader where x 0 (t) ∈ R n , v 0 (t) ∈ R n , and u 0 (t) ∈ R n is the position, velocity, and acceleration vector of the trackingleader described in OXYZ . The nonzero acceleration u 0 (t) is considered to be unknown to all the other agents of the multiple tactical missile system, which is to be determined in Section III.B. The dynamics of the N formation-leaders and M followers can be given by where x i (t) ∈ R n , v i (t) ∈ R n , and u i (t) ∈ R n denote the position, velocity, and acceleration vector of the ith tactical missile defined in OXYZ , respectively. u i (t)(i = 1, 2, · · · , M + N ) are the designated protocols (see the analysis and design of the protocols in Section III.C) that should be approached by formation-leaders and followers. The relationship between the acceleration and overload is determined by where n i (t) is the overload vector for agent i.
In this study, a 3D many-to-one engagement scenario with formation-containment tracking requirement is considered, and the following assumptions are adopted.
Assumption 1: The nonzero overload n i (t) is bounded and there exists a maximum magnitude of the overload n m such that n i (t) ≤ n m due to the finite actuation capability for physical actuators during flight.
Assumption 2: Considering that long-term motion characteristics are mainly determined by long period mode in flight, short period process and the effect of rotation in the control process is neglected in the guidance realization.
Assumption 3: The directed graph among the trackingleader and the formation-leaders in G contain a spanning tree, in which the tracking-leader is the root node. There is at least one directed path between each follower and some formationleader.
Under Assumption 3, the Laplacian matrix L ∈ R (N +M +1)×(N +M +1) associated with the directed graph G has the following form Therefore Eq. (17) and (18) can be rewritten aṡ Definition 1: The formation-leader i is said to achieve the expected time-varying formation tracking specified by where h i (t) ∈ R 3 (i = 1, 2, · · · , N ) characterize the desired time-varying formation configuration, and h i (t) = [h ix (t), h iv (t)] T denotes piecewise continuously differentiable formation vector for the formation-leader i. Definition 2: The follower j is said to achieve containment if where λ jk (j ∈ F; k ∈ E) represent nonnegative constants satisfying N k=1 λ jk = 1 for the follower j. Definition 3: The multiple tactical missile system Eq. (21) is said to achieve formation-containment tracking if Eq. (22) and (23) hold for all the formation-leaders and followers, respectively.
Not that all the results hereafter can be easily extended to the higher dimensional case with n ≥ 2 using the Kronecker product. In the three-dimensional space, n is set to be 3.
Lemma 1 [28]: Consider the following nonlinear systeṁ where ϒ (0) = 0, ξ ∈ R n , U 0 ∈ R n is defined as an open-loop field at the origin of the system. Given a positive definite Then the origin of the system is a finite-time stable equilibrium, and the pause time t r depends on the initial value ξ (0) = ξ 0 , its upper bound is Lemma 2 [16]: Under Assumption 3, there exists a diago- The formation-containment tracking issues for multiple tactical missiles considered in this current paper are mainly focused as follows: (1) how to construct a hierarchical framework for cooperative attack with formation-containment tracking requirement; (2) how to implement a terminal guidance with a desired terminal impact angle in finite time against stationary target within 3D space; (3) under what conditions formation-containment tracking can be accomplished with the tracking-leader of uncertain control input; (4) how to develop the formation-containment tracking control protocol and 6-DOF missile control system associated with nonlinear uncertainties.

A. INTEGRATED COOPERATIVE GUDIANCE AND CONTROL FRAMEWORK
For the proposed integrated cooperative guidance and control framework, different tactical missiles are firstly divided into attack missiles, reconnaissance missiles and decoy missiles according to their respective combat missions (see Table 1). Without loss of generality, a schematic diagram of cooperative attack with formation-containment tracking requirement is shown in Fig.3, in which multiple tactical missiles communicate with each other through data links, and achieve cooperative guidance and control along the attack direction using the nearest neighbors' information.
Instead of completing the combat mission independently by high-performance multi-purpose missiles or platforms, the capabilities are scattered to multiple missiles, which jointly form a combat system to complete the mission. The battlefield orientation of the formation-containment tracking framework is that the attack missile converges to a given impact angle in the course of engagement, and strikes the heavily defensed targets in finite time to guarantee better lethal performances of the warheads. Reconnaissance missiles with detection equipment are expected to generate specific formation configurations, which expands the  cooperative search radius of the target. Then the detected target location is broadcast to attack missile that will adjust its fight online to realize fixed-point attack. Decoy missiles suppress and disperse the fire power of the enemy in the confrontation area. Consequently, the missile group will survive the threat and penetrate the enemy's anti-air defense systems successfully, forming an overwhelming advantage against the surrounding combat opponents.
Then the flowchart of the proposed cooperative framework listed in Fig.4 is implemented through a step-by-step layer such as the guidance layer, coordination layer and control layer. Missile and target guidance engagement geometry, kinematic dynamics and actuator control loop are applied to depict different layers.
Remark 1: Overloads in diverse directions are coordinated variables in the proposed framework, which link different layers together. In the course of attack missile engagement, the overload design in Section III.B has been transformed from the LOS coordinate system (OX 4 Y 4 Z 4 ) to the ballistic coordinate system (OX 2 Y 2 Z 2 ) in Section II.C during the engagement process of attack missile. The cooperative scheme in Section III.C is implemented in the inertial coordinate system (OXYZ ). Then the ballistic coordinate system is introduced into the three-loop missile autopilot and thrust adjustment design procedure in Section III.D. Accordingly, coordinate transformations from the ballistic coordinate sys- tem to the inertial system, and on the contrary, from the inertial system to the ballistic coordinate system are implemented once, respectively. For the conciseness and clarity, the subscript of different coordinate systems will be omitted throughout this paper.

B. TERMINAL SLIDE MODE VARIABLE STRUCTURE GUIDANCE LAW WITH IMPACT ANGLE
Considering the 3D guidance system (9), define the sliding mode surface as follows where is the longitudinal and lateral motion sliding mode surface, respectively. β 1 > 0 and β 2 > 0 are any positive constants. Positive odd numbers p 1 > 0, q 1 > 0, p 2 > 0 and q 2 > 0 are chosen to satisfy 1 < p i q i < 2. Design the terminal guidance law as where The stability analysis of the guidance loop is shown in Appendix.
Then we can get the following theorem. Theorem 1: In the process of terminal guidance for attack missile, the guidance law (28) can ensure that the 3D guidance system (9) converges to the sliding mode surface s = 0 in finite time and the system converges to zero on the sliding mode surface in finite time.
Remark 2: The guidance law based on sliding mode control method is applied widely in the existing literatures, and the LOS angle and angle rate is usually chosen as the sliding mode surface (see [29][30][31][32][33][34][35][36]). Theoretically speaking, the smaller the parameters β i (i = 1, 2) and q i /p i (i = 1, 2) are, the shorter the time of system (9) converges to the sliding mode surface s = 0, which indicates a better control performance. However, when the state of system has reached the sliding mode surface s = 0, the smaller the parameters β i (i = 1, 2) are, the longer the time of system (9) converges to the equilibrium point, which may result in an unacceptable performance or even destabilize the system. Therefore, β i (i = 1, 2) should select a moderate value to trade off between control performances in convergence to and during the sliding mode surface.
Remark 3: Attack missile maintains an axial overload of zero through the tangential balance between thrust and drag force (further discussed in Section III.D), which is beneficial to ensure the warhead security. Convergence to the desired impact angle within a finite time is important in most practical guidance applications [29], which is conductive to the killing effects of the warheads in time-varying modern warfare.

C. FORMATION-CONTAINMENT TRACKING ANALYSIS AND DESIGN
In this section, velocity, positions and acceleration vectors take the form defined in the inertial coordinate system. Consider the following time-varying formation-containment tracking protocol where u i (t) = u i (t) (x) , u i (t) (y) , u i (t) (z) T is the acceleration along axial, vertical and lateral directions. K p = k p1 , k p2 (p = 1, 2), are constant gain matrices. τ and σ denote positive constants, and f (s(t)), s(t) ∈ R n is a nonlinear function to be ascertained later.
The following theorem provides sufficient conditions for system (21) with the proposed protocols to accomplish formation-containment tracking. Let T . Under protocols (29) and (30), multi-agent system (21) can be transformed as follows.
Theorem 2: The multiple tactical missile system under protocols (29) and (30) accomplish formation-containment tracking if the following conditions hold simultaneously.
(i) For all i ∈ E,ḣ ix (t) = h iv (t), that is, (ii) Let K 1 = ρ 1 B T P, and K 2 = ρ 2 B T P. Choose P is the positive definite solution to the following algebraic Riccati equation where s(t) ∈ R m represents (iv) Specify τ and σ in advance, and the values are relatively large with Define (t) and (t) as

The Formation-Leaders Can Achieve Formation Tracking:
Construct the Lyapunov candidate function as follows: where D E is a positive diagonal matrix described in Lemma2.
Taking the derivative of V E (t) along the trajectory of Eq. (40) obtainsV Based on Eq. (36), one has Thus, one can further has Under Assumption 1, one can obtain From Eq. (46) and (47), it holds thaṫ Substituting Eq. (34) and (37) into (48) leads tȯ Since Eq. (33) and (35) hold, it follows from (49) thaṫ Therefore, it can be concluded that (t) → 0 as t → ∞, that is, Eq. (22) is required. Furthermore, The Followers Can Achieve Containment: Construct the Lyapunov candidate function as follows: where D F is a positive diagonal matrix specified in Lemma2. The time derivative of V F (t) along the trajectories of Eq. (41) satisfieṡ One gets from Eq. (44) and (45) that Recalling that ff ( i (t)) ≤ 1, therefore Moreover, Noting that Eq. (34) holds, from Eq. (53), (54), and (55), it can be obtained thaṫ Under the condition that lim t→∞ (t) = 0, then it can be concluded that (t) → 0 as t → ∞. In other words, definition (23) is accomplished. The conclusion of Theorem 2 can be drawn.
Remark 4: The guidance law implemented by the trackingleader devised in Section III.B acts as an unknown input to formation-leaders. Note that the second term of Eq. (29) based on the sliding mode control provides an effective approach to make up for the tracking-leader's unknown input. To avoid undesirable chattering of the closed-loop system, the term s(t)/ s(t) in Eq. (36) can be replaced by s(t)/( s(t) +ω), whereω is a positive constant. Due to the unknown control input of the leader, the existing methods in [9,10] are no longer applicable for swarm systems in the absence of explicit knowledge of the formation-leaders' control input. Unlike the containment case in [11], there exist external inputs for each follower, thus the containment methods cannot be used directly in the formation-containment tracking case in this paper. The second term of Eq. (29) is utilized to suppress the external inputs from the formationleaders.
Remark 5: From the above analysis, reconnaissance missile outweighs attack missile on the maneuver capability, while decoy missile is more mobile than reconnaissance missile. In the overall design stage, the performance parameters should be designated consistent with the overload requirements of different types of missiles. More powerful maneuver capability can be guaranteed through relatively larger pneumatic actuators, lighter weight or engine with more specific impulse.

D. CONTROL SYSTEMS DESIGN FOR SALVO ATTACK OF TACTICAL MISSILES
Given the expected overload n xc , n yc , n zc in Section III.B and III.C, the rudder and axial overload controllers based on genetic algorithm and DOB combined with the transition process are proposed to drive the plant to the three direction reference overloads. The change of coordinate from the inertial system to the ballistic coordinate system is adopted through Eq. (4) in this section. Furthermore, to facilitate control systems design, the missile control model is decomposed into four channels, i.e., yawing channel, pitching channel, rolling channel and velocity channel.

Channel 1 Controller Design for Yawing Channel:
After the decomposition, we first propose the conventional two-loop overload autopilot with proportion-integral (PI) correction, which is arranged transition process in yawing channel. The following transition process during the transition period T 0 is constructed with respect to n zc devised in the preceding section.
The block diagram of the controller for yawing channel is shown in Fig.5. The controller in yawing channel is designed as Deflection angle of rudder δ y in Eq. (58) can control the actual β to track the required sideslip angle and finally steer n z ton zc , wheren zc is the output of a transition process in the form ofn zc = n zc trns (T 0 , t) and the command for the controller to keep up with. K p1 , K i1 and K d1 are controller parameters.
Note that the missile may be statically unstable for a sizable flight envelope, which makes the controller design domain narrow and difficult to find a workable solution. To cope with it, genetic algorithm is employed to find and optimize the optimal control parameters in this paper. The parameters selecting problem can be formulated in the following form: where J yaw denotes performance index, which can be calculated by where e yaw (t) =n zc − n z is the tracking error of the lateral overload. To ensure the transition process of the lateral overload fast and smoothly without overshoot, penalty terms for tracking error e yaw and deflection angle δ y are introduced. The weight coefficients for the performance index J yaw are selected as w yaw1 , w yaw2 and w yaw3 . Remark 6: The penalty terms for tracking error and deflection angle have the merits of the optimal flight control law such as good tracking, energy conserving. The last term of the second equation of Eq. (60) is utilized to alleviate the control effort when the tracking overshoot situation appears. Based on the above analysis, it is indicated that solving the genetic algorithm problem (59) with performance index (60) yields that each missile has good autopilot dynamic response capability, ensuring the control actuation system to track the guidance command without delay.
Channel 2 Controller Design for Pitching Channel: Similar to the yawing channel case, the controller in pitching channel is derived as Deflection angle of elevator δ z in Eq. (61) can control the actual α to track the required attack angle and finally steer n y ton yc , wheren yc is the output of a transition process in the form ofn yc = n yc trns (T 0 , t) and the command for the controller to track. K p2 , K i2 and K d2 are controller parameters to be determined appropriately. Fig.6 depicts the block diagram of the controller for pitching channel.  The following performance index on the base of genetic algorithm is also given as where e pitch (t) =n yc − n y is the tracking error of the vertical overload. w pitch1 , w pitch2 and w pitch3 are the weight coefficients for the performance index J pitch .

Channel 3 Controller Design For Rolling Channel:
For Skid-to-Turn (STT) missile, its rolling channel is generally held to be stable, hence, the command of the roll angle γ c is set to be zero. Design the control law as follows where actual and expected roll angel is γ , γ c = 0, respectively. δ x is the deflection angle of aileron. Selecting parameters K p3 , K d3 is obtained via the proposed approach. The block diagram of the controller for rolling channel is demonstrated by Fig.7. Design the following performance index for rolling channel where e roll (t) = γ c −γ , w roll1 , w roll2 and w roll3 are the weight coefficients. VOLUME 8, 2020

2) DESIGN OF AXIAL OVERLOAD CONTROL LOOP
Channel 4 Controller Design for Velocity Channel: Assumption 4: The fuel equivalence ratio affects the thrust T through an approximated linear relationship as shown in Eq. (65). T remains to be zero when is set to be zero. After is increased to 1, T changes to the maximum.
where T m is the maximum thrust of the powered missile.
In the light of assumption 4, the fuel equivalence ratio is controlled by the PI controller, which is designed as where K p and K i are control parameters to be chosen, T c is given reference thrust coincident with the above-mentioned expected axial overload (19), (28), (29) and (30). And it is worth mentioning that the axial overload for attack missile is set to be zero. Fig.8 presents the block diagram of the controller for velocity channel.
Considering the fact that the expected axial overload n xc is generated from the expressions of D in Eq. (28) and u i (t) in Eq. (29) and (30), it follows from the first equation of Eq. (12) that Obviously, if we obtain the exact information of the expected axial overload n xc , drag force X and the current flight attitude, then the reference thrust T c can be derived accordingly. However, due to modeling uncertainties and unpredictable disturbances in practical applications, drag force cannot be estimated accurately. For the uncertain drag force, we regard drag force as the lumped disturbance with d x = − X m , and there exists a positive constant d such that ḋ x ≤ d . Therefore, the velocity dynamics with external disturbance d x can be formulated aṡ Next, a DOB is proposed to estimate the external disturbance in the following theorem. Theorem 3: For the dynamics in velocity channel described by Eq. (68), the following DOB is constructed to estimate d x as where e V is the observation error of the velocity V , z 11 and z 12 are the observation value of V and d x , ε > 0 is a constant. The DOB for d x can guarantee that z 11 → V (t) and z 12 → d x (t) as t → ∞ if this following case holds. α 1 , α 2 are positive constants and the characteristic polynomial λ 2 + α 1 λ + α 2 meets Hurwitz conditions. Proof: Define the observing error as E = [e 1 , e 2 ] T , where e 1 = e V ε , e 2 = d x − z 12 . We can get form Eq. (69) that And the observing error equation can be written as where A e = −α 1 1 −α 2 0 and B e = 0 1 . The characteristic equation of matrix A e is λ 2 + α 1 λ + α 2 = 0. Based on the DOB parameter α 1 , α 2 and the Hurwitz conditions the characteristic polynomial λ 2 + α 1 λ + α 2 meets, there exist symmetric positive definite matrices P e and Q e , which satisfy the following algebraic Riccati equation Define the Lyapunov function as V 0 = εE T P e E > 0. Thus, there iṡ where λ min Q e is the minimum nonzero eigenvalue of matrix Q e . The positive constant ε is sufficiently small such that ε < λ min( Q e ) E 2 d P e B e , then it holds thatV 0 < 0, which means the estimation error for velocity channel based on DOB is asymptotically stable. This completes the proof of Theorem 3.
After obtaining the estimation of z 12 from DOB, replace d x with z 12 and then Eq. (67) is given by To comprehensively demonstrate the above-mentioned analysis, the overall control block diagram for velocity channel is shown in Fig.9. The final designed controller is composed by three equations, i.e., Eq. (65), (66), and (74).

Remark 7:
The turbojet engine thrust modeling is simplified as a linear model, and many details of the mathematical model are ignored. To deal with it, the effect of engine dynamics modeling error can be incorporated into the expression of d x . Although acceleration tracking controller design for multiple missiles is presented in [24], [25], the proposed cooperative controller is not applicable in that the explicit knowledge of the drag coefficient and the consideration of turbojet engine thrust were absent.

IV. NUMERICAL SIMULATIONS
In this section, the formation-containment tracking guidance and control of multiple tactical missiles to attack a stationary warship for berthing refueling in three dimensions is investigated to illustrate the effectiveness of the proposed strategies. The multiple tactical missile system is required to perform a salvo attack, where attack missile is designed to break through the interception of the anti-air defense system and strike the target with terminal impact angle constraint to achieve maximum lethality, and constitutes the formationcontainment tracking motion modes together with reconnaissance and decoy missiles.

A. ASCENARIO
Assume that there is a STT missile group with 1 attack missile, 3 reconnaissance missiles and 2 decoy missiles, which are set as the tracking-leader (labeled by 0), the formationleaders (labeled by 1-3) and the followers (labeled by 4,5) from the theoretical level, respectively. The simulations are conducted in the MATLAB2017b on a personal computer with the simulation step 0.002s and miss distance 10m. The missile model parameters are mass m = 735kg, moment of inertia J = diag (25, 350, 350) kg · m 2 , reference area S = 1m 2 , reference lengthc = 1m, and maximum thrust T m = 8kN . Refer to http://dx.doi.org/10.21227/ h3k0-b213 for the detailed data associated with the aerodynamic coefficients, and local aerodynamics can be obtained through linear interpolation. The initial conditions of 6 missiles and the warship are displayed in Table 2.
The directed communication graph of 6 missiles with 0-1 weights is demonstrated in Fig.10.
The attack mission of missile 0 is implemented in two different phases. In the first phase, during t ∈ [0s, 100s), missile 0 is at a trimmed cruise conditions: v = 270m/s, and h = 4000m. In the second phase, during t ∈ [100s, 117.5s],  missile 0 adjusts its boresight to aim at its target in the engagement scenario. An impact angle of q εd = −80 o and q βd = −60 o is selected to meet the requirement of the warheads to acquire better killing effects when hitting the stationary target.

B. BSIMULATION RESUTS
In the simulation, the guidance law parameters are designed to be β 1 = β 2 = 0.065, p 1 = p 2 = 5, q 1 = q 2 = 3, η 1 = η 2 = 0.008.Note that sign function sgn (s i ) may cause undesirable chattering, saturation function sat(s i ) is employed to eliminate the chattering. The effectiveness of the proposed approach in Section III.B can be illustrated in Fig.11. As is shown in Fig.11 (b), it can be concluded that missile 0 can attack the target precisely and the maximum miss distance is 0.1197m. Fig.11 (c) shows that the LOS angle q ε in vertical plane and q β in horizontal plane converge to about −80 o Based on Theorem 2, the control parameters in Eq. (34), (37) can be chosen as ρ 1 = 1.5, ρ 2 = 0.8, τ = 6, and η = 8. Fig.12 depicts the state trajectory snapshots of the 6 missiles at different time instants, where the tactical missiles from 0 to 5 are marked by black pentagram, green asterisk, blue asterisk, glaucous asterisk, red circle, and purple circle respectively, and the convex hull spanned by the states of missiles 1-3 is denoted by blue solid lines. The evolution of the formation-containment tracking geometric configuration can be clearly seen from Fig.12 (a-d).
According to the aerodynamic feature of the missile, the following aerodynamic coefficients around the trim condition are obtained using the linear interpolated lookup table.  To optimize the control parameters, genetic algorithm listed in Section III.D is applied with roughly given initial search values. Taking yawing channel for missile 0 as an example, let w yaw1 = 0.05, w yaw2 = 1, w yaw3 = 1, and the search scopes for parameters in yawing controller (57)  are set to be K p1 ∈ [0, 0.15], K i1 ∈ [0, 0.2] and K d1 ∈ [0, 2]. By setting the numbers of iterations and individuals, the optimal control parameters in yawing channel are converged asymptotically via the genetic algorithm such that K p1 = 0.1215, K i1 = 0.0531, K d1 = 1.6100. Similar to the yawing channel, it can also be calculated that the optimal control parameters in pitching controller (61) and rolling  controller (63) are K p2 = 0.1172, K i2 = 0.0209, K d2 = 1.5306, K p3 = 0.1109, K d3 = 0.0239. As exhibited in Fig.13 (a-c), after the control parameters adjustment, performance indexes in yawing, pitching and rolling channels reach and stay on the minimum conditions, respectively. Fig.13 (d-f) reveal that satisfactory system responses can be guaranteed within the solved controller parameter set. From Fig.14, it can be deduced that the actual overload in Y-direction and Z-direction can track respective command ideally due to the proposed controllers (58), (61), (63). Additionally, vertical and lateral overload command can be tracked through reasonable rudder deflection angles. Fig.15 demonstrates the estimation of the uncertain drag force of missile 0 by DOB (69). It can be concluded from Fig.15 (a) that the DOB can estimate the uncertain drag force with high precision, which is used to adjust the throttle setting of the turbojet engine through the controller (66). Thus, the overload of attack missile in X direction keeps zero according to the tangential balance between thrust and drag force, as displayed in Fig.15 (b). Fig.15 (c) depicts that the amplitude and transition of the throttle setting are quite reasonable and stable despite the uncertainty of the drag force estimation.
The trajectory of the position for the 6 missiles within 117.55s is shown in Fig.16. Fig.17 displays that the curves of the formation and containment errors approach to zero over time for each missile. The oscillations of formation tracking errors and containment errors at the moment t = 100s are induced by the switching from cruise phase to homing guidance. The effectiveness of the proposed DOB, the good tracking of the three direction reference overloads and the smooth transitions of flight parameters are exhibited in Fig.18. From Figs. 11-18, one can obtain the following inferences: (i) the states of the 3 formation-leaders preserve a parallel triangle formation while making tracks for the state trajectory of the tracking-leader, which realizes attack against a stationary target; (ii) the states of followers converge to the convex hull formed by those of the formation-leaders; (iii) the proposed DOB could estimate the uncertain drag force with high precision, and the desired reference overload can be effectively kept up with under this cooperative scheme through reasonable actuator actions. Therefore, the multiple tactical missile system can realize cooperative flight with the designated formation-containment tracking geometry.

V. CONCLUSIONS
XINGGUANG XU was born in 1988. He received the B.S. degree from Northwestern Polytechnical University, in 2011, started graduate study at Harbin Institute of Technology, and the M.S. degree from the third institute of China Aerospace Science and Industry Corporation, in 2014. He is currently pursuing the Ph.D. degree in navigation, guidance, and control with Beihang University (BUAA). He is also an Engineer with the Beijing Institute of Mechanical and Electrical Engineering. His research interests include aircraft design, fault-tolerant control, and cooperative control of multiagent systems.
CHANGRONG CHEN was born in 1994. He received the B.S. degree from the Harbin Institute of Technology, in 2017, and the M.S. degree from the third institute of China Aerospace Science and Industry Corporation, in 2020. His research interests include aircraft control systems design, adaptive control, and cooperative control of multiagent systems.
ZHANG REN received the B.S. degree in aircraft guidance, the M.S. degree in navigation, and the Ph.D. degree in simulation from Northwestern Polytechnical University, in 1982, 1985, and 1994, respectively. He held the Visiting Professor position at the University of California at Riverside and Louisiana State University, USA, respectively. He is currently a Professor with the School of Automation Science and Electronic Engineering, Beihang University. His research interests include aircraft guidance, navigation and control, fault detection, isolation and recovery, and cooperative control of multiagent systems.
SHUSHENG LI was born in 1962. He is currently a Professor and a Ph.D. Supervisor with the third institute of China Aerospace Science and Industry Corporation. He is also a National Outstanding Contribution Expert in the field of tracks and controls and granted a special allowance from the State Council. His research interests include aircraft design, fault diagnosis, tracks and controls, and design for testability.