Integrated Design of Multi-Constrained Snake Maneuver Surge Guidance Control for Hypersonic Vehicles in the Dive Segment

Abstract

Focusing on the large maneuver penetration of the hypersonic glide vehicle with multiple constraints and uncertain disturbance, a robust integrated guidance and control law, which can achieve the snake-shape maneuver, is designed. A snake-shape maneuver acceleration command, in the framework of sine function, determined by the altitude, target declination of the line of sight and the missile-target distance, is discussed. The integrated guidance and control law includes the terminal guidance law with multiple constraints, attitude control law and angular velocity control law. In the terminal guidance law design, the sliding mode control is adopted while the adaptive technique is applied to estimate the disturbance. The selected sliding mode surface has variable gain determined by the estimated time-to-go. With the designed terminal guidance law, using the snake-shape maneuver acceleration command as the bias item, the angular rate of the line of sight will converge to zero and the line of sight angle will converge to the expected value, simultaneously. The attitude control law and angular velocity control law are designed to track the expected attack and bank angles. The stability of the whole system is proved with the application of the Lyapunov theorem. The effectiveness and robustness of the proposed integrated guidance and control law is verified by simulation.

1. Introduction

The growing anti-missile technology has a greater threat for the penetration mission of the hypersonic glide vehicles (HGVs) [1,2,3]. The large maneuver penetration in dive trajectory is an effective strategy to decrease the interception probability of the anti-missile technology and improve the fight efficiency [4,5]. Theoretically, the large maneuver penetration of the HGV will force the interception weapons to cost more kinetic energy by making the line of sight (LOS) angular velocity and the required overload of the interception weapons change periodically. Therefore, it is beneficial for HGVs to achieve the large maneuver during the flight. Meanwhile, the guidance law with impact angle constraint has become a key solution to improving the damage capability of HGV [6,7,8].

At present, large maneuver penetration trajectories include spiral maneuver, jump maneuver and snake-shape maneuver [9]. There are three strategies to achieve the large maneuver penetration.

Firstly, the method based on the motion trajectory for a virtual target is discussed in Refs. [10,11,12,13]. With this strategy, at each moment, a virtual target point on the expected trajectory is defined and the guidance law is designed to track the expected trajectory of the virtual target point. A logarithmic-type spiral trajectory is adopted as a virtual target trajectory in Ref. [10] and an adaptive proportional guidance law is designed to achieve tracking to the virtual target trajectory. Focusing on the HGV’s spiral trajectory optimization problem with impact angle constraint, different spiral trajectories with different impact angles are obtained in Ref. [11] with the application of the continuous second-order cone programming (SOCP) method, in which the optimal target is to achieve the maximum falling velocity. Ref. [12] transforms the trajectory optimization problem of the HGV penetrating two interceptor weapons and reaching the target area to a SOCP problem and a penetration strategy, which only requires the initial LOS angular velocity information of the interceptors. It follows that a maneuver trajectory to penetrate two interceptor weapons is described. Some optimized trajectories are applied in deep neural network (DNN) training in Ref. [13] and the trained DNN is used as a maneuver command generator to achieve the real-time control. The problem of penetrating to two interceptor weapons is solved.

Another way to achieve penetration is to add an extra weighted maneuver overload on a traditional homing guidance law. As depicted in Ref. [14], a screw angular velocity is defined and a spiral maneuver overload command is generated. It is essentially a three-dimensional bias proportional guidance law. In Ref. [15], a sine-type overload command, designed on longitude and lateral planes, respectively, is treated as a bias item on a traditional overload command which is used to zero the LOS angular velocity. In Ref. [16], the principle of HGV penetration is firstly analyzed and the penetration model of the HGV and the interceptor is established. It transforms the HGV penetration problem into a trajectory optimization problem and the solver can be used to generate the overload command directly.

The periodic control based on the LOS angular velocity can also be applied to achieve the maneuver motion. A sliding mode guidance law, considering the auto-pilot dynamic performance simultaneously, is designed in Ref. [17] such that the LOS angular velocity can track certain sine-type signals and the snake-shaped maneuver is achieved. In Ref. [18], a kind of spiral trajectory, satisfying angle constraint, is designed, which can be transformed as the maneuver commands of the LOS angle and angular velocity. Furthermore, based on the feedback linearization technique, a sliding mode controller is designed to achieve the spiral maneuver guidance. Considering the terminal speed constraint, Ref. [19] establishes the relationship between the terminal speed and the maneuver amplitude based on predictive correction and the terminal multi-constraints maneuver guidance is achieved.

Although maneuver penetration is achieved in the mentioned references, there still exists some limitations. Refs. [10,12,13] do not consider the impact angle constraint after the penetration. The falling angle can be controlled in Ref. [11]; however, the azimuth angle constraint is not considered. The trajectory discussed in Ref. [14] seems to be fixed due to its non-adjustable maneuver amplitude. The spiral motion depicted in Ref. [15] is achieved based on the fixed lines in inertial reference, which limits the ability of the HGV large-range maneuver and target tracking. The feedback linearization in Ref. [18] decreases the robustness of the system due to the nonlinearities and uncertainties of the HGV. Furthermore, the penetration will make the coupling of the translation and rotation stronger and further increase the uncertainty of the system. However, the attitude control problems are ignored in Refs. [10,11,12,13,16,17].

Recently, integrated guidance and control (IGC) technique has attracted much attention due to its strong ability to deal with the coupling of the guidance and control systems of the missile. To solve the uncertainty problem and the full state constraints of the system, a dynamic surface control based IGC technique is proposed in Ref. [20]. A fuzzy adaptive fault tolerant IGC is discussed in Ref. [21] to deal with the actuator fault and multiple coupling uncertainties. An adaptive dynamic surface control with radial basis function neural network (RBFNN) estimating the uncertainties and external disturbances is proposed in Ref. [22]. A three-dimensional robust IGC, given in Ref. [23], has been proposed considering the impact angle constraint and input saturation, with the application of the extend state observer. However, it does not consider the penetration problem in dive phase. Ref. [24] proposes a high-order sliding mode observer to estimate the unknown uncertainties and states. The results of the observer can be used to compensate the disturbance. However, the above-mentioned techniques have great difficulties in engineering practice due to their complex structure.

In summary, as for the HGV’s large maneuver penetration in dive phase and the impact angle constraints, an innovative snake-shape maneuver IGC method is proposed. To begin with, a novel snake-shape maneuver model, whose amplitude is self-adjustable according to its height, is designed. With the height of the HGV becoming smaller, the maneuver amplitude becomes smaller to decrease the kinetic energy cost and guarantee the terminal speed requirement. Furthermore, the LOS declination constraint is considered in the snake-shape trajectory design which helps to achieve the snake-shape penetration and terminal impact angle constrained control. The proposed IGC, based on adaptive sliding mode control, can effectively deal with the strong coupling, uncertainties and disturbances of system caused by the large angle maneuver, which will increase the robustness of the system.

2. Preliminaries

In this section, the snake-shape maneuver model is firstly designed in the frame of sine-type function. Then, the nonlinear six degree-of-freedoms (DOF) dynamic model of the HGV is firstly given followed by the three-dimensional relative motion of the HGV and the target. Combining these two models, the integrated guidance and control model (IGCM), with strict-feedback form is then derived, fully considering the coupling of the aerodynamic and channels.

2.1. Snake-Shape Maneuver Model

In this manuscript, the snake-shape maneuver model is achieved in the frame of the sine-type function, which can be described as in Equation (1).

$$\begin{array}{l}{m}_{command}=A(h)\sin \left[\varpi ({\eta }_f,R)\right]\\ A(h)=\left\{\begin{array}{cc}{A}_0\left(h-h_s\right)/\left(h_0-h_s\right)& h\ge h_s\\ 0& h<h_s\end{array}\right.\\ \varpi ({\eta }_f,R)=\frac{2\pi \left(p+{\eta }_f/2{\pi }^2\right)}{R_0}\left|R_0-R\right|\end{array}$$

(1)

In Equation (1), $h$ is the current height of the HGV and $R$ is the relative distance of the HGV and the target. $h_0$ is the initial height when the snake-shape maneuver starts and $h_s$ is the terminal height when the snake-shape maneuver stops. $m_{command}$ is the command control force to achieve the snake-shape maneuver. The amplitude of the control force, $A(h)$, is determined by its current height $h$. The snake-shape maneuver is achieved only when $h\ge h_s$. The constant $A_0$ determines the amplitude of the control force at the initial height $h_0$. Note that usually the atmosphere density become larger while $A(h)$ will become smaller, with $h$ becoming smaller. Therefore, the designed snake-shape maneuver model will cost less energy. Let $R_0$ be the initial distance of $R$, the phase of the control force, $\varpi ({\eta }_f,p)$, is determined by the target LOS declination ${\eta }_f$. $p$ is a positive integer, indicating the cycle number of the trajectory.

2.2. Six DOF Dynamics of HGV

The adopted air-to-surface HGV in this assumed to be a rigid airframe in an plane-symmetric configuration [25]. The six-DOF system dynamics of the HGV is combined by translational and rotational dynamics, which can both be described by six differential equations. The translational dynamics of the can be given by [26]:

$$\left\{\begin{array}{l}\dot{x}=V\cos \theta \cos \phi \\ \dot{y}=V\sin \theta \\ \dot{z}=-V\cos \theta \sin \phi \\ m\dot{V}=-D-mg\sin \theta \\ mV\dot{\theta }=L\cos {\gamma }_v-N\sin {\gamma }_v-mg\cos \theta \\ -mV\cos \theta \dot{\phi }=L\sin {\gamma }_v+N\cos {\gamma }_v\end{array}\right.$$

(2)

where $x$, $y$ and $z$ are components of the position vector of the HGV in the inertial reference frame. $m$ is the mass and $V$ is the velocity value of HGV. $\theta$ and $\phi$ are the flight path angle and heading angle. ${\gamma }_v$ is bank angle. $g$ is the gravity acceleration. $D$, $L$ and $N$ are the aerodynamic drag, lift and lateral forces, which can be given by

$$\left\{\begin{array}{c}D=C_{D}qS\\ L=C_{L}qS\\ N=C_{N}qS\end{array}\right.$$

(3)

in which $C_D$, $C_L$ and $C_N$ are the aerodynamic drag, lift and lateral force coefficients, respectively. $q=0.5\rho V^2$ is the dynamic pressure and $S$ is the reference size. $\rho$ denotes the atmosphere density.

The angle of attack $\alpha$, sideslip angle $\beta$, bank angle ${\gamma }_v$ and the three-axis angular velocities, ${\omega }_x$, ${\omega }_y$, ${\omega }_z$ are adopted to describe the rotational dynamic model of the HGV, which can be given by:

$$\left\{\begin{array}{l}\dot{\alpha }={\omega }_z-\tan \beta ({\omega }_x\cos \alpha -{\omega }_y\sin \alpha )+\frac{mg\cos \theta \cos {\gamma }_v-L}{mV\cos \beta }\\ \dot{\beta }={\omega }_x\sin \alpha +{\omega }_y\cos \alpha +\frac{mg\cos \theta \sin {\gamma }_v+N}{mV}\\ {\dot{\gamma }}_v=\sec \beta ({\omega }_x\cos \alpha -{\omega }_y\sin \alpha )+\frac{L(\tan \beta +\tan \theta \sin {\gamma }_v)+N\tan \theta \cos {\gamma }_v-mg\cos \theta \tan \beta \cos {\gamma }_v}{mV}\\ {\dot{\omega }}_x=\frac{I_{yy}{M}_x+I_{xy}{M}_y}{I_{xx}{I}_{yy}-I_{xy}^2}+\frac{I_{xy}(I_{zz}-I_{xx}-I_{yy})}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_x{\omega }_z+\frac{I_{xy}^2+I_{yy}^2-I_{yy}{I}_{zz}}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_y{\omega }_z\\ {\dot{\omega }}_y=\frac{I_{xy}{M}_x+I_{xx}{M}_y}{I_{xx}{I}_{yy}-I_{xy}^2}+\frac{I_{xy}(I_{xx}+I_{yy}-I_{zz})}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_x{\omega }_z+\frac{I_{xx}{I}_{zz}-I_{xx}^2-I_{xy}^2}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_y{\omega }_z\\ {\dot{\omega }}_z=\frac{M_z}{I_{zz}}+\frac{I_{xy}}{I_{zz}}\left({\omega }_x^2-{\omega }_y^2\right)+\frac{I_{xx}-I_{yy}}{I_{zz}}{\omega }_x{\omega }_y\end{array}\right.$$

(4)

where $I_{xx}$, $I_{yy}$, $I_{zz}$ are the inertia moment of the HGV around $x$, $y$, $z$ axes of the body reference frame and $I_{xy}$ is the inertia product around the $xy$ plane. $M_x$, $M_y$, $M_z$ are the three-axis aerodynamic control torque elements which can be modeled as

$$\left\{\begin{array}{l}{M}_x=qSl{C}_{mx}\\ M_y=qSl{C}_{my}\\ M_z=qSl{C}_{mz}\end{array}\right.$$

(5)

in which $l$ is the reference length.

$C_{mx}$, $C_{my}$ and $C_{mz}$ are rolling, yaw and pitch torque coefficients, respectively. The aerodynamic force and torque coefficients, $C_L$, $C_N$, $C_{mx}$, $C_{my}$ and $C_{mz}$ are usually obtained by wind tunnel experiment or computational fluid dynamics (CFD) simulation. For convenience, these coefficients are curve fitted by a smooth formula, given in Appendix A, which is determined by Mach number $Ma$, angle of attack $\alpha$, sideslip angle $\beta$, and three-axis fin deflections ${\delta }_x$, ${\delta }_y$ and ${\delta }_z$. These aerodynamic coefficients can be generally denoted generally, according to Appendix A, as

$$\left\{\begin{array}{l}{C}_D=C_{D0}(Ma)+C_{D1}(Ma,\alpha ,\beta )+C_{D2}^x{\delta }_x\\ C_L=C_{L0}(Ma)+C_{L1}(Ma,\alpha ,\beta )+C_{L2}^z{\delta }_z\\ C_N=C_{N0}(Ma)+C_{N1}(Ma,\alpha ,\beta )+C_{N2}^y{\delta }_y\\ C_{mx}=C_{mx0}(Ma)+C_{mx1}(Ma,\alpha ,\beta )+C_{mx2}^x{\delta }_x+C_{mx2}^y{\delta }_y+C_{mx2}^z{\delta }_z\\ C_{my}=C_{my0}(Ma)+C_{my1}(Ma,\alpha ,\beta )+C_{my2}^x{\delta }_x+C_{my2}^y{\delta }_y+C_{my2}^z{\delta }_z\\ C_{mz}=C_{mz0}(Ma)+C_{mz1}(Ma,\alpha ,\beta )+C_{mz2}^x{\delta }_x+C_{mz2}^y{\delta }_y+C_{mz2}^z{\delta }_z\end{array}\right.$$

(6)

in which $C_{i0}(i=L,N,M,mx,my,mz)$ are the items only related to the constant item and Mach number Ma. $C_{i1}(i=L,N,M,mx,my,mz)$ are the items related to the angle of attack $\alpha$ and sideslip angle $\beta$, while $C_{i2}^j(i=L,N,M,mx,my,mz;j=x,y,z)$ denote the corresponding aerodynamic derivatives with respect to the fin deflections ${\delta }_x$, ${\delta }_y$ and ${\delta }_z$. According to the aerodynamic theorem, the angle of attack $\alpha$ and sideslip angle $\beta$ play major roles in producing the aerodynamic force while the fin deflections the aerodynamic torque.

2.3. Missile-Target Relative Motion Model

To begin with, we give the definition of the line of sight (LOS) reference frame to describe the relative motion of the HGV and the target. As described in Figure 1, the origin point of the LOS reference frame lies in the mass center of the HGV. x_s-axis points to the target. y_s-axis is perpendicular to x_s-axis, lying on the vertical plane. z_s-axis completes the right-hand reference frame.

To describe the strike process to the target of the HGV, the relative motion between the hypersonic vehicle and the target will be given in the LOS reference frame. The three-dimensional relative motion can be given as [26]

$$\left[\begin{array}{c}\ddot{R}\\ \ddot{\epsilon }\\ \ddot{\eta }\end{array}\right]=\left[\begin{array}{c}R{\dot{\epsilon }}^2+R{\dot{\eta }}^2\cos \epsilon \cos \epsilon \\ -2\dot{R}\dot{\epsilon }/R-{\dot{\eta }}^2\cos \epsilon \sin \epsilon \\ 2\dot{\epsilon }\dot{\eta }\tan \epsilon -2\dot{R}\dot{\eta }/R\end{array}\right]+\left[\begin{array}{ccc}1& 0& 0\\ 0& -1/R& 0\\ 0& 0& 1/R\cos \epsilon \end{array}\right]\left[\begin{array}{c}{a}_{hv}^{xs}\\ a_{hv}^{ys}\\ a_{hv}^{zs}\end{array}\right]$$

(7)

where R is the relative distance between the HGV and the target. $\epsilon$ is the LOS elevation angle and $\eta$ is the LOS azimuth angle. $a_{hv}^{xs}$, $a_{hv}^{ys}$, $a_{hv}^{zs}$ are the components of the acceleration of the HGV in the LOS reference frame. Let ${\mathit{C}}_v^g$ be the rotation matrix of the ground inertial reference frame relative to the ballistic reference frame while ${\mathit{C}}_g^s$ the rotation matrix of the LOS reference frame relative to the ground inertial reference frame. Then, the acceleration components, $a_{hv}^{xv}$, $a_{hv}^{yv}$, $a_{hv}^{zv}$ in ballistic reference frame, can be calculated according to the following relationship

$$\left[\begin{array}{l}{a}_{hv}^{xs}\\ a_{hv}^{ys}\\ a_{hv}^{zs}\end{array}\right]={\mathit{C}}_v^s\left[\begin{array}{l}{a}_{hv}^{xv}\\ a_{hv}^{yv}\\ a_{hv}^{zv}\end{array}\right]=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right]\left[\begin{array}{l}{a}_{hv}^{xv}\\ a_{hv}^{yv}\\ a_{hv}^{zv}\end{array}\right]$$

(8)

where ${\mathit{C}}_v^s={\mathit{C}}_g^s{\mathit{C}}_v^g$ is the rotation matrix of the ballistic reference frame relative to the LOS reference frame and $T_{ij}(i=1,2,3;j=1,2,3)$ is the corresponding element of ${\mathit{C}}_v^s$. According to Figure 1, one can obtain that

$${\mathit{C}}_g^s=\left[\begin{array}{ccc}\cos \epsilon \cos \eta & \sin \epsilon & -\cos \epsilon \sin \eta \\ -\sin \epsilon \cos \eta & \cos \epsilon & \sin \epsilon \sin \eta \\ \sin \eta & 0& \cos \eta \end{array}\right]$$

(9)

and ${\mathit{C}}_v^g$ can be given by [26]

$${\mathit{C}}_v^g=\left[\begin{array}{ccc}\cos \theta \cos \phi & -\sin \theta \cos \phi & \sin \phi \\ \sin \theta & \cos \theta & 0\\ -\cos \theta \sin \phi & \sin \theta \sin \phi & \cos \phi \end{array}\right]$$

(10)

Combining Equations (9) and (10), T_ij can be calculated as

$$\left\{\begin{array}{l}{T}_{11}=\cos \epsilon \cos \eta \cos \theta \cos \phi +\sin \epsilon \sin \theta +\cos \epsilon \sin \eta \cos \theta \sin \phi \\ T_{12}=-\cos \epsilon \cos \eta \sin \theta \cos \phi +\sin \epsilon \cos \theta -\cos \epsilon \sin \eta \sin \theta \sin \phi \\ T_{13}=\cos \epsilon \cos \eta \sin \phi -\cos \epsilon \sin \eta \cos \phi \\ T_{21}=-\sin \epsilon \cos \eta \cos \theta \cos \phi +\cos \epsilon \sin \theta -\sin \epsilon \sin \eta \cos \theta \sin \phi \\ T_{22}=\sin \epsilon \cos \eta \sin \theta \cos \phi +\cos \epsilon \cos \theta +\sin \epsilon \sin \eta \sin \theta \sin \phi \\ T_{23}=-\sin \epsilon \cos \eta \sin \phi +\sin \epsilon \sin \eta \cos \phi \\ T_{31}=\sin \eta \cos \theta \cos \phi -\cos \eta \cos \theta \sin \phi \\ T_{32}=-\sin \eta \sin \theta \cos \phi +\cos \eta \sin \theta \sin \phi \\ T_{33}=\sin \eta \sin \phi +\cos \eta \cos \phi \end{array}\right.$$

(11)

Since the sideslip angle of the HGV approximates zero when BTT maneuver is adopted, the lateral force N on the HGV can be ignored, which will be treated as uncertainty in control law design. An equivalent principle of the guidance law can be summarized as that the LOS elevation angle velocity $\stackrel{.}{\epsilon }$ and LOS azimuth angle $\stackrel{.}{\eta }$ converge to zero. Therefore, we mainly focus on the dynamics of the LOS elevation angle $\epsilon$ and LOS azimuth angle $\eta$. Thus, rearranging Equation (7), considering Equations (2), (8) and (11), will lead to

$$\left[\begin{array}{l}\ddot{\epsilon }\\ \ddot{\eta }\end{array}\right]=f_{\epsilon \eta }+\left[\begin{array}{cc}-1/R& 0\\ 0& 1/R\cos \epsilon \end{array}\right]\left[\begin{array}{cc}{T}_{22}& T_{23}\\ T_{32}& T_{33}\end{array}\right]\left[\begin{array}{l}L\cos {\gamma }_v/m\\ L\sin {\gamma }_v/m\end{array}\right]+{\Delta }_s$$

(12)

where the nonlinear function $f_{\epsilon \eta }$ is

$$f_{\epsilon \eta }=\left[\begin{array}{l}-2\dot{R}\dot{\epsilon }/R-{\dot{\eta }}^2\cos \epsilon \sin \epsilon \\ -2\dot{R}\dot{\eta }/R+2\dot{\epsilon }\dot{\eta }\tan \epsilon \end{array}\right]$$

(13)

and the uncertain disturbance ${\Delta }_s$ is

$${\Delta }_s=\left[\begin{array}{cc}-1/R& 0\\ 0& 1/R\cos \epsilon \end{array}\right]\left[\begin{array}{ccc}{T}_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right]\left[\begin{array}{c}-D/m-g\sin \theta \\ -N\sin {\gamma }_v/m-g\cos \theta \\ N\cos {\gamma }_v/m\end{array}\right]$$

(14)

2.4. Integrated Model of Guidance and Control

To facilitate deriving the IGCM and control law, we hoped to rewrite the dynamics of the HGV in vector form. As for the attitude dynamics, recalling the first three equations of Equation (4), it can be rewritten as

$$\left[\begin{array}{c}\dot{\alpha }\\ \dot{\beta }\\ \dot{\gamma }\end{array}\right]=f_{\alpha \beta \gamma }+\left[\begin{array}{ccc}-\cos \alpha \tan \beta & \sin \alpha \tan \beta & 1\\ \sin \alpha & \cos \alpha & 0\\ \cos \alpha \sec \beta & -\sin \alpha \sec \beta & 0\end{array}\right]\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]$$

(15)

in which

$$f_{\alpha \beta \gamma }=\left[\begin{array}{c}\frac{mg\cos \theta \cos {\gamma }_v-L}{mV\cos \beta }\\ \frac{N+m\mathrm{gcos}\theta \sin {\gamma }_v}{mV}\\ \frac{L(\tan \beta +\tan \theta \sin {\gamma }_v)+N\tan \theta \cos {\gamma }_v-m\mathrm{gcos}\theta \tan \beta \cos {\gamma }_v}{mV}\end{array}\right]$$

(16)

Similarly, the last three equations of Equation (4) can also be rearranged as

$$\left[\begin{array}{c}{\dot{\omega }}_x\\ {\dot{\omega }}_y\\ {\dot{\omega }}_z\end{array}\right]=f_{\omega xyz}+qSl\left[\begin{array}{ccc}\frac{I_{yy}}{I_{xx}{I}_{yy}-I_{xy}^2}& \frac{I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}& 0\\ \frac{I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}& \frac{I_{xx}}{I_{xx}{I}_{yy}-I_{xy}^2}& 0\\ 0& 0& \frac{1}{I_{zz}}\end{array}\right]\left[\begin{array}{ccc}{C}_{mx2}^x& C_{mx2}^y& C_{mx2}^z\\ C_{my2}^x& C_{my2}^y& C_{my2}^z\\ C_{mz2}^x& C_{mz2}^y& C_{mz2}^z\end{array}\right]\left[\begin{array}{l}{\delta }_x\\ {\delta }_y\\ {\delta }_z\end{array}\right]+\Delta r$$

(17)

where $\Delta r$ represents the uncertainty of the IGCM caused by aerodynamic coupling and external disturbance which can be given as

$$\Delta r=qSl\left[\begin{array}{ccc}\frac{I_{yy}}{I_{xx}{I}_{yy}-I_{xy}^2}& \frac{I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}& 0\\ \frac{I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}& \frac{I_{xx}}{I_{xx}{I}_{yy}-I_{xy}^2}& 0\\ 0& 0& \frac{1}{I_{zz}}\end{array}\right]\left(\left[\begin{array}{l}{C}_{mx0}(Ma)+C_{mx1}(Ma,\alpha ,\beta )\\ C_{my0}(Ma)+C_{my1}(Ma,\alpha ,\beta )\\ C_{mz0}(Ma)+C_{mz1}(Ma,\alpha ,\beta )\end{array}\right]\right)$$

(18)

The nonlinear function $f_{\omega xyz}$ represents the coupling components of the three-axis channel coupling with the formulation as

$$f_{\omega xyz}=\left[\begin{array}{l}{\left(I_{xx}{I}_{yy}-I_{xy}^2\right)}^{-1}\left(I_{xy}(I_{zz}-I_{xx}-I_{yy}){\omega }_x{\omega }_z+(I_{xy}^2+I_{yy}^2-I_{yy}{I}_{zz}){\omega }_y{\omega }_z\right)\\ {\left(I_{xx}{I}_{yy}-I_{xy}^2\right)}^{-1}\left(I_{xy}(I_{xx}+I_{yy}-I_{zz}){\omega }_x{\omega }_z+(I_{xy}^2+I_{yy}^2-I_{yy}{I}_{zz}){\omega }_y{\omega }_z\right)\\ I_{zz}^{-1}\left(I_{xy}({\omega }_x^2-{\omega }_y^2)+(I_{xx}-I_{yy}){\omega }_x{\omega }_y\right)\end{array}\right]$$

(19)

It is shown that Equations (15) and (17) describe the three-dimensional attitude motion of the HGV with the strict-feedback form considering the coupling of the aerodynamic and channels.

According to parallel approaching method, zeroing the LOS angle velocity will lead to zero miss distance, i.e., the LOS direction will keep constant ultimately [27]. Additionally, considering the impact angle constraint to enhance the attack effectiveness, it is usually required that the LOS angle velocity $\dot{\epsilon }$ and $\dot{\eta }$ will eventually converge to zero while the LOS angle $\epsilon$ and $\eta$ the desired LOS angle ${\epsilon }_f$ and ${\eta }_f$ [8]. Therefore, to describe the integrated dynamics of the HGV succinctly, define ${\mathit{x}}_1={\left[\begin{array}{cc}\epsilon -{\epsilon }_f& \eta -{\eta }_f\end{array}\right]}^T$, ${\mathit{x}}_2={\left[\begin{array}{cc}\dot{\epsilon }& \dot{\eta }\end{array}\right]}^T$, ${\mathit{x}}_3={\left[\begin{array}{ccc}\alpha & \beta & {\gamma }_v\end{array}\right]}^T$ and ${\mathit{x}}_4={\left[\begin{array}{ccc}{\omega }_x& {\omega }_y& {\omega }_z\end{array}\right]}^T$ in which ${\epsilon }_f$ and ${\eta }_f$ are the target LOS elevation and azimuth angles, respectively. $\mathit{u}={\left[\begin{array}{ccc}{\delta }_x& {\delta }_y& {\delta }_z\end{array}\right]}^T$ is the control input. Rewriting Equations (12), (15) and (17), the following three-dimensional IGCM can be obtained as

$$\left\{\begin{array}{l}{\dot{\mathit{x}}}_1={\mathit{x}}_2\\ {\dot{\mathit{x}}}_2={\mathit{f}}_{\epsilon ,\eta }+{\mathit{g}}_1{\mathit{x}}_3^{\#}+{\Delta }_s\\ {\dot{\mathit{x}}}_3={\mathit{f}}_{\alpha \beta \gamma }+{\mathit{g}}_2{\mathit{x}}_4\\ {\dot{\mathit{x}}}_4={\mathit{f}}_{\omega xyz}+{\mathit{g}}_3\mathit{u}+{\Delta }_r\end{array}\right.$$

(20)

where ${\mathit{g}}_1$, ${\mathit{g}}_2$ and ${\mathit{g}}_3$ are control gain matrices defined by

$${\mathit{g}}_1=\left[\begin{array}{cc}-1/mR& 0\\ 0& 1/mR\cos \epsilon \end{array}\right]\left[\begin{array}{cc}{T}_{22}& T_{23}\\ T_{32}& T_{33}\end{array}\right]$$

(21)

$${\mathit{g}}_2=\left[\begin{array}{ccc}-\cos \alpha \tan \beta & \sin \alpha \tan \beta & 1\\ \sin \alpha & \cos \alpha & 0\\ \cos \alpha \sec \beta & -\sin \alpha \sec \beta & 0\end{array}\right]$$

(22)

$${\mathit{g}}_3=q{S}_r{L}_r\left[\begin{array}{ccc}{\left(I_{xx}{I}_{yy}-I_{xy}^2\right)}^{-1}{I}_{yy}& {\left(I_{xx}{I}_{yy}-I_{xy}^2\right)}^{-1}{I}_{xy}& 0\\ {\left(I_{xx}{I}_{yy}-I_{xy}^2\right)}^{-1}{I}_{xy}& {\left(I_{xx}{I}_{yy}-I_{xy}^2\right)}^{-1}{I}_{xx}& 0\\ 0& 0& I_{zz}^{-1}\end{array}\right]\left[\begin{array}{ccc}{m}_{x,{\delta }_x}& m_{x,{\delta }_y}& m_{x,{\delta }_z}\\ m_{y,{\delta }_x}& m_{y,{\delta }_y}& m_{y,{\delta }_z}\\ m_{z,{\delta }_x}& m_{z,{\delta }_y}& m_{z,{\delta }_z}\end{array}\right]$$

(23)

and

$${\mathit{x}}_3^{\#}=\left[\begin{array}{l}L\cos {\gamma }_v/m\\ L\sin {\gamma }_v/m\end{array}\right]$$

(24)

According to Equation (20), the established IGC model has a strict-feedback form with unmatched uncertainty Δ_s and matched uncertainty ${\Delta }_r$. Apparently, if LOS angular rate vector ${\mathit{x}}_2$ converges to zero, HGV will fly towards the target. Additionally, HGV will strike the target with the desired target LOS inclination and declination angles when ${\mathit{x}}_1$ converges to zero. Before formulating the control law, the following assumptions are conducted.

Assumption 1. 

The uncertainties ${\Delta }_s$ and ${\Delta }_r$ in the IGC model of HGV are unknown but bounded. The upper bounds of these uncertainties are also unknown. That is to say, there exist unknown positive constants  $M_s$ and  $M_r$ such that $‖{\Delta }_s‖\le M_s$ and $‖{\Delta }_r‖\le M_r$, where $‖·‖$ represents the Euclidean norm.

Assumption 2. 

The control gain matrices ${\mathit{g}}_1$, ${\mathit{g}}_2$ and ${\mathit{g}}_3$ are smooth and bounded.

Assumption 1 is evident when the states of Equation (20) is bounded, which always holds for bounded control input. As for Assumption 2, the smoothness and boundedness of these matrices can be guaranteed by smoothness and boundedness of the states of Equation (20).

3. The Maneuver Penetration IGC Design Considering the Terminal Impact Angle Constraint

The goals of the maneuver penetration IGC design considering the terminal impact angle constraint can be summarized as follows: firstly, the HGV will hit the target accurately with a target terminal impact angle; secondly, all states of Equation (20) are bounded, stable and controllable; finally, the HGV will maintain high-robustness during the whole flight procedure. Mathematically, according to Equation (20) the above mentioned control goals can also be described to design the control law for angles of fin deflections $\mathit{u}$ such that ${\lim }_{t\to t_f}{\mathit{x}}_1=0$ and ${\lim }_{t\to t_f}{\mathit{x}}_2=0$, and all signals in Equation (20) are bounded, when considering the uncertainties.

As depicted in Figure 2, the control law is accomplished through three-loop control with the following steps. Firstly, design the terminal guidance law considering impact angle constraints. According to the relative motion of the HGV and the target, design the target angle of attack ${\alpha }_c$ and bank angle ${\gamma }_{vc}$, which can provide the desired acceleration given in the guidance law; then, design the attitude control law. Treating the angular velocity vector $\mathsf{\omega }$ as the virtual control input ${\mathsf{\omega }}_c$, track the target angle of attack ${\alpha }_c$ and bank angle ${\gamma }_{vc}$, to provide the command aerodynamic force. Finally, design the angular velocity control law. Design the control law for angles of fin deflections $\mathit{u}$ such that the actual angular velocity will coincide with the expected value ${\mathsf{\omega }}_c$.

3.1. The Guidance Law Design with Snake-Shaped Trajectory Considering Terminal Impact Angle Constraint

Consider the first two equations of Equation (20) as

$$\left\{\begin{array}{l}{\dot{\mathit{x}}}_1={\mathit{x}}_2\\ {\dot{\mathit{x}}}_2={\mathit{f}}_{\epsilon ,\eta }+{\mathit{g}}_1{\mathit{x}}_3^{\#}+{\Delta }_s\end{array}\right.$$

(25)

where ${\mathit{x}}_3^{\#}$ can be treated as the virtual control input. The sliding mode control technique is adopted such that ${\lim }_{t\to t_f}{\mathit{x}}_1=0$ and ${\lim }_{t\to t_f}{\mathit{x}}_2=0$ considering uncertainties ${\Delta }_s$. The sliding mode surface is selected as

$${\mathit{S}}_1=\left[\begin{array}{l}\dot{\epsilon }\\ \dot{\eta }\end{array}\right]+\frac{k_{ls}}{T_{go}}\left[\begin{array}{l}\epsilon -{\epsilon }_f\\ \eta -{\eta }_f\end{array}\right]={\mathit{x}}_2+\frac{k_{ls}}{T_{go}}{\mathit{x}}_1$$

(26)

where $k_{ls}>0$ is the selected sliding mode gain and $T_{go}=-r/\dot{r}$ is the estimated time-to-go. It can be analyzed that ${\mathit{S}}_1=0$ will lead to ${\mathit{x}}_1=\exp (-k_{ls}(t-t_0)/T_{go}){\mathit{x}}_1(t_0)$. That is to say, when the states of Equation (22) come to the sliding mode surface, it will eventually converge to the equilibrium, i.e., $\dot{\epsilon }\to 0$, $\dot{\eta }\to 0$, $\epsilon \to {\epsilon }_f$ and $\eta \to {\eta }_f$. Note that in Equation (22), the variable coefficient factor which varies inversely with the time-to-go is adopted to control the convergence rate. It is apparent that a greater $T_{go}$ will lead to a slower convergence rate. Therefore, in the first stage of the flight, ${\mathit{S}}_1$ plays a major role in driving the LOS angular rate ${\mathit{x}}_2$ to be zero while in the terminal stage of the flight, with the $T_{go}$ becoming smaller, ${\mathit{S}}_1$ acts in achieving the terminal impact angle constraint.

If the aerodynamic force is treated as a virtual input firstly, adopting the exponential reaching law, the adaptive control law can be designed as

$$\left[\begin{array}{l}{L}_e\cos {\gamma }_{ve}\\ L_e\sin {\gamma }_{ve}\end{array}\right]=m{\mathit{g}}_1^{-1}\left(-\frac{{\epsilon }_1}{T_{go}}{\mathit{S}}_1-\frac{K_1}{T_{go}}\frac{{\mathit{S}}_1}{‖{\mathit{S}}_1‖}-\frac{k_{ls}}{T_{go}}{\mathit{x}}_2-{\mathit{f}}_{\epsilon ,\eta }-{\widehat{M}}_s\frac{{\mathit{S}}_1}{‖{\mathit{S}}_1‖}\right)$$

(27)

where intermediate variables $L_e$ and ${\gamma }_{ve}$ represent the expected lift force and bank angle, respectively, without considering the snake-shape maneuver. Positive constants ${\epsilon }_1$ and $K_1$ are to be designed control parameters. ${\widehat{M}}_s$ is the estimated value of the upper bound $M_s$ of the disturbance ${\Delta }_s$, which is updated by the following adaptive law:

$${\dot{\widehat{M}}}_s={\zeta }_1‖{\mathit{S}}_1‖$$

(28)

in which ${\zeta }_1$ is a positive control parameter to be selected.

To analyze the stability of the system, a candidate Lyapunov function, $V_{11}$, for the subsystem Equation (25) is selected as

$$V_{11}=\frac{1}{2}{S}_1^T{S}_1+\frac{1}{2{\zeta }_1}{({\widehat{M}}_s-M_s)}^2$$

(29)

Direct computation yields that

$$\begin{array}{ll}{\dot{V}}_{11}& ={\mathit{S}}_1^T{\dot{\mathit{S}}}_1+\frac{1}{{\zeta }_1}({\widehat{M}}_s-M_s){\dot{\widehat{M}}}_s\\ & ={\mathit{S}}_1^T\left(-\frac{{\epsilon }_1}{T_{go}}{\mathit{S}}_1-\frac{K_1}{T_{go}}\frac{{\mathit{S}}_1}{‖{\mathit{S}}_1‖}-{\widehat{M}}_s\frac{{\mathit{S}}_1}{‖{\mathit{S}}_1‖}+{\Delta }_s\right)+({\widehat{M}}_s-M_s)‖{\mathit{S}}_1‖\\ & \le {\mathit{S}}_1^T\left(-\frac{{\epsilon }_1}{T_{go}}{\mathit{S}}_1-\frac{K_1}{T_{go}}\frac{{\mathit{S}}_1}{‖{\mathit{S}}_1‖}\right)+M_s‖{\mathit{S}}_1‖-{\widehat{M}}_s‖{\mathit{S}}_1‖+({\widehat{M}}_s-M_s)‖{\mathit{S}}_1‖\\ & \le -\frac{{\epsilon }_1}{T_{go}}{‖{\mathit{S}}_1‖}^2-\frac{K_1}{T_{go}}‖{\mathit{S}}_1‖\le 0\end{array}$$

(30)

It will eventually come to the conclusion that the subsystem Equation (25) is asymptotically stable under the virtual control law Equations (27) and (28), with the application of Barbalat’s Lemma [28].

The control law, given as Equations (27) and (28), is designed based on zeroing LOS angular rate ${\mathit{x}}_2$ with a certain terminal impact angle ${\epsilon }_f$ and ${\eta }_f$. Therefore, the trajectory of the HGV is straight. To achieve the snake-shaped trajectory, an additional item is considered in Equation (27). The new control law is constructed as

$$\left[\begin{array}{l}{L}_c\cos {\gamma }_{vc}\\ L_c\sin {\gamma }_{vc}\end{array}\right]=\left[\begin{array}{l}{L}_e\cos {\gamma }_{ve}\\ L_e\sin {\gamma }_{ve}\end{array}\right]+m{\mathit{m}}_c$$

(31)

where $L_c$ and ${\gamma }_{vc}$ are the target lift force and bank angle, considering the snake-shape maneuver problem. The bias item ${\mathit{m}}_c$ is

$${\mathit{m}}_c=\left\{\begin{array}{cc}{[0 A_0\frac{h-h_s}{h_0-h_s}\sin \left(\frac{2\pi (p+{\eta }_f/2{\pi }^2)}{R_0}\left|R_0-R\right|\right)]}^T& h\ge h_s\\ {\left[\begin{array}{cc}0& 0\end{array}\right]}^T& h<h_s\end{array}\right.$$

(32)

Due to the nonlinear relationship between the angle of attack $\alpha$ and the bank angle ${\gamma }_v$, the numerical method will be applied to obtain the target angle of attack $\alpha$ and the bank angle ${\gamma }_v$. According to the wind tunnel experiment, we can obtain the table of the lift coefficient $c_L(\alpha ,Ma)$ with respect to the angle of attack $\alpha$ and the Mach number $Ma$. To have a faster convergence rate, Newton-Raphson’s method can be applied in finding the target angle of attack and bank angle (Algorithm 1).

Algorithm 1: Determine the target angle of attack ${\alpha }_c$ and the target bank angle${\gamma }_{vc}$
Input: the expected aerodynamic forces $F_{cx}=L_c\cos {\gamma }_{vc}$, $F_{cy}=L_c\sin {\gamma }_{vc}$ according to Equation (24); current Mach number $Ma$; lift coefficient table $c_L(\alpha ,Ma)$; iterative tolerance ε > 0; dynamic pressure $q$, reference size $S$. Output: the target angle of attack ${\alpha }_c$ and the bank angle ${\gamma }_{vc}$
	1. Initial: ${\alpha }_0=0$;
	2. Calculate the expected aerodynamic lift force coefficient $c_{Lc}=\sqrt{F_{cx}^2+F_{cy}^2}/qS$;
	3. According to Mach number $Ma$ and current angle of attack ${\alpha }_k$, find the aerodynamic lift force coefficient $c_L({\alpha }_k,Ma)$ and its differential coefficient ${c^{\prime }}_L({\alpha }_k,Ma)$;
	4. Update angle of attack: ${\alpha }_{k+1}={\alpha }_k-\frac{c_L({\alpha }_k,Ma)-c_{Lc}}{{c^{\prime }}_L({\alpha }_k,Ma)}$;
	5. Determine if $\left|{\alpha }_{k+1}-{\alpha }_k\right|\le ϵ$, if not, return to Step 1;
	6. Get the target angle of attack ${\alpha }_c={\alpha }_{k+1}$;
	7. Calculate the target bank angle ${\gamma }_{vc}=\arctan (F_{cy},F_{cx})$;
	8. End.

3.2. The Attitude Control Law Design

In Section 3.1, the guidance law to achieve the snake-shaped trajectory is designed and the target angle of attack ${\alpha }_c$ and the bank angle ${\gamma }_{vc}$ are obtained. For the third equation of Equation (20), treating the angular velocity ${\mathit{x}}_4$ as the virtual control input, it is expected to design an attitude control law ${\mathit{x}}_{4c}={\left[\begin{array}{ccc}{\omega }_{xc}& {\omega }_{yc}& {\omega }_{zc}\end{array}\right]}^T$ such that ${\mathit{x}}_3$ will coincide with ${\mathit{x}}_{3c}={\left[\begin{array}{ccc}{\alpha }_c& 0& {\gamma }_{vc}\end{array}\right]}^T$.

$${\dot{\mathit{x}}}_3={\mathit{f}}_{\alpha \beta \gamma }+{\mathit{g}}_2{\mathit{x}}_4$$

(33)

Define the tracking error as

$${\mathit{S}}_2={\mathit{x}}_3-{\mathit{x}}_{3c}={\left[\begin{array}{ccc}\alpha -{\alpha }_c& \beta & {\gamma }_v-{\gamma }_{vc}\end{array}\right]}^T$$

(34)

If the virtual control law is designed as

$${\mathit{x}}_{4c}={\mathit{g}}_2^{-1}\left(-{\epsilon }_2{\mathit{S}}_2-{\mathit{f}}_{\alpha \beta \gamma }+{\dot{\mathit{x}}}_{3c}\right)$$

(35)

where ${\epsilon }_2$ is the control parameter to be designed, then the tracking error ${\mathit{S}}_2$ will converge to zero asymptotically.

3.3. The Angular Velocity Control Law Design

In this section, the control law for fin deflections will be designed to track the virtual angular velocity ${\mathit{x}}_{4c}$ which has been discussed in Section 3.2. For the last equation of Equation (20), define the tracking error as

$${\mathit{S}}_3={\mathit{x}}_4-{\mathit{x}}_{4c}={\left[\begin{array}{ccc}{\omega }_x-{\omega }_{xc}& {\omega }_y-{\omega }_{yc}& {\omega }_z-{\omega }_{zc}\end{array}\right]}^T$$

(36)

whose dynamics can be given by

$${\dot{\mathit{S}}}_3={\dot{\mathit{x}}}_4-{\dot{\mathit{x}}}_{4c}={\mathit{f}}_{\omega xyz}+{\mathit{g}}_3\mathit{u}+{\Delta }_r-{\dot{\mathit{x}}}_{4c}$$

(37)

To achieve the tracking and eliminate the influence of the uncertainties, the adaptive control law is designed as

$$\mathit{u}={\left[\begin{array}{ccc}{\delta }_x& {\delta }_y& {\delta }_z\end{array}\right]}^T={\mathit{g}}_3^{-1}\left(-{\mathit{K}}_3\frac{{\mathit{S}}_3}{‖{\mathit{S}}_3‖}-{\epsilon }_3{\mathit{S}}_3-{\mathit{f}}_{\omega xyz}-{\widehat{M}}_r\frac{{\mathit{S}}_3}{‖{\mathit{S}}_3‖}+{\dot{\mathit{x}}}_{4c}\right)$$

(38)

where ${\mathit{K}}_3,{\mathsf{\epsilon }}_3\in {ℝ}^{3\times 3}$ are two positive definite matrices to be selected and ${\widehat{M}}_r$ is the estimate value of the upper bound $M_r$ of the uncertainty ${\Delta }_r$. The update law of ${\widehat{M}}_r$ is

$${\dot{\widehat{M}}}_r={\zeta }_2{\mathit{S}}_3$$

(39)

in which ${\zeta }_2$ is a positive control parameter to be selected.

The stability of the subsystem Equation (37) can be proved by selecting a Lyapunov function as

$$V_{13}=\frac{1}{2}{\mathit{S}}_3^T{\mathit{S}}_3+\frac{1}{2{\zeta }_2}{({\widehat{M}}_r-M_r)}^2$$

(40)

The time differentiation of $V_{13}$ along the trajectories of Equation (37) can be given as

$$\begin{array}{ll}{\dot{V}}_{13}& ={\mathit{S}}_3^T{\mathit{S}}_3+\frac{1}{{\zeta }_2}({\widehat{M}}_r-M_r){\dot{\widehat{M}}}_r\\ & ={\mathit{S}}_3^T\left(-{\mathit{K}}_3\frac{{\mathit{S}}_3}{‖{\mathit{S}}_3‖}-{\mathsf{\epsilon }}_3{\mathit{S}}_3-{\widehat{M}}_r\frac{{\mathit{S}}_3}{‖{\mathit{S}}_3‖}+{\Delta }_r\right)+({\widehat{M}}_r-M_r){\mathit{S}}_3\\ & =-{\lambda }_{\min }({\mathit{K}}_3)‖{\mathit{S}}_3‖-{\lambda }_{\min }({\mathsf{\epsilon }}_3){‖{\mathit{S}}_3‖}^2\le 0\end{array}$$

(41)

Similar to the analysis in Section 3.1, with the application of Barbalat’s Lemma [28], one can conclude that the designed control law Equations (38) and (39) can guarantee the asymptotical stability of the subsystem (37).

Remark 1: 

The stability analysis is proved for the three subsystems. According to Section 3.1, Section 3.2 and Section 3.3, each subsystem is all asymptotically stable with the application of the control law in each stable. Therefore, the stability of the outer loop subsystem can be guaranteed only when the inner loop provides the virtual command exactly, i.e., no boundary layer error happens, which cannot be satisfied due to the dynamic performance of the inner loop subsystem. The stability of the whole system is analyzed in Appendix C.

Remark 2: 

The existence of the switch items in Equations (27) and (38) makes the control law unsmooth, which may lead to the chattering problem. To avoid this, the functions  ${\mathit{S}}_i/‖{\mathit{S}}_i‖$ will all be replaced by  ${\mathit{S}}_i/\left(‖{\mathit{S}}_i‖+{\delta }_i\right)$, where  ${\delta }_i$ is a small constant. The property that  ${\mathit{S}}_i/\left(‖{\mathit{S}}_i‖+{\delta }_i\right)$ is continuous and smooth implies that the control law will not encounter chattering.

4. Numerical Simulation

4.1. Feasibility Simulation

To demonstrate the effectiveness of the proposed snake-shape maneuver integrated guidance and control (SMIGC) technique, a numerical simulation is conducted for an HGV whose profile parameters and aerodynamic parameters can be found in Appendix A and

Table A1. The general parameters of HGV.

Variable Symbol	Value	Meaning
$m$	2200 kg	Mass of HGV
$S$	6.33 m2	Reference area
$L$	3.6 m	Reference length
$I_{xx}$	345 kg·m2	Moment if inertia around $x$ axis
$I_{yy}$	4917.1 kg·m2	Moment if inertia around $y$ axis
$I_{zz}$	4956.6 kg·m2	Moment if inertia around $z$ axis
$I_{xy}$	1302.4 kg·m2	Product of inertia

. Assume that the HGV locates at $(L_0,B_0,H_0)=\left(96°,46°,30 \mathrm{km}\right)$ initially where $L_0$, $B_0$ and $H_0$ are the initial values of the longitude, latitude and height of the HGV. The longitude, latitude and height of the target are $(L_t,B_t,H_t)=\left(95°,45°,0\right)$. The initial velocity of HGV is 2200 m/s. The target LOS angles are $({\epsilon }_f,{\eta }_f)=(-70°,-20°)$. The initial attitude and angular velocity are all assumed to be zero. The control parameters are shown in

Table 1. Control parameters.

Parameter	Value	Parameter	Value
Kls	diag [5, 15]	K1	diag [0.01, 0.01]
K2	diag [0.001, 0.01, 0.001]	K3	diag [0.01, 0.001, 0.01]
ε1	diag [10, 10]	ε2	diag [15, 32, 10]
ε3	diag [20, 20, 55]	A0	1
P	2

. In addition, a comparison with a recent angle constrained integrated guidance and control (ACIGC) [23] is conducted.

To avoid a saturation problem in the simulation, some parameters are limited in certain ranges as illustrated in

Table 2. Parameter limitations.

Parameter	Value	Parameter	Value
Fin angle limitation	[−20°, 20°]	Fin angle rate limitation	[−100°/s, 100°/s]
Angle of attack limitation	[−4°, 20°]	Sideslip angle limitation	[−5°, 5°]
Angle of attack rate limitation	[−20°/s, 20°/s]	Sideslip angle rate limitation	[−20°/s, 20°/s]
Bank angle rate limitation	[−60°/s, 60°/s]

. The simulation will be terminated when the height of the HGV is lower than zero, i.e., $H<0$. The simulation results are shown in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

From the three-dimensional trajectory of the HGV, as depicted in Figure 3, it is shown that the HGV achieves snake-shaped penetration and hits the target successfully under the proposed SMIGC technique. It can also be seen from Figure 3 that the amplitude of the sine-type maneuver on horizontal plane becomes smaller with the relative distance becoming smaller. It helps one verify the effectiveness of the proposed method. Accordingly, the HGV will also hit the target successfully, under the ACIGC, however, its trajectory is straighter due to lack of the penetration ability. Figure 4a shows the relative distance of the HGV and the target under the SMIGC and ACIGC. It can be found that the HGV with SMIGC will spend more time hitting the target whereas the landing accuracy can still reach to 0.5 m. Figure 4b indicates that the SMIGC of the HGV will cost more kinetic energy, and therefore it has a lower terminal velocity. Figure 5 describes the time histories of the flight path and flight heading angles of the HGV. It can be seen that the existence of the snake-shape maneuver on the lateral plane will cause periodic changes in the heading angle and the terminal values of the flight path, and that flight heading angles will converge to the desired ones. The time histories of the angle of attack, sideslip angle and bank angle are shown in Figure 6 and Figure 7a. They show that, with the application of the SMIGC and ACIGC, the actual angles can track the expected angles quickly, whereas SMIGC requires larger overload when implementing the penetration, hence, the HGV with SMIGC has larger angle of attack in the early stage of the flight. On the other hand, the sideslip angle remains at the neighborhood of zero while a periodic tile angle is required to provide the lateral overload. Furthermore, during the last stage of the flight, the bank angle will exceed 90° and reach to 120° such that the main lifting surface can provide force to dive. Figure 7b and Figure 8 show the time history of the three-axis angular velocity while Figure 9a describes the time histories of the attitude angles and fin deflection angles, which indicates the effectiveness of the proposed control technique.

4.2. Monte Carlo Simulation

To verify the robustness of the system under the proposed control law, the Monte Carlo simulation was conducted. In this simulation, the perturbation parameters included the aerodynamic parameters, the mass, the moment of inertia, the character size, the character length, the atmosphere density and the wind speed. These perturbation parameters complied with normal distribution $N(\mu ,{\sigma }^2)$, where $\mu$, $\sigma$ represents the mean value and standard deviation of the corresponding perturbation parameters, respectively. The mean value of these perturbation parameters are given as in Appendix A and Appendix B, while the range of these parameters is listed in

Table 3. The range of the perturbation parameters.

Perturbation Parameters	Error Range ($3\mathit{\sigma }$)	Perturbation Parameters	Error Range ($3\mathit{\sigma }$)
Aerodynamic parameters	20%	Mass	2%
Moment of inertia	1%	character size	2%
Character length	1%	atmosphere density	20%
Wind speed	15 m/s

. The simulation results of the 800 shooting tests are described in Figure 10 and Figure 11.

It can be seen from Figure 10a that at each shooting test of the 800 simulations, under the proposed IGC technique, the snake-shape maneuver can all be achieved. From Figure 10b, considering the influence of the perturbation parameters, the landing accuracy, measured by Circular Error Probability (CEP), is better than 1 m. Figure 11 describes the errors of the inclination and deflection of the trajectory. The control accuracy of the trajectory inclination and deflection angles reaches 1°. In summary, according to the simulation results of the 800 shouting tests, the system shows great robustness to the perturbation of the parameters listed in Table 3.

5. Conclusions

Focusing on the large penetration problem of the HGV, this paper proposes a robust IGC method which can help achieve the snake-shape maneuver. Firstly, a novel snake-shape maneuver acceleration command is designed, in which amplitude and phase are determined by height and missile-target distance. The maneuver command will be zero when the height of the HGV is lower than the selected height, indicating that the maneuver will not influence the hitting accuracy. Secondly, based on the missile-target relative motion model and the designed snake-shape maneuver acceleration command, a guidance law, in the framework of sliding mode control, is designed considering the terminal angle constraints of LOS and the snake-shape maneuver. In what follows, the expected angle of attack and bank angle are obtained. Finally, the attitude controller and angular velocity are designed to track expected angle of attack and bank angle, with the application of the adapt control technique to deal with the external disturbance. With the help of the Lyapunov theorem, the global asymptotical stability of the whole system is proved. The effectiveness and robustness of the proposed method are also ensured by feasibility simulation and Monte Carlo simulation.

Compared with some existing works, the proposed method has the following novelties. Firstly, a simple but effective snake-shape maneuver command is designed, which balances the contradictions between the large penetration and the hitting accuracy. Secondly, it deals with the snake-shape maneuver control problem, satisfying the terminal LOS angle constraints simultaneously. Last but not least, the whole system has strong robustness under the proposed control law, with the help of the adaptive and sliding mode control techniques. Although the proposed IGC can only achieve the snake-shape maneuver on the lateral plane instead of any other selected plane, it still shows great potential to achieve the penetration for hypersonic vehicle or other BTT missile.