# Integrated Guidance and Control for Hypersonic Morphing Missile Based on Variable Span Auxiliary Control

Authors: Cunyu Bao; Peng Wang; Guojian Tang
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

The morphing aircraft can improve the flight performance of hypersonic vehicles by satisfying the flight requirements of large airspace and large velocity field. In this paper, for the hypersonic variable span missile, the dynamic model and aerodynamic model are established on the variable span characteristic. The adaptive dynamic surface control back-stepping method is used to establish the integrated guidance and control (IGC) with terminal angular constraint in the dive phase of the hypersonic variable span missile. The span variety is used to assist the lift control to achieve fast and stable control for the centroid motion. The simulation results demonstrate that the feasibility and robustness of the IGC method of the hypersonic variable span missile is better than the invariable span.


---

## Body

## 1. Introduction

In the dive phase, problems such as large Mach number variation range and large overload for a hypersonic vehicle are obvious. The translational and rotational aircraft dynamics are characterized by being fast time-varying, being nonlinear, having a strong coupling, and being uncertain [1]. The traditional aircraft guidance and the control system mainly separate the control and guidance subsystems, without utilizing the coupling information between subsystems based on engineering experience or singular perturbation theory. The integrated guidance control (IGC) system can make full use of the coupling information between the subsystems to improve the performance of the entire system. For the attitude control of hypersonic vehicles, Wang et al. proposed a hierarchical predictive control method [2] and developed a rapid generation method of multitarget entry trajectory [3] that could solve the problem of effective attitude controller and reentry guidance for the general hypersonic vehicle model. Hou and Duan proposed an integrated design method for guidance and control based on an adaptive dynamic surface which avoids the differential explosion problem of the traditional back-stepping method [4]. Wang et al. designed a method of the guidance control system based on banking turning control and dynamic surface control, and the uncertain term was estimated by the state observer and compensated by the control law [5]. Zhao et al. proposed a method of IGC for hypersonic vehicles that satisfies the terminal angular constraint based on the full-coupling model [6]. Then, Zhao et al. studied a method of IGC based on L2 gain interference consistency, which simplifies the structure and calculation of the IGC method [7]. Wang et al. studied the guidance control problem of hypersonic vehicles with composite constraints in the dive phase. The integrated model based on the sliding mode dynamic surface can realize the three-dimensional angle constraint of landing velocity [8].For a hypersonic waverider aircraft, the intelligent morphing technology can effectively solve the problem that aerodynamic performance deteriorates sharply when deviating from the design state [9, 10]. When intelligent morphing technology is applied to the flying missile for the battlefield environment and combat mission changes, the accuracy, combat effectiveness, and cost-effectiveness ratio can greatly improve [11]. In 2004, Jae-Sung et al. studied the aerodynamic and aeroelastic properties of a variable span cruise missile similar to the shape of the Tomahawk and showed that the symmetrical increase in the span of the two wings can effectively reduce the drag and greatly improve the flight range [12]. In the modeling of morphing aircraft, Seigler et al. established a dynamic model of large-scale morphing flight for flight control problems and discussed various control methods [13]. Yang et al. developed dynamic modeling and dynamic response analysis of the variable sweep and variable span aircraft and proved that the combination of a suitable variant mode and variant velocity can reduce the workload of the flight control system [14]. Yue et al. analyzed the longitudinal multibody dynamics of the Z-wing morphing aircraft and showed that the variant of aerodynamic characteristics is the main factor affecting the dynamic characteristics of the aircraft during the folding process [15]. A joint simulation method based on Missile DATCOM and MATLAB proposed by Zhang Fair is applied to the study of variable wing dynamic characteristics of typical axisymmetric missiles [16]. Nobleheart et al. designed a time-varying system for a single-neural network adaptive controller for morphing aircraft [17]. In terms of the control of the morphing aircraft, Dong et al. studied the smooth-switching linear variable parameter (LPV) H∞ control problem for a class of variable span aircraft to ensure the stability and robustness of the aircraft system [18]. Gandhi et al. proposed a control scheme that used model identification to realize morphing control and keep the aircraft stable [19]. In terms of the ballistic design of the morphing missile, Huang et al. used the Radau pseudospectral method (RPM) to optimize the gliding trajectory of the variable sweep aircraft with the target of largest range [20]. Wei-Ming et al. designed a trajectory-attitude dual-loop adaptive sliding mode controller for the variable sweep missile to ensure the stability of the ballistic tracking to the optimal scheme when the wing is scheduled to be actuated [21].Most of the above researches were aimed at the IGC system of the hypersonic vehicle and the stability control of the morphing aircraft. There were few reports on the control research of the hypersonic morphing missile and how to use the span variant to assist the ballistic control. The research in this paper is aimed at the IGC for a class of hypersonic variable span missiles. Based on the adaptive dynamic surface method, an IGC method in the dive phase with terminal angular constraint is designed and the span variant is skillfully applied to the auxiliary control of the IGC to enhance the control accuracy and reduce the workload of the control system. The variable span auxiliary control strategy was hardly studied for hypersonic missiles before. The results of simulation of the new method proposed in the paper show that the total flight time and the miss distance of the variable span are less than invariable span shown. The variety of span could adjust the lift-to-drag ratio of the missile so that it is easier to realize the terminal angular constraint, and the variable span reduces the workload for the attitude control system. At the same time, the robustness of a variable span missile is stronger than that of an invariable span missile.In Section2 in the paper, the geometric model, aerodynamic model, and 6 DOF motion model of the hypersonic deformation missile are established. In Section 3, the IGC model in the dive phase of the hypersonic morphing missile is established and the variable extension control method adapts to the IGC method. Section 4 displays the numerical simulation results and analyses of the IGC method. Section 5 is the conclusion for the paper.


## 2. The Motion Model of Hypersonic Variable Span Missile

### 2.1. The Shape and Parameters of Hypersonic Variable Span Missile

The outline of the shape of the variable-span hypersonic missile selected in this paper is shown in Figure1.Figure 1
Main outline of a hypersonic morphing missile.The shape parameters and morphing modes are shown in Figure2. The lift of the morphing missile was provided by a pair of horizontal front wings, and the four-rear wing was a +-shaped layout. The span variant of the wing is synchronously changed on both sides of the wings without differential changes. The synchronous or differential rotation of the four tails controls the rotational motion of the missile.Figure 2
The shape parameters and morphing modes of the hypersonic morphing missile.(1) Geometric Parameters of the Missile. l1 is the length of the warhead. l2 is the length of the total missile. d is the diameter of the main body and bottom of the warhead(2) Variable Span Geometry. l3 is the distance from the front end of the wing root to the nose of the warhead. c1 is the chord length of the wing root. c2 is the chord length of the wing tip. b1 is the span of the wing. χ1 is the sweep angle of the leading edge of the wing(3) Geometric Parameters of the Tails. l4 is the distance from the front end of the tail root to the nose of warhead. c3 is the chord length of the tail root. c4 is the chord length of the tail tip. b3 is the span of the tail, and χ2 is the sweep angle of the leading edge of the tailb1 is the span of the invariant wing, and b2 is the maximum span of the variable span wing. When b is assumed as the span of the wing in a variable span state, the morphing rate of the wing span is defined as
(1)ξ=b−b1b2−b1.The missile characteristic parameters are shown: in Table1, the mass of the missile is m, and the mass of the wing is m1 and m2 and m2. The moments of inertia are Ix, Iy, and Iz. The reference area is S0. The longitudinal reference length is c¯, and the lateral–directional reference length is b¯.Table 1
Missile characteristic parameter.
Name
Value
Name
Value
Name
Value
l1
3.43 m
c2
0.92 m
χ1
60∘
l2
8.00 m
b1
1.63 m
χ2
51.86∘
l3
4.70 m
b2
2.55 m
S0
5.66 m
l4
9.73 m
c3
1.70 m
c¯
1.148 m
d
1.14 m
c4
0.85 m
b¯
0.766 m
c1
2.12 m
b3
0.91 m
m
5000 kg
mi,i=1,2
100 kg
Ix
15064 kg·m2
Iy,IZ
150646 kg·m2


### 2.2. Hypersonic Variable Span Missile Dynamic Model

#### 2.2.1. Centroid Dynamic Model

Due to the small range of the dive phase, the Earth’s rotation can be ignored. The centroid motion equation of the variable span missile in the ballistic coordinate frame is [22]
(2)av=gHx+−D−F′m,aθ=gHy+LcosγV−NsinγV−Fsy′m,aσ=gHz+LsinγV+NcosγV−Fsz′m.In Eq. (2), a=avaθaσT is the acceleration vector in the ballistic coordinate frame. v is the flight velocity, θ is the flight path angle, and σ is the heading angle. In Eq. (2), the component of gravitational acceleration in the ballistic coordinate frame is [22]
(3)gHxgHygHz=−μr3xcosθcosσ+y+R0sinθ−zcosθsinσ−xsinθcosσ+y+R0cosθ+zsinθsinσxsinσ+zcosσIn Eq. (3), μ is the gravity constant of the Earth, R0 is the radius of the Earth, and r is the mode of the Earth’s radial diameter r, which is expressed in the launch coordinate frame as r=xyzT. In Eq. (3), the expression of Fs′=Fsx′Fsy′Fsz′T is [22]
(4)Fs′=∑i=12Fsi′=HVTBVT∑i=12Fsi.In Eq. (4), HV and BV are the transformation matrixes of the velocity coordinate frame to the ballistic coordinate frame and the body coordinate frame, respectively. Fsi is the additional force due to the span variety and the expression is [22]
(5)Fsi=δ2siδt2+2ω×δsiδt+δωδt×si+ω×ω×si.In Eq. (5), i=1,2 denotes the two sides of the wing, respectively. ω is the angular rates of the missile’s rational movement, and si is the position vector of the wing mass center relative to the missile mass center. The expression of Fsi in the body coordinate frame is [22]
(6)FsxiFsyiFszi=mis¨xi+2ωyṡzi−2ωzṡyi−ωy2+ωz2sxi+ωxωy−ω̇zsyi+ωxωz+ω̇yszis¨yi+2ωzṡxi−2ωxṡzi−ωx2+ωz2syi+ωyωz−ω̇xszi+ωxωy+ω̇zsxis¨zi+2ωxṡyi−2ωyṡxi−ωx2+ωy2szi+ωxωz−ω̇ysxi+ωyωz+ω̇xsyi.The wing is trapezoidal, and the span of the wing varies in theXZ plane of the body coordinate frame. Based on the geometric knowledge, si could be expressed in the body coordinate frame as
(7)sxi=l0−l3−c12+2c2+c16c1+c2c1−c2−2b−dtanχ1,syi=0,szi=2c2+c13c1+c2b+2c1+c26c1+c2d.According to Eq. (7), the first and second derivative of si with time is
(8)ṡxi=−2c2+c1tanχ1ḃ3c1+c2,ṡyi=0,ṡzi=2c2+c13c1+c2ḃ,s¨xi=−2c2+c1tanχ1b¨3c1+c2,s¨yi=0,s¨zi=2c2+c13c1+c2b¨.Each variable component of Eq. (6) is expressed in the body coordinate frame. In Eq. (2), L,D,N are the lift, drag, and lateral force of the missile, respectively, as
(9)L=qS0C¯L,D=qS0C¯D,N=qS0C¯N.In Eq. (9), q=0.5ρv2 is the dynamic pressure. ρ is the atmospheric density which can be calculated according to the standard atmospheric model. C¯L,C¯D,C¯N are the lift coefficient, the drag coefficient, and the lateral force coefficient, respectively. S0 is the reference area when ξ=0. The missile reference area varies during the morphing process. However, for the convenience of calculations, Eq. (11) is regarded as invariant and the influence of the reference area variant caused by the span variant is classified into the equivalent aerodynamic coefficient as
(10)C¯jξ=SξCjξS0,j=L,D,N.The equivalent aerodynamic coefficient calculation model is in Section2.3.


#### 2.2.2. Rotational Dynamic Model

According to the literature, the rotational motion model of the variable span missile can be(11)ω̇x=Ix−1Mx+MSx+Ix−1Iy−Izωyωz−Ix−1İxωx,ω̇y=Iy−1My+MSy+Iz−1Iz−Ixωzωx−Iy−1İyωy,ω̇z=Iz−1Mz+MSz+Iz−1Ix−Iyωxωy−Iz−1İzωz.In Eq. (11), the additional moment due to the span variant is expressed as [22]
(12)MS=∑i=12si×migi−∑i=12misi×dvodt+misi×d2sidt2.In Eq. (12), g is the gravity vector and vo is the velocity vector of the missile centroid. The component expressions in the body coordinate frame are [22]
(13)MS=−∑i=12miμr3sxisyiszi×BGxy+Rez+misxisyiszi×v̇xv̇yv̇z+mi−ωzsyi+ωysziωzsxi−ωxszi−ωysxi+ωxsyi×vxvyvz+sxisyiszi×Fsxi′Fsyi′Fszi′.In Eq. (13), BG is the transformation matrix of the launch coordinate frame to the body coordinate frame, and the other components of each variable are expressed in the body coordinate frame.In Eq. (11), M=MxMyMzT are the roll, yaw, and pitch moment, respectively, of the missile as
(14)Mx=qS0bm¯x,My=qS0bm¯y,Mz=qS0cm¯z.In (14), mx,my,mz are the roll, yaw, and pitch moment coefficients, respectively. Similar to Eq. (11), the variant in the reference area caused by the span variant is converted into the equivalent aerodynamic torque coefficient as
(15)m¯jξ=SξmjξS0,j=x,y,z.The equivalent aerodynamic moment coefficient calculation model is in Section2.3.


### 2.3. Aerodynamic Model of Hypersonic Variable Span Missile

Due to the lack of models and aerodynamic research of hypersonic variable span missiles, the mature hypersonic morphing aircraft is unavailable. This paper uses the Missile DATCOM software of the US Air Force Laboratory (AFL) to estimate the variable span missile aerodynamic model. The software is used to calculate the aerodynamic data under different combinations of morphing rate, Mach number, attack angle, side slip angle, and roll, yaw, and pitch fin deflection angles. Then, the model is identified by the least squares method. The aerodynamic model of the variable span missile is(16)C¯L=CL,Ma,ξα⋅αL+CL,Ma,ξ,αδx⋅δxL+CL,Ma,ξ,αδx⋅δzL,C¯D=CD,Ma,ξ0+CD,Ma,ξα⋅αD+CD,Ma,ξ,αδx⋅δxD+CL,Ma,ξ,αδy⋅δyD+CL,Ma,ξ,αδz⋅δzD,C¯N=CN,Ma,αβ⋅βN+CL,Ma,αδy⋅δyN,m¯x=mx,Ma,αβ⋅βmx+mx,Ma,αδx⋅δxmx,m¯y=my,Ma,αβ⋅βmy+my,Ma,αδy⋅δymy,m¯z=mz,Ma,ξα⋅αmz+mz,Ma,ξ,αδz⋅δzmz.In Eq. (16), the subscript of each coefficient represents the variable contained in the coefficient matrix. According to the model identified by the obtained aerodynamic data, the structural form of each variable in Eq. (17) is
(17)αL=1αα3T,δxL=δx2,δzL=δzδz2T,αD=1α2T,δxD=δx2,δyD=δy2,δxD=δz,βN=β,δyN=δy,βmx=β,δxmx=δx,βmy=β,δymy=δy,αmz=1αα3T,δzmz=δzδz2T,CP,Ma,ξq=CP,Maq⋅1ξT.In Eq. (17), CP,Ma,ξq is the coefficient matrix of different aerodynamic coefficients (P=L,D,N,mx,my,mz) relative to different variables (q=α,β,δx,δy,δz). CP,Maq is the coefficient matrix of CP,Ma,ξq relative to the morphing rate. For the designing of the IGC method, the aerodynamic load model needs to be processed as
(18)L=qS0CLα⋅α+wL,Mx=qS0bmxβ⋅βmx+mxδx⋅δx+wMx,My=qS0bmyβ⋅βmy+myδy⋅δy+wMy,Mz=qS0cmzαα+mzδz⋅δz+wMz.In Eq. (18), α=1αα3T. CLα is the partial derivative vector of the lift coefficient for the attack angle. The uncertainty wL includes the influence of pitch fin deflections on the lift. mxδx,myδy, and mzδz are the partial derivatives of the roll moment coefficient for the roll fin deflections, the yaw moment coefficient for the yaw fin deflections, and the pitch moment coefficient for the pitch fin deflections, respectively. The uncertainty wMz includes the influence of the quadratic pitch fin deflections on the pitching moment.


## 2.1. The Shape and Parameters of Hypersonic Variable Span Missile

The outline of the shape of the variable-span hypersonic missile selected in this paper is shown in Figure1.Figure 1
Main outline of a hypersonic morphing missile.The shape parameters and morphing modes are shown in Figure2. The lift of the morphing missile was provided by a pair of horizontal front wings, and the four-rear wing was a +-shaped layout. The span variant of the wing is synchronously changed on both sides of the wings without differential changes. The synchronous or differential rotation of the four tails controls the rotational motion of the missile.Figure 2
The shape parameters and morphing modes of the hypersonic morphing missile.(1) Geometric Parameters of the Missile. l1 is the length of the warhead. l2 is the length of the total missile. d is the diameter of the main body and bottom of the warhead(2) Variable Span Geometry. l3 is the distance from the front end of the wing root to the nose of the warhead. c1 is the chord length of the wing root. c2 is the chord length of the wing tip. b1 is the span of the wing. χ1 is the sweep angle of the leading edge of the wing(3) Geometric Parameters of the Tails. l4 is the distance from the front end of the tail root to the nose of warhead. c3 is the chord length of the tail root. c4 is the chord length of the tail tip. b3 is the span of the tail, and χ2 is the sweep angle of the leading edge of the tailb1 is the span of the invariant wing, and b2 is the maximum span of the variable span wing. When b is assumed as the span of the wing in a variable span state, the morphing rate of the wing span is defined as
(1)ξ=b−b1b2−b1.The missile characteristic parameters are shown: in Table1, the mass of the missile is m, and the mass of the wing is m1 and m2 and m2. The moments of inertia are Ix, Iy, and Iz. The reference area is S0. The longitudinal reference length is c¯, and the lateral–directional reference length is b¯.Table 1
Missile characteristic parameter.
Name
Value
Name
Value
Name
Value
l1
3.43 m
c2
0.92 m
χ1
60∘
l2
8.00 m
b1
1.63 m
χ2
51.86∘
l3
4.70 m
b2
2.55 m
S0
5.66 m
l4
9.73 m
c3
1.70 m
c¯
1.148 m
d
1.14 m
c4
0.85 m
b¯
0.766 m
c1
2.12 m
b3
0.91 m
m
5000 kg
mi,i=1,2
100 kg
Ix
15064 kg·m2
Iy,IZ
150646 kg·m2


## 2.2. Hypersonic Variable Span Missile Dynamic Model

### 2.2.1. Centroid Dynamic Model

Due to the small range of the dive phase, the Earth’s rotation can be ignored. The centroid motion equation of the variable span missile in the ballistic coordinate frame is [22]
(2)av=gHx+−D−F′m,aθ=gHy+LcosγV−NsinγV−Fsy′m,aσ=gHz+LsinγV+NcosγV−Fsz′m.In Eq. (2), a=avaθaσT is the acceleration vector in the ballistic coordinate frame. v is the flight velocity, θ is the flight path angle, and σ is the heading angle. In Eq. (2), the component of gravitational acceleration in the ballistic coordinate frame is [22]
(3)gHxgHygHz=−μr3xcosθcosσ+y+R0sinθ−zcosθsinσ−xsinθcosσ+y+R0cosθ+zsinθsinσxsinσ+zcosσIn Eq. (3), μ is the gravity constant of the Earth, R0 is the radius of the Earth, and r is the mode of the Earth’s radial diameter r, which is expressed in the launch coordinate frame as r=xyzT. In Eq. (3), the expression of Fs′=Fsx′Fsy′Fsz′T is [22]
(4)Fs′=∑i=12Fsi′=HVTBVT∑i=12Fsi.In Eq. (4), HV and BV are the transformation matrixes of the velocity coordinate frame to the ballistic coordinate frame and the body coordinate frame, respectively. Fsi is the additional force due to the span variety and the expression is [22]
(5)Fsi=δ2siδt2+2ω×δsiδt+δωδt×si+ω×ω×si.In Eq. (5), i=1,2 denotes the two sides of the wing, respectively. ω is the angular rates of the missile’s rational movement, and si is the position vector of the wing mass center relative to the missile mass center. The expression of Fsi in the body coordinate frame is [22]
(6)FsxiFsyiFszi=mis¨xi+2ωyṡzi−2ωzṡyi−ωy2+ωz2sxi+ωxωy−ω̇zsyi+ωxωz+ω̇yszis¨yi+2ωzṡxi−2ωxṡzi−ωx2+ωz2syi+ωyωz−ω̇xszi+ωxωy+ω̇zsxis¨zi+2ωxṡyi−2ωyṡxi−ωx2+ωy2szi+ωxωz−ω̇ysxi+ωyωz+ω̇xsyi.The wing is trapezoidal, and the span of the wing varies in theXZ plane of the body coordinate frame. Based on the geometric knowledge, si could be expressed in the body coordinate frame as
(7)sxi=l0−l3−c12+2c2+c16c1+c2c1−c2−2b−dtanχ1,syi=0,szi=2c2+c13c1+c2b+2c1+c26c1+c2d.According to Eq. (7), the first and second derivative of si with time is
(8)ṡxi=−2c2+c1tanχ1ḃ3c1+c2,ṡyi=0,ṡzi=2c2+c13c1+c2ḃ,s¨xi=−2c2+c1tanχ1b¨3c1+c2,s¨yi=0,s¨zi=2c2+c13c1+c2b¨.Each variable component of Eq. (6) is expressed in the body coordinate frame. In Eq. (2), L,D,N are the lift, drag, and lateral force of the missile, respectively, as
(9)L=qS0C¯L,D=qS0C¯D,N=qS0C¯N.In Eq. (9), q=0.5ρv2 is the dynamic pressure. ρ is the atmospheric density which can be calculated according to the standard atmospheric model. C¯L,C¯D,C¯N are the lift coefficient, the drag coefficient, and the lateral force coefficient, respectively. S0 is the reference area when ξ=0. The missile reference area varies during the morphing process. However, for the convenience of calculations, Eq. (11) is regarded as invariant and the influence of the reference area variant caused by the span variant is classified into the equivalent aerodynamic coefficient as
(10)C¯jξ=SξCjξS0,j=L,D,N.The equivalent aerodynamic coefficient calculation model is in Section2.3.


### 2.2.2. Rotational Dynamic Model

According to the literature, the rotational motion model of the variable span missile can be(11)ω̇x=Ix−1Mx+MSx+Ix−1Iy−Izωyωz−Ix−1İxωx,ω̇y=Iy−1My+MSy+Iz−1Iz−Ixωzωx−Iy−1İyωy,ω̇z=Iz−1Mz+MSz+Iz−1Ix−Iyωxωy−Iz−1İzωz.In Eq. (11), the additional moment due to the span variant is expressed as [22]
(12)MS=∑i=12si×migi−∑i=12misi×dvodt+misi×d2sidt2.In Eq. (12), g is the gravity vector and vo is the velocity vector of the missile centroid. The component expressions in the body coordinate frame are [22]
(13)MS=−∑i=12miμr3sxisyiszi×BGxy+Rez+misxisyiszi×v̇xv̇yv̇z+mi−ωzsyi+ωysziωzsxi−ωxszi−ωysxi+ωxsyi×vxvyvz+sxisyiszi×Fsxi′Fsyi′Fszi′.In Eq. (13), BG is the transformation matrix of the launch coordinate frame to the body coordinate frame, and the other components of each variable are expressed in the body coordinate frame.In Eq. (11), M=MxMyMzT are the roll, yaw, and pitch moment, respectively, of the missile as
(14)Mx=qS0bm¯x,My=qS0bm¯y,Mz=qS0cm¯z.In (14), mx,my,mz are the roll, yaw, and pitch moment coefficients, respectively. Similar to Eq. (11), the variant in the reference area caused by the span variant is converted into the equivalent aerodynamic torque coefficient as
(15)m¯jξ=SξmjξS0,j=x,y,z.The equivalent aerodynamic moment coefficient calculation model is in Section2.3.


## 2.2.1. Centroid Dynamic Model

Due to the small range of the dive phase, the Earth’s rotation can be ignored. The centroid motion equation of the variable span missile in the ballistic coordinate frame is [22]
(2)av=gHx+−D−F′m,aθ=gHy+LcosγV−NsinγV−Fsy′m,aσ=gHz+LsinγV+NcosγV−Fsz′m.In Eq. (2), a=avaθaσT is the acceleration vector in the ballistic coordinate frame. v is the flight velocity, θ is the flight path angle, and σ is the heading angle. In Eq. (2), the component of gravitational acceleration in the ballistic coordinate frame is [22]
(3)gHxgHygHz=−μr3xcosθcosσ+y+R0sinθ−zcosθsinσ−xsinθcosσ+y+R0cosθ+zsinθsinσxsinσ+zcosσIn Eq. (3), μ is the gravity constant of the Earth, R0 is the radius of the Earth, and r is the mode of the Earth’s radial diameter r, which is expressed in the launch coordinate frame as r=xyzT. In Eq. (3), the expression of Fs′=Fsx′Fsy′Fsz′T is [22]
(4)Fs′=∑i=12Fsi′=HVTBVT∑i=12Fsi.In Eq. (4), HV and BV are the transformation matrixes of the velocity coordinate frame to the ballistic coordinate frame and the body coordinate frame, respectively. Fsi is the additional force due to the span variety and the expression is [22]
(5)Fsi=δ2siδt2+2ω×δsiδt+δωδt×si+ω×ω×si.In Eq. (5), i=1,2 denotes the two sides of the wing, respectively. ω is the angular rates of the missile’s rational movement, and si is the position vector of the wing mass center relative to the missile mass center. The expression of Fsi in the body coordinate frame is [22]
(6)FsxiFsyiFszi=mis¨xi+2ωyṡzi−2ωzṡyi−ωy2+ωz2sxi+ωxωy−ω̇zsyi+ωxωz+ω̇yszis¨yi+2ωzṡxi−2ωxṡzi−ωx2+ωz2syi+ωyωz−ω̇xszi+ωxωy+ω̇zsxis¨zi+2ωxṡyi−2ωyṡxi−ωx2+ωy2szi+ωxωz−ω̇ysxi+ωyωz+ω̇xsyi.The wing is trapezoidal, and the span of the wing varies in theXZ plane of the body coordinate frame. Based on the geometric knowledge, si could be expressed in the body coordinate frame as
(7)sxi=l0−l3−c12+2c2+c16c1+c2c1−c2−2b−dtanχ1,syi=0,szi=2c2+c13c1+c2b+2c1+c26c1+c2d.According to Eq. (7), the first and second derivative of si with time is
(8)ṡxi=−2c2+c1tanχ1ḃ3c1+c2,ṡyi=0,ṡzi=2c2+c13c1+c2ḃ,s¨xi=−2c2+c1tanχ1b¨3c1+c2,s¨yi=0,s¨zi=2c2+c13c1+c2b¨.Each variable component of Eq. (6) is expressed in the body coordinate frame. In Eq. (2), L,D,N are the lift, drag, and lateral force of the missile, respectively, as
(9)L=qS0C¯L,D=qS0C¯D,N=qS0C¯N.In Eq. (9), q=0.5ρv2 is the dynamic pressure. ρ is the atmospheric density which can be calculated according to the standard atmospheric model. C¯L,C¯D,C¯N are the lift coefficient, the drag coefficient, and the lateral force coefficient, respectively. S0 is the reference area when ξ=0. The missile reference area varies during the morphing process. However, for the convenience of calculations, Eq. (11) is regarded as invariant and the influence of the reference area variant caused by the span variant is classified into the equivalent aerodynamic coefficient as
(10)C¯jξ=SξCjξS0,j=L,D,N.The equivalent aerodynamic coefficient calculation model is in Section2.3.


## 2.2.2. Rotational Dynamic Model

According to the literature, the rotational motion model of the variable span missile can be(11)ω̇x=Ix−1Mx+MSx+Ix−1Iy−Izωyωz−Ix−1İxωx,ω̇y=Iy−1My+MSy+Iz−1Iz−Ixωzωx−Iy−1İyωy,ω̇z=Iz−1Mz+MSz+Iz−1Ix−Iyωxωy−Iz−1İzωz.In Eq. (11), the additional moment due to the span variant is expressed as [22]
(12)MS=∑i=12si×migi−∑i=12misi×dvodt+misi×d2sidt2.In Eq. (12), g is the gravity vector and vo is the velocity vector of the missile centroid. The component expressions in the body coordinate frame are [22]
(13)MS=−∑i=12miμr3sxisyiszi×BGxy+Rez+misxisyiszi×v̇xv̇yv̇z+mi−ωzsyi+ωysziωzsxi−ωxszi−ωysxi+ωxsyi×vxvyvz+sxisyiszi×Fsxi′Fsyi′Fszi′.In Eq. (13), BG is the transformation matrix of the launch coordinate frame to the body coordinate frame, and the other components of each variable are expressed in the body coordinate frame.In Eq. (11), M=MxMyMzT are the roll, yaw, and pitch moment, respectively, of the missile as
(14)Mx=qS0bm¯x,My=qS0bm¯y,Mz=qS0cm¯z.In (14), mx,my,mz are the roll, yaw, and pitch moment coefficients, respectively. Similar to Eq. (11), the variant in the reference area caused by the span variant is converted into the equivalent aerodynamic torque coefficient as
(15)m¯jξ=SξmjξS0,j=x,y,z.The equivalent aerodynamic moment coefficient calculation model is in Section2.3.


## 2.3. Aerodynamic Model of Hypersonic Variable Span Missile

Due to the lack of models and aerodynamic research of hypersonic variable span missiles, the mature hypersonic morphing aircraft is unavailable. This paper uses the Missile DATCOM software of the US Air Force Laboratory (AFL) to estimate the variable span missile aerodynamic model. The software is used to calculate the aerodynamic data under different combinations of morphing rate, Mach number, attack angle, side slip angle, and roll, yaw, and pitch fin deflection angles. Then, the model is identified by the least squares method. The aerodynamic model of the variable span missile is(16)C¯L=CL,Ma,ξα⋅αL+CL,Ma,ξ,αδx⋅δxL+CL,Ma,ξ,αδx⋅δzL,C¯D=CD,Ma,ξ0+CD,Ma,ξα⋅αD+CD,Ma,ξ,αδx⋅δxD+CL,Ma,ξ,αδy⋅δyD+CL,Ma,ξ,αδz⋅δzD,C¯N=CN,Ma,αβ⋅βN+CL,Ma,αδy⋅δyN,m¯x=mx,Ma,αβ⋅βmx+mx,Ma,αδx⋅δxmx,m¯y=my,Ma,αβ⋅βmy+my,Ma,αδy⋅δymy,m¯z=mz,Ma,ξα⋅αmz+mz,Ma,ξ,αδz⋅δzmz.In Eq. (16), the subscript of each coefficient represents the variable contained in the coefficient matrix. According to the model identified by the obtained aerodynamic data, the structural form of each variable in Eq. (17) is
(17)αL=1αα3T,δxL=δx2,δzL=δzδz2T,αD=1α2T,δxD=δx2,δyD=δy2,δxD=δz,βN=β,δyN=δy,βmx=β,δxmx=δx,βmy=β,δymy=δy,αmz=1αα3T,δzmz=δzδz2T,CP,Ma,ξq=CP,Maq⋅1ξT.In Eq. (17), CP,Ma,ξq is the coefficient matrix of different aerodynamic coefficients (P=L,D,N,mx,my,mz) relative to different variables (q=α,β,δx,δy,δz). CP,Maq is the coefficient matrix of CP,Ma,ξq relative to the morphing rate. For the designing of the IGC method, the aerodynamic load model needs to be processed as
(18)L=qS0CLα⋅α+wL,Mx=qS0bmxβ⋅βmx+mxδx⋅δx+wMx,My=qS0bmyβ⋅βmy+myδy⋅δy+wMy,Mz=qS0cmzαα+mzδz⋅δz+wMz.In Eq. (18), α=1αα3T. CLα is the partial derivative vector of the lift coefficient for the attack angle. The uncertainty wL includes the influence of pitch fin deflections on the lift. mxδx,myδy, and mzδz are the partial derivatives of the roll moment coefficient for the roll fin deflections, the yaw moment coefficient for the yaw fin deflections, and the pitch moment coefficient for the pitch fin deflections, respectively. The uncertainty wMz includes the influence of the quadratic pitch fin deflections on the pitching moment.


## 3. IGC Method in the Dive Phase with Variable Span Auxiliary Control

### 3.1. IGC Model in the Dive Phase

In this paper, the IGC design of the dive phase of the hypersonic variable span missile is studied. The fixed point on the ground is selected as the target, and the IGC model in the dive phase with the terminal angular constraint is established.


#### 3.1.1. Relative Motion of Variable Span Missile and Target

As shown in Figure3, the East–North–Up (ENU) coordinate frame is formed by a plane tangential to the Earth’s surface and is attached to the target. The East axis is labeled X, the North axis is labeled Y, and the up axis is labeled Z. Define the line-of-sight (LOS) coordinate frame so that its origin is on the target where the Sx-axis points to the vehicle, the Sy-axis is in the horizontal plane of the target, and the Sz-axis is determined using the right-hand rule. ηD is the angle between the direction of velocity and the LOS, and γD is the longitudinal velocity azimuth in the dive plane. rT is the distance between the center of mass of the variable span missile and the target point. λD and λT are the angles of elevation and azimuth of LOS.Figure 3
ENU coordinate frame and LOS coordinate frame.The relative dynamic model of the variable span missile and the target point by the line-of-sight angles are expressed as [5]
(19)r¨T=rTλ̇D2+rTλ̇T2cos2λD+arT,λ¨D=−2ṙTλ̇DrT−λ̇T2sinλDcosλD+aλDrT,λ¨T=−2ṙTλ̇TrT+2λ̇Dλ̇TtanλD−aλTrTcosλD.The relationship between the components in the LOS coordinate frame and the trajectory coordinate frame of the missile acceleration vector can be expressed as [6]
(20)aλD=SH2,1aV+SH2,2aθ+SH2,3aσ,aλT=SH3,1aV+SH3,2aθ+SH3,3aσ,where SHi,ji,j=1,2,3are the elements in the conversion matrix of the trajectory coordinate frame to the LOS coordinate frame. i represents the row, and j represents the column.To study the effect of the span variant on the IGC method, the hypersonic missile adopts the BTT control strategy in the dive phase to show the influence of the span variant more clearly. Therefore, the lateral force can be neglected. The influence of the lateral force in Equation (2) is regarded as the uncertainty term. Substituting Eq. (18) into Eq. (2) yields [6]
(21)aV=gHx−D+Fsx′m,aθ=gHy−1mFsy′+1mqS0CLcosγV+wθ,aσ=gHz−1mFsz′+1mqS0CLsinγV+wσ,where the uncertain items wθ and wσ can be expressed as [6]
(22)wθ=cosγVmwL−sinγVmN,wσ=sinγVmwL+cosγVmN.For missiles and fixed targets on the ground, the relative motion between them can be expressed as [6]
(23)ṙT=−vcosηD,rTλ̇D=vsinηD,then
(24)rTλ̇D=vsinηD=vsinλD+γDF=0,where γDF is the terminal angle. Because of λD+γDF<π/2 at the terminal time, it can be obtained as
(25)λD+γDF=0.The terminal angle can be constrained by Eq. (25). By substituting Eq. (20) and Eq. (21) into Eq. (19), the relative dynamics of the target and the hypersonic variable span missile are denoted as
(26)ẋF=x01,ẋ0=f0+g0x1∗+w0.In Eq. (26),
(27)xF=λD+γDF,x0=λ̇Dλ̇TT=x01x02T,x1∗=CLcosγVCLsinγVT,w0=wλ¨Dwλ¨TT,f0=−2ṙTλ̇DrT−λ̇T2sinλDcosλD−2ṙTλ̇TrT+2λ̇Dλ̇TtanλD+1rTSH2,1SH2,2SH2,3−SH3,1cosλD−SH3,2cosλD−SH3,3cosλDaVgHy−1mFsy′gHz−1mFsz′,g0=qS0mrTSH2,2SH2,3−SH3,2cosλD−SH3,3cosλD,w0=1rTSH2,2SH2,3−SH3,2cosλD−SH3,3cosλDwθwσ.


#### 3.1.2. Rotational Motion Model

The rotational motion model by the attack, sideslip, and bank angle is [7]
(28)α̇β̇γ̇V=−cosαtanβsinαtanβ1sinαcosα0cosαsecβ−sinαsecβ0ωxωyωz−secβcosγVsinγV−tanβcosγVaθV+−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVaθVcosθ.Then, the rotational motion model that is available for control can be denoted as(29)ẋ1=f1+g1x2+w1.In Eq. (29),
(30)x1=αβγVT,x2=ωxωyωzT,w1=wα̇wβ̇wγ̇VT,(31)f1=−secβcosγVsinγV−tanβcosγVaθV+−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVaσVcosθ,g1=−cosαtanβsinαtanβ1sinαcosα0cosαsecβ−sinαsecβ0,w1=−secβcosγVsinγV−tanβcosγVwθ̇−−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVwσ̇.In Eq. (31), the uncertainty can be expressed as [5]
(32)wθ̇=wθv,wσ̇=−wσvcosθ.


#### 3.1.3. Processing of the Rotational Dynamic Model

When ignoring the Earth’s rotation, the launch coordinate frame is equivalent to the ground inertial system. By substituting Eq. (18) into Eq. (11), the simplified state equations with the input of fin deflections can be denoted as [7]
(33)ω̇x=qS0bmxδxIxδx+qS0bmxββ+MSx−Iz−IyIxωyωz+wω̇x,ω̇y=qS0bmyδxIyδy+qS0bmyββ+MSy−Ix−IzIyωxωz+wω̇y,ω̇z=qS0cmzδzIzδz+qS0cmzαα+MSz−Iy−IxIzωxωy+wω̇z.Then, the rotational dynamic model that is available for control can be denoted as(34)ẋ2=f2x1,x2+g2tu+w¯2.In Eq. (34) [5],
(35)u=δaδeδrT,w2=wω̇xwω̇ywω̇zT,(36)f2=qS0bmxββ+MSx−Iz−IyIxωyωzqS0bmyββ+MSy−Ix−IzIyωxωzqS0cmzαα+MSz−Iy−IxIzωxωy,g2=qS0⋅diagbmxδxIx,bmyδyIy,cmzδzIz,w2=−İxωxIx−İyωyIywMz−İzωzIz.The variant rate of the three-axis moment of inertia caused by the variable span is difficult to calculate accurately, so it is regarded as an uncertainty in the Eq. (36).By the combination of Eq. (26), (29), and (34), an IGC model for hypersonic variable span missiles with terminal angular constraints can be obtained as
(37)ẋF=x01,ẋ0=f0+g0x1∗+w0,ẋ1=f1+g1x2+w1,ẋ2=f2+g2u+w2,y=xFx0TT.In Eq. (37), w0, w1, and w2 are bounded function vectors for uncertain positions. Let us assume that there is a set of unknown constants ei that satisfies [7]
(38)wi≤ei,i=0,1,2.The method of solvedx1 by x1∗ will be introduced in Section 2.2.


### 3.2. The IGC Method in the Dive Phase with Variable Span Control

The IGC model in the dive phase guidance control with terminal angular constraint is shown in Eq. (37). The adaptive dynamic surface back-stepping control method is designed to make the system stable and make the output as close as possible to zero and make the system have strong robustness for uncertain factors.


#### 3.2.1. Design of the Adaptive Dynamic Surface Back-Stepping IGC Method

To design an adaptive dynamic surface back-stepping control method, Eq. (37) is required to nonsingular [23] as
(39)detg1=αsecβ≠0.Therefore, let us assume to makeg1 nonsingularity as
(40)α,β∈Ω:α,β0<α<π2,β<π2,in the whole dive phase. And Ω is the bounded closed set in R2, so Ω is a tight set.The design steps of the adaptive dynamic surface back-stepping IGC method are as follows [6].Step 1.
Design the dynamic surface [5]
(41)s0=λ̇D+kFvrxFλ̇TT=s01s02T,where kF is the angular error coefficient, of which the magnitude determines the weight of the terminal angular error in the dynamic surface.
The rate of change of the dynamic surfaces0 is designed according to the exponential approach law.
(42)ṡ0=−vk01rs01−ε01rsats01,d01−vk02rs02−ε02rsats02,d02,where k01 and k02 are positive gain constants, ε01 andε02 are the gain-of-saturation function, and d01 and d02 are the boundary layer thickness. sats,d is a saturation function which defined as [7]
(43)sats,d=1,s>d,sd,s≤d,−1,s<−d.
When the missile is far from the target point andr is large, the rate of approaching law in (42) will be slow so that the missile has a small overload in the initial stage. When the missile approaches the target point and r is small, the rate of approaching law in (42) increases so that λ̇D does not diverge to improve the accuracy of hitting. The selection of the saturation function can effectively eliminate the chattering of the first virtual control variable.
Deriving the dynamic surfaces0 and then combining (42) can get the first virtual control input x1v as
(44)x1v∗=g0−1−vrk0s0−1rε0sats0,d0−f0′−f0.
In Eq. (44), the result of deriving the terminal angle is
(45)f0′=kFvrx01+kFv̇r−ṙvr2xF0.
In Eq. (44), k0=diagk01,k02, ε0=diagε01,ε02, and d0=d01d02T.
According tox1v∗=CLvcosγVvCLvsinγVvT, it can be translated as
(46)γVv=tan−1x1v∗2x1v∗1,CLv=x1v∗1cosγVv.
According to the aerodynamic model of the variable span missile, the structural form of the lift coefficient caused by the attack angle is(47)CLv=CL0+CLααv+CLα3αv3.
It can be solved as(48)αv=CLv−CL0−CLα3αv,t−13CLα.
In Eq. (48), take αv as the main factor affecting for lift, and αv,t−13 is the previous moment attack angle.
Then,(49)x1v=αvγVvT.
In order to exert the advantage of the variable span hypersonic missile, the morphing rate and attack angle of the wing are coordinated controls. The roadmap of the IGC method for a hypersonic morphing missile based on variable span auxiliary control is shown in Figure4, while the IGC method with invariable span is shown in Figure 5.
The principle of coordinated control shown in Figure4 is to translate the change of the lift coefficient into a change of the attack angle and the morphing rate, so that the span variant can withstand the change demand of the lift command and reduce the change demand of the attack angle control. While in the traditional IGC method with invariable span, the change demand of the lift command is only reached by the attack angle. It would improve the stability and velocity of missile IGC mission. That is,
(50)ĊLv=ĊLαξαv+CLαξα̇v.
Solve Eq. (50) to get
(51)ĊLαξ=ĊLv−CLαξt−1α̇v−ĊL0,t−1ξt−1−ĊL,t−1α3ξt−1αv3αv.
In Eq. (51), take αv as the main factor affecting for lift; ξt−1,ĊL0,t−1,ĊL,t−1α3 is the previous moment value of each variable.
According to the model of the aerodynamic coefficient of the missile,(52)CLαξ=1ξ⋅CLα,2×3ξ,where CLαξ is the coefficient matrix of CLα with respect to the morphing rate, and the expression is
(53)CLαξ=CLα,1,1ξCLα,1,2ξCLα,1,3ξCLα,2,1ξCLα2,2ξCLα,2,3ξ.
Let us substitute Eq. (53) into (3.35) and deriving on both sides to get
(54)ĊLαξ=ξ̇⋅CLα2,2ξ.
Combining Eq. (51) and Eq. (54), the variation law of the morphing rate can be obtained as
(55)ξ̇c=ĊLv−CLαξt−1α̇v−ĊL0,t−1ξt−1−ĊL,t−1α3ξt−1αv3αv⋅CLα2,2ξ.
In order to prevent differential explosion, it will pass through the first-order filter as(56)τξξ¨+ξ̇=ξ̇cξ̇0=ξ̇c0.Step 2.
Design the dynamic surface(57)s1=x1−x1d∗,where x1d∗=x1d10x1d2T, x1d is the output value of the following first-order filter as
(58)τ1ẋ1d+x1d=x1vx1d0=x1v0,where τ1=diagτ11,0,τ13 is the time constant of the filter.
Design the second virtual control input as(59)x2v=−g1−1f1+e^1s1+k1s1−ẋ1d∗+g0Ts010g0Ts02T,where ẋ1d∗=ẋ1d10ẋ1d2T, k1=diagk11,k12,k13, and e^1 is an estimate of the upper bound e1 of the uncertainty w1, which can be obtained by the following adaptive law
(60)e^̇1=υ1s1Ts1−μ1e^1,e^10=0,where υ1 and μ1 are positive constants.Step 3.
Design the dynamic surface(61)s2=x2−x2d.
Similarly,x2d is the output value of the following first-order filter as
(62)τ2ẋ2d+x2d=x2vx2d0=x2v0,where τ2=diagτ21,τ22,τ23 is the time constant of the filter. The control input of the system is designed as
(63)u=−g2−1f2+e^2s2+k2s2−ẋ2d+g1Ts1,where k2=diagk21,k22,k23 and e^2 is an estimate of the upper bound e2 of the uncertainty w2, which can be obtained by the following adaptive law as
(64)e^̇2=υ2s2Ts2−μ2e^2,e^20=0,where υ2 and μ2 are positive constants.
In summary, the adaptive dynamic surface back-stepping IGC method for variable span control of the system (3.20) can be expressed as(65)s0=x01+kFvRxFx02T,x1v=g0−1−vrk0s0−1rε0sats0,d0−f0′−f0,ξ̇c=ĊLv−CLαξα̇vαvTαvαvT−1CLα,2ξTCLα,2ξCLα,2ξT−1,τξξ¨+ξ̇=ξ̇cξ̇0=ξ̇c0,τ1ẋ1d+x1d=x1v,x1d0=x1v0,x1d∗=x1d10x1d2T,ẋ1d∗=ẋ1d10ẋ1d2T,s1=x1−x1d∗,x2v=−g1−1f1+e^1s1+k1s1−ẋ1d∗+g0Ts010g0Ts02T,e^̇1=υ1s1Ts1−μ1e^1,e^10=0,τ2ẋ2d+x2d=x2v,x2d0=x2v0,s2=x2−x2d,u=−g2−1f2+e^2s2+k2s2−ẋ2d+g1Ts1,e^̇2=υ2s2Ts2−μ2e^2,e^20=0.Figure 4
The roadmap of the IGC method with variable span auxiliary control.Figure 5
The roadmap of the IGC method with invariable span.


#### 3.2.2. Stability Analysis

The stability of the closed-loop system obtained using the control law in Eq. (65) is analyzed below.Define the filter errorb1 and b2, respectively, as
(66)bi=xid−xiv,i=1,2.Then, the derivative of the filtering error is(67)ḃi=τi−1xiv−xid−ẋiv=−τi−1bi−ẋiv,i=1,2.AssumingBi=−ẋiv,i=1,2, according to the assumption in Section 3.2.1, Bi is a tight set, then it has an upper bound as [7]
(68)Bi≤Mi,i=1,2.We define the upper bound estimation error of the norm as [7]
(69)e~i=ei−e^i,i=1,2.Then, its derivative is(70)e~̇i=−e^̇i,i=1,2.Substituting the parameter adaptive law (60) and Eq. (64) into Eq. (70) yields
(71)e~̇i=−υisiTsi−μie^i=−υisiTsi−μiei+μie~i,i=1,2.s1∗ can be defined as
(72)s1∗=x1∗−x1d=x1∗−b1+x1v.The derivative ofs0 is
(73)ṡ0=f0+g0s1∗+b1+x1v+w0+f0′.Substituting Eq. (44) into Eq. (73) yields
(74)ṡ0=g0s1∗+b1+w0−vrk0s0−1rε0sats0,d0.The derivative ofs1 is
(75)ṡ1=f1+g1s2+x2v+b2+w1−ẋ1d∗.Substituting Eq. (59) into Eq. (75) yields
(76)ṡ1=g1s2+b2+w1−k1s1−e^1s1−g0Ts010g0Ts02T.The derivative ofs2 is
(77)ṡ2=f2+g2u+w2−ẋ2d.Now the expression of the control inputu can be obtained by substituting Eq. (77) into Eq. (63) as
(78)ṡ2=w2−k2s2−e^2s2−g1Ts1In summary, the entire closed-loop system with the control law can be composed of the dynamic surface vector, the filtered error vector, and the upper bound estimation error of the uncertainty as(79)ṡ0=g0s1∗+b1+w0−vrk0s0−1rε0sats0,d0,ṡ1=g1s2+b2+w1−k1s1−e^1s1−g0Ts010g0Ts02T,ṡ2=w2−k2s2−e^2s2−g1Ts1,ḃi=−τi−1bi−ẋiv,i=1,2,e~̇i=−υisiTsi−μiei+μie~i,i=1,2.The Lyapunov function is defined as(80)V=∑i=02Vsi+∑i=12Vbi+∑i=12Vei,where Vsi=1/2siTsi,i=0,1,2, Vbi=1/2biTbi,i=1,2, and Vei=1/2υie~i2,i=1,2.The derivatives of the Lyapunov function (80) with respect to time are
(81)V̇s0=s0Tṡ0=s0Tg0s1∗+b1+s0Tw0−vrs0Tk0s0−1rs0Tε0sats0,d0,(82)V̇s1=s1Tṡ1=s1Tg1s2+b2+s1Tw1−s1Tk1s1−e^1s1Ts1−s1∗Tg0Ts0,(83)V̇s2=s2Tṡ2=s2Tw2−s2Tk2s2−e^2s2Ts2−s2Tg1Ts1,(84)V̇bi=biTḃi,i=1,2,(85)V̇ei=1υie~ie~̇i,i=1,2.According to Young’s inequality, we can get [23]
(86)s0Tg0b1≤s0Tg0b1≤s0g0b1≤s02g022+12b12.Sinceki,i=0,1,2 is a diagonal matrix, the following inequality is established.
(87)−siTkisi≤−kiminsiTsi=−kiminsi2,i=0,1,2,where kimin is the minimum value of the diagonal elements in each gain matrix ki. Since the uncertainty w0 satisfies the following inequality,
(88)w0≤e0.e0 is a positive constant. So, the following inequality is established [7]:
(89)s0Tw0≤e0s01+s02,and
(90)s0Tε0sats0,d0=ε01s01+ε02s02s01>d01ands02>d02,ε01s012d01+ε02s02s01≤d01ands02>d02,ε01s012d01+ε02s022d02s01≤d01ands02≤d02,ε01s01+ε02s022d02s01>d01ands02≤d02.The follow inequality existed as(91)ε01s012d01≤ε01s01,ε02s022d02≤ε02s02,if
(92)ε01d01≥e0maxrs01,ε02d02≥e0maxrs02.Then,(93)ε01s012d01≥e0s012maxrs01≥e0s01r,ε02s022d02≥e0s022maxrs02≥e0s02r.Then(94)re0s01+e0s02≤ε01s012d01+ε02s022d02.Substituting Eq. (91) into Eq. (94) yields
(95)re0s01+e0s02≤ε01s012d01+ε02s022d02.Substituting Eq. (91) into Eq. (95) yields
(96)re0s01+e0s02≤ε01s01+ε02s02,ε01s012d01+ε02s02,ε01s012d01+ε02s022d02,ε01s01+ε02s022d02.Combining Eq. (90) with Eq. (96) yields
(97)re0s01+e0s02≤s0Tε0sats0,d0.Combining Eq. (89) and Eq. (97) yields
(98)rs0Tw0≤s0Tε0sats0,d0.Combining (88) and (90), the following inequalities are true:
(99)s0Tw0−1rs0Tε0sats0,d0≤0.Substituting Eq. (86), Eq. (87), and (99) into Eq. (81) yields
(100)V̇s0≤s0Tg0s1∗+b1−vrs0Tk0s0≤s0Tg0s1∗+s02g022+12b12−vrk0mins02.Then, using Young’s inequality again, we can get [7]
(101)s1Tg1b2≤s1Tg1b2≤s12g122+12b22,and
(102)s1Tw1≤s1Tw1≤s1w1≤e1s1≤e14+e1s12.Substituting Eq. (87), Eq. (101), and Eq. (102) into Eq. (82) yields
(103)V̇s1≤s1Tg1s2+s12g122+12b22+e14+e1s12−k1mins12−e^1s12−s1∗Tg0Ts0≤s1Tg1s2+s12g122+12b22+e~1s12−k1mins12+e14−s1∗Tg0Ts0.Considering that(104)s2Tw2≤e24+e2s22.Substituting Eq. (87) and Eq. (104) into Eq. (83) yields
(105)V̇s2≤−s2Tg1Ts1+e24+e~2s22−k2mins22.Substituting Eq. (67) into Eq. (84) yields
(106)V̇bi=−biTτi−1bi−biTẋiv≤−1τimaxbi2+biẋiv≤−1τimaxbi2+biMi≤Mi2−1τimaxbi2+14.Substituting Eq. (71) into Eq. (85) yields
(107)V̇ei=−e~isiTsi−μiei+μie~i≤−e~isi2−μie~i2+μi2ei2+e~i2=−e~isi2−μie~i22+μi2ei2.The derivatives of the Lyapunov function like Eq. (100), (103), (105), (106), and (107) are summed to obtain [5]
(108)V̇=∑i=02V̇si+∑i=12V̇bi+∑i=12V̇ei≤s0Tg0s1∗+s02g022+12b12−vrk0mins02+s1Tg1s2+s12g122+12b22+e~1s12−k1mins12+e14−s1∗Tg0Ts0−s2Tg1Ts1+e24+e~2s22−k2mins22+∑i=12Mi2−1τimaxbi2+12−∑i=12e~isi2−∑i=12μie~i22+∑i=12μi2ei2=g022−vrk0mins02+g122−k1mins12−k2mins22+∑i=1212+Mi2−1τimaxbi2−∑i=12μie~i22+e14+e24+12+∑i=12μi2ei2.Let(109)c=e14+e24+12+∑i=12μi2ei2.Then,c is a positive constant. Because g0 and g1 are a tight set, that is,
(110)g0≤G0,g1≤G1.Therefore, when parameter selection in the control law Eq. (65) and Eq. (92) satisfies the following inequalities:
(111)k0min≥maxrvγ+G022,k1min≥γ+G122,k2min≥γ2,τimax≤21+γ+2Mi2,i=1,2,μi≥γυi,i=1,2,where γ is a positive constant. The derivatives of the Lyapunov function (80) with respect to time satisfies
(112)V̇≤−γV+c.Then, the function satisfies the following inequality(113)Vt≤V0e−γt+cγ1−e−γt,when
(114)V≥cγ.V̇ is negative. We select state variables as
(115)x=s0Ts1Ts2Tb1Tb2Te~1υ1e~2υ2T.Then,(116)x=2V.That is, whenx≥2c/γ, V̇ was negative. For any given constant υi,i=1,2, as long as the selected element of ki,i=0,1,2 and μi,i=1,2 is large enough and the element of τi,i=1,2 is small enough, it can make γ large enough the same as 2c/γ which is small enough, so that the final bounds of s0, s1, and s2 are sufficiently small. So the boundaries of xF, λ̇D, λ̇T, and β can be small enough. It indicates that the missile can eventually hit the target point at the required terminal angular and the side-slip angle remains near zero during the flight.


## 3.1. IGC Model in the Dive Phase

In this paper, the IGC design of the dive phase of the hypersonic variable span missile is studied. The fixed point on the ground is selected as the target, and the IGC model in the dive phase with the terminal angular constraint is established.


### 3.1.1. Relative Motion of Variable Span Missile and Target

As shown in Figure3, the East–North–Up (ENU) coordinate frame is formed by a plane tangential to the Earth’s surface and is attached to the target. The East axis is labeled X, the North axis is labeled Y, and the up axis is labeled Z. Define the line-of-sight (LOS) coordinate frame so that its origin is on the target where the Sx-axis points to the vehicle, the Sy-axis is in the horizontal plane of the target, and the Sz-axis is determined using the right-hand rule. ηD is the angle between the direction of velocity and the LOS, and γD is the longitudinal velocity azimuth in the dive plane. rT is the distance between the center of mass of the variable span missile and the target point. λD and λT are the angles of elevation and azimuth of LOS.Figure 3
ENU coordinate frame and LOS coordinate frame.The relative dynamic model of the variable span missile and the target point by the line-of-sight angles are expressed as [5]
(19)r¨T=rTλ̇D2+rTλ̇T2cos2λD+arT,λ¨D=−2ṙTλ̇DrT−λ̇T2sinλDcosλD+aλDrT,λ¨T=−2ṙTλ̇TrT+2λ̇Dλ̇TtanλD−aλTrTcosλD.The relationship between the components in the LOS coordinate frame and the trajectory coordinate frame of the missile acceleration vector can be expressed as [6]
(20)aλD=SH2,1aV+SH2,2aθ+SH2,3aσ,aλT=SH3,1aV+SH3,2aθ+SH3,3aσ,where SHi,ji,j=1,2,3are the elements in the conversion matrix of the trajectory coordinate frame to the LOS coordinate frame. i represents the row, and j represents the column.To study the effect of the span variant on the IGC method, the hypersonic missile adopts the BTT control strategy in the dive phase to show the influence of the span variant more clearly. Therefore, the lateral force can be neglected. The influence of the lateral force in Equation (2) is regarded as the uncertainty term. Substituting Eq. (18) into Eq. (2) yields [6]
(21)aV=gHx−D+Fsx′m,aθ=gHy−1mFsy′+1mqS0CLcosγV+wθ,aσ=gHz−1mFsz′+1mqS0CLsinγV+wσ,where the uncertain items wθ and wσ can be expressed as [6]
(22)wθ=cosγVmwL−sinγVmN,wσ=sinγVmwL+cosγVmN.For missiles and fixed targets on the ground, the relative motion between them can be expressed as [6]
(23)ṙT=−vcosηD,rTλ̇D=vsinηD,then
(24)rTλ̇D=vsinηD=vsinλD+γDF=0,where γDF is the terminal angle. Because of λD+γDF<π/2 at the terminal time, it can be obtained as
(25)λD+γDF=0.The terminal angle can be constrained by Eq. (25). By substituting Eq. (20) and Eq. (21) into Eq. (19), the relative dynamics of the target and the hypersonic variable span missile are denoted as
(26)ẋF=x01,ẋ0=f0+g0x1∗+w0.In Eq. (26),
(27)xF=λD+γDF,x0=λ̇Dλ̇TT=x01x02T,x1∗=CLcosγVCLsinγVT,w0=wλ¨Dwλ¨TT,f0=−2ṙTλ̇DrT−λ̇T2sinλDcosλD−2ṙTλ̇TrT+2λ̇Dλ̇TtanλD+1rTSH2,1SH2,2SH2,3−SH3,1cosλD−SH3,2cosλD−SH3,3cosλDaVgHy−1mFsy′gHz−1mFsz′,g0=qS0mrTSH2,2SH2,3−SH3,2cosλD−SH3,3cosλD,w0=1rTSH2,2SH2,3−SH3,2cosλD−SH3,3cosλDwθwσ.


### 3.1.2. Rotational Motion Model

The rotational motion model by the attack, sideslip, and bank angle is [7]
(28)α̇β̇γ̇V=−cosαtanβsinαtanβ1sinαcosα0cosαsecβ−sinαsecβ0ωxωyωz−secβcosγVsinγV−tanβcosγVaθV+−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVaθVcosθ.Then, the rotational motion model that is available for control can be denoted as(29)ẋ1=f1+g1x2+w1.In Eq. (29),
(30)x1=αβγVT,x2=ωxωyωzT,w1=wα̇wβ̇wγ̇VT,(31)f1=−secβcosγVsinγV−tanβcosγVaθV+−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVaσVcosθ,g1=−cosαtanβsinαtanβ1sinαcosα0cosαsecβ−sinαsecβ0,w1=−secβcosγVsinγV−tanβcosγVwθ̇−−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVwσ̇.In Eq. (31), the uncertainty can be expressed as [5]
(32)wθ̇=wθv,wσ̇=−wσvcosθ.


### 3.1.3. Processing of the Rotational Dynamic Model

When ignoring the Earth’s rotation, the launch coordinate frame is equivalent to the ground inertial system. By substituting Eq. (18) into Eq. (11), the simplified state equations with the input of fin deflections can be denoted as [7]
(33)ω̇x=qS0bmxδxIxδx+qS0bmxββ+MSx−Iz−IyIxωyωz+wω̇x,ω̇y=qS0bmyδxIyδy+qS0bmyββ+MSy−Ix−IzIyωxωz+wω̇y,ω̇z=qS0cmzδzIzδz+qS0cmzαα+MSz−Iy−IxIzωxωy+wω̇z.Then, the rotational dynamic model that is available for control can be denoted as(34)ẋ2=f2x1,x2+g2tu+w¯2.In Eq. (34) [5],
(35)u=δaδeδrT,w2=wω̇xwω̇ywω̇zT,(36)f2=qS0bmxββ+MSx−Iz−IyIxωyωzqS0bmyββ+MSy−Ix−IzIyωxωzqS0cmzαα+MSz−Iy−IxIzωxωy,g2=qS0⋅diagbmxδxIx,bmyδyIy,cmzδzIz,w2=−İxωxIx−İyωyIywMz−İzωzIz.The variant rate of the three-axis moment of inertia caused by the variable span is difficult to calculate accurately, so it is regarded as an uncertainty in the Eq. (36).By the combination of Eq. (26), (29), and (34), an IGC model for hypersonic variable span missiles with terminal angular constraints can be obtained as
(37)ẋF=x01,ẋ0=f0+g0x1∗+w0,ẋ1=f1+g1x2+w1,ẋ2=f2+g2u+w2,y=xFx0TT.In Eq. (37), w0, w1, and w2 are bounded function vectors for uncertain positions. Let us assume that there is a set of unknown constants ei that satisfies [7]
(38)wi≤ei,i=0,1,2.The method of solvedx1 by x1∗ will be introduced in Section 2.2.


## 3.1.1. Relative Motion of Variable Span Missile and Target

As shown in Figure3, the East–North–Up (ENU) coordinate frame is formed by a plane tangential to the Earth’s surface and is attached to the target. The East axis is labeled X, the North axis is labeled Y, and the up axis is labeled Z. Define the line-of-sight (LOS) coordinate frame so that its origin is on the target where the Sx-axis points to the vehicle, the Sy-axis is in the horizontal plane of the target, and the Sz-axis is determined using the right-hand rule. ηD is the angle between the direction of velocity and the LOS, and γD is the longitudinal velocity azimuth in the dive plane. rT is the distance between the center of mass of the variable span missile and the target point. λD and λT are the angles of elevation and azimuth of LOS.Figure 3
ENU coordinate frame and LOS coordinate frame.The relative dynamic model of the variable span missile and the target point by the line-of-sight angles are expressed as [5]
(19)r¨T=rTλ̇D2+rTλ̇T2cos2λD+arT,λ¨D=−2ṙTλ̇DrT−λ̇T2sinλDcosλD+aλDrT,λ¨T=−2ṙTλ̇TrT+2λ̇Dλ̇TtanλD−aλTrTcosλD.The relationship between the components in the LOS coordinate frame and the trajectory coordinate frame of the missile acceleration vector can be expressed as [6]
(20)aλD=SH2,1aV+SH2,2aθ+SH2,3aσ,aλT=SH3,1aV+SH3,2aθ+SH3,3aσ,where SHi,ji,j=1,2,3are the elements in the conversion matrix of the trajectory coordinate frame to the LOS coordinate frame. i represents the row, and j represents the column.To study the effect of the span variant on the IGC method, the hypersonic missile adopts the BTT control strategy in the dive phase to show the influence of the span variant more clearly. Therefore, the lateral force can be neglected. The influence of the lateral force in Equation (2) is regarded as the uncertainty term. Substituting Eq. (18) into Eq. (2) yields [6]
(21)aV=gHx−D+Fsx′m,aθ=gHy−1mFsy′+1mqS0CLcosγV+wθ,aσ=gHz−1mFsz′+1mqS0CLsinγV+wσ,where the uncertain items wθ and wσ can be expressed as [6]
(22)wθ=cosγVmwL−sinγVmN,wσ=sinγVmwL+cosγVmN.For missiles and fixed targets on the ground, the relative motion between them can be expressed as [6]
(23)ṙT=−vcosηD,rTλ̇D=vsinηD,then
(24)rTλ̇D=vsinηD=vsinλD+γDF=0,where γDF is the terminal angle. Because of λD+γDF<π/2 at the terminal time, it can be obtained as
(25)λD+γDF=0.The terminal angle can be constrained by Eq. (25). By substituting Eq. (20) and Eq. (21) into Eq. (19), the relative dynamics of the target and the hypersonic variable span missile are denoted as
(26)ẋF=x01,ẋ0=f0+g0x1∗+w0.In Eq. (26),
(27)xF=λD+γDF,x0=λ̇Dλ̇TT=x01x02T,x1∗=CLcosγVCLsinγVT,w0=wλ¨Dwλ¨TT,f0=−2ṙTλ̇DrT−λ̇T2sinλDcosλD−2ṙTλ̇TrT+2λ̇Dλ̇TtanλD+1rTSH2,1SH2,2SH2,3−SH3,1cosλD−SH3,2cosλD−SH3,3cosλDaVgHy−1mFsy′gHz−1mFsz′,g0=qS0mrTSH2,2SH2,3−SH3,2cosλD−SH3,3cosλD,w0=1rTSH2,2SH2,3−SH3,2cosλD−SH3,3cosλDwθwσ.


## 3.1.2. Rotational Motion Model

The rotational motion model by the attack, sideslip, and bank angle is [7]
(28)α̇β̇γ̇V=−cosαtanβsinαtanβ1sinαcosα0cosαsecβ−sinαsecβ0ωxωyωz−secβcosγVsinγV−tanβcosγVaθV+−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVaθVcosθ.Then, the rotational motion model that is available for control can be denoted as(29)ẋ1=f1+g1x2+w1.In Eq. (29),
(30)x1=αβγVT,x2=ωxωyωzT,w1=wα̇wβ̇wγ̇VT,(31)f1=−secβcosγVsinγV−tanβcosγVaθV+−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVaσVcosθ,g1=−cosαtanβsinαtanβ1sinαcosα0cosαsecβ−sinαsecβ0,w1=−secβcosγVsinγV−tanβcosγVwθ̇−−secβcosθsinγVcosθcosγVsinθ+tanβcosθsinγVwσ̇.In Eq. (31), the uncertainty can be expressed as [5]
(32)wθ̇=wθv,wσ̇=−wσvcosθ.


## 3.1.3. Processing of the Rotational Dynamic Model

When ignoring the Earth’s rotation, the launch coordinate frame is equivalent to the ground inertial system. By substituting Eq. (18) into Eq. (11), the simplified state equations with the input of fin deflections can be denoted as [7]
(33)ω̇x=qS0bmxδxIxδx+qS0bmxββ+MSx−Iz−IyIxωyωz+wω̇x,ω̇y=qS0bmyδxIyδy+qS0bmyββ+MSy−Ix−IzIyωxωz+wω̇y,ω̇z=qS0cmzδzIzδz+qS0cmzαα+MSz−Iy−IxIzωxωy+wω̇z.Then, the rotational dynamic model that is available for control can be denoted as(34)ẋ2=f2x1,x2+g2tu+w¯2.In Eq. (34) [5],
(35)u=δaδeδrT,w2=wω̇xwω̇ywω̇zT,(36)f2=qS0bmxββ+MSx−Iz−IyIxωyωzqS0bmyββ+MSy−Ix−IzIyωxωzqS0cmzαα+MSz−Iy−IxIzωxωy,g2=qS0⋅diagbmxδxIx,bmyδyIy,cmzδzIz,w2=−İxωxIx−İyωyIywMz−İzωzIz.The variant rate of the three-axis moment of inertia caused by the variable span is difficult to calculate accurately, so it is regarded as an uncertainty in the Eq. (36).By the combination of Eq. (26), (29), and (34), an IGC model for hypersonic variable span missiles with terminal angular constraints can be obtained as
(37)ẋF=x01,ẋ0=f0+g0x1∗+w0,ẋ1=f1+g1x2+w1,ẋ2=f2+g2u+w2,y=xFx0TT.In Eq. (37), w0, w1, and w2 are bounded function vectors for uncertain positions. Let us assume that there is a set of unknown constants ei that satisfies [7]
(38)wi≤ei,i=0,1,2.The method of solvedx1 by x1∗ will be introduced in Section 2.2.


## 3.2. The IGC Method in the Dive Phase with Variable Span Control

The IGC model in the dive phase guidance control with terminal angular constraint is shown in Eq. (37). The adaptive dynamic surface back-stepping control method is designed to make the system stable and make the output as close as possible to zero and make the system have strong robustness for uncertain factors.


### 3.2.1. Design of the Adaptive Dynamic Surface Back-Stepping IGC Method

To design an adaptive dynamic surface back-stepping control method, Eq. (37) is required to nonsingular [23] as
(39)detg1=αsecβ≠0.Therefore, let us assume to makeg1 nonsingularity as
(40)α,β∈Ω:α,β0<α<π2,β<π2,in the whole dive phase. And Ω is the bounded closed set in R2, so Ω is a tight set.The design steps of the adaptive dynamic surface back-stepping IGC method are as follows [6].Step 1.
Design the dynamic surface [5]
(41)s0=λ̇D+kFvrxFλ̇TT=s01s02T,where kF is the angular error coefficient, of which the magnitude determines the weight of the terminal angular error in the dynamic surface.
The rate of change of the dynamic surfaces0 is designed according to the exponential approach law.
(42)ṡ0=−vk01rs01−ε01rsats01,d01−vk02rs02−ε02rsats02,d02,where k01 and k02 are positive gain constants, ε01 andε02 are the gain-of-saturation function, and d01 and d02 are the boundary layer thickness. sats,d is a saturation function which defined as [7]
(43)sats,d=1,s>d,sd,s≤d,−1,s<−d.
When the missile is far from the target point andr is large, the rate of approaching law in (42) will be slow so that the missile has a small overload in the initial stage. When the missile approaches the target point and r is small, the rate of approaching law in (42) increases so that λ̇D does not diverge to improve the accuracy of hitting. The selection of the saturation function can effectively eliminate the chattering of the first virtual control variable.
Deriving the dynamic surfaces0 and then combining (42) can get the first virtual control input x1v as
(44)x1v∗=g0−1−vrk0s0−1rε0sats0,d0−f0′−f0.
In Eq. (44), the result of deriving the terminal angle is
(45)f0′=kFvrx01+kFv̇r−ṙvr2xF0.
In Eq. (44), k0=diagk01,k02, ε0=diagε01,ε02, and d0=d01d02T.
According tox1v∗=CLvcosγVvCLvsinγVvT, it can be translated as
(46)γVv=tan−1x1v∗2x1v∗1,CLv=x1v∗1cosγVv.
According to the aerodynamic model of the variable span missile, the structural form of the lift coefficient caused by the attack angle is(47)CLv=CL0+CLααv+CLα3αv3.
It can be solved as(48)αv=CLv−CL0−CLα3αv,t−13CLα.
In Eq. (48), take αv as the main factor affecting for lift, and αv,t−13 is the previous moment attack angle.
Then,(49)x1v=αvγVvT.
In order to exert the advantage of the variable span hypersonic missile, the morphing rate and attack angle of the wing are coordinated controls. The roadmap of the IGC method for a hypersonic morphing missile based on variable span auxiliary control is shown in Figure4, while the IGC method with invariable span is shown in Figure 5.
The principle of coordinated control shown in Figure4 is to translate the change of the lift coefficient into a change of the attack angle and the morphing rate, so that the span variant can withstand the change demand of the lift command and reduce the change demand of the attack angle control. While in the traditional IGC method with invariable span, the change demand of the lift command is only reached by the attack angle. It would improve the stability and velocity of missile IGC mission. That is,
(50)ĊLv=ĊLαξαv+CLαξα̇v.
Solve Eq. (50) to get
(51)ĊLαξ=ĊLv−CLαξt−1α̇v−ĊL0,t−1ξt−1−ĊL,t−1α3ξt−1αv3αv.
In Eq. (51), take αv as the main factor affecting for lift; ξt−1,ĊL0,t−1,ĊL,t−1α3 is the previous moment value of each variable.
According to the model of the aerodynamic coefficient of the missile,(52)CLαξ=1ξ⋅CLα,2×3ξ,where CLαξ is the coefficient matrix of CLα with respect to the morphing rate, and the expression is
(53)CLαξ=CLα,1,1ξCLα,1,2ξCLα,1,3ξCLα,2,1ξCLα2,2ξCLα,2,3ξ.
Let us substitute Eq. (53) into (3.35) and deriving on both sides to get
(54)ĊLαξ=ξ̇⋅CLα2,2ξ.
Combining Eq. (51) and Eq. (54), the variation law of the morphing rate can be obtained as
(55)ξ̇c=ĊLv−CLαξt−1α̇v−ĊL0,t−1ξt−1−ĊL,t−1α3ξt−1αv3αv⋅CLα2,2ξ.
In order to prevent differential explosion, it will pass through the first-order filter as(56)τξξ¨+ξ̇=ξ̇cξ̇0=ξ̇c0.Step 2.
Design the dynamic surface(57)s1=x1−x1d∗,where x1d∗=x1d10x1d2T, x1d is the output value of the following first-order filter as
(58)τ1ẋ1d+x1d=x1vx1d0=x1v0,where τ1=diagτ11,0,τ13 is the time constant of the filter.
Design the second virtual control input as(59)x2v=−g1−1f1+e^1s1+k1s1−ẋ1d∗+g0Ts010g0Ts02T,where ẋ1d∗=ẋ1d10ẋ1d2T, k1=diagk11,k12,k13, and e^1 is an estimate of the upper bound e1 of the uncertainty w1, which can be obtained by the following adaptive law
(60)e^̇1=υ1s1Ts1−μ1e^1,e^10=0,where υ1 and μ1 are positive constants.Step 3.
Design the dynamic surface(61)s2=x2−x2d.
Similarly,x2d is the output value of the following first-order filter as
(62)τ2ẋ2d+x2d=x2vx2d0=x2v0,where τ2=diagτ21,τ22,τ23 is the time constant of the filter. The control input of the system is designed as
(63)u=−g2−1f2+e^2s2+k2s2−ẋ2d+g1Ts1,where k2=diagk21,k22,k23 and e^2 is an estimate of the upper bound e2 of the uncertainty w2, which can be obtained by the following adaptive law as
(64)e^̇2=υ2s2Ts2−μ2e^2,e^20=0,where υ2 and μ2 are positive constants.
In summary, the adaptive dynamic surface back-stepping IGC method for variable span control of the system (3.20) can be expressed as(65)s0=x01+kFvRxFx02T,x1v=g0−1−vrk0s0−1rε0sats0,d0−f0′−f0,ξ̇c=ĊLv−CLαξα̇vαvTαvαvT−1CLα,2ξTCLα,2ξCLα,2ξT−1,τξξ¨+ξ̇=ξ̇cξ̇0=ξ̇c0,τ1ẋ1d+x1d=x1v,x1d0=x1v0,x1d∗=x1d10x1d2T,ẋ1d∗=ẋ1d10ẋ1d2T,s1=x1−x1d∗,x2v=−g1−1f1+e^1s1+k1s1−ẋ1d∗+g0Ts010g0Ts02T,e^̇1=υ1s1Ts1−μ1e^1,e^10=0,τ2ẋ2d+x2d=x2v,x2d0=x2v0,s2=x2−x2d,u=−g2−1f2+e^2s2+k2s2−ẋ2d+g1Ts1,e^̇2=υ2s2Ts2−μ2e^2,e^20=0.Figure 4
The roadmap of the IGC method with variable span auxiliary control.Figure 5
The roadmap of the IGC method with invariable span.


### 3.2.2. Stability Analysis

The stability of the closed-loop system obtained using the control law in Eq. (65) is analyzed below.Define the filter errorb1 and b2, respectively, as
(66)bi=xid−xiv,i=1,2.Then, the derivative of the filtering error is(67)ḃi=τi−1xiv−xid−ẋiv=−τi−1bi−ẋiv,i=1,2.AssumingBi=−ẋiv,i=1,2, according to the assumption in Section 3.2.1, Bi is a tight set, then it has an upper bound as [7]
(68)Bi≤Mi,i=1,2.We define the upper bound estimation error of the norm as [7]
(69)e~i=ei−e^i,i=1,2.Then, its derivative is(70)e~̇i=−e^̇i,i=1,2.Substituting the parameter adaptive law (60) and Eq. (64) into Eq. (70) yields
(71)e~̇i=−υisiTsi−μie^i=−υisiTsi−μiei+μie~i,i=1,2.s1∗ can be defined as
(72)s1∗=x1∗−x1d=x1∗−b1+x1v.The derivative ofs0 is
(73)ṡ0=f0+g0s1∗+b1+x1v+w0+f0′.Substituting Eq. (44) into Eq. (73) yields
(74)ṡ0=g0s1∗+b1+w0−vrk0s0−1rε0sats0,d0.The derivative ofs1 is
(75)ṡ1=f1+g1s2+x2v+b2+w1−ẋ1d∗.Substituting Eq. (59) into Eq. (75) yields
(76)ṡ1=g1s2+b2+w1−k1s1−e^1s1−g0Ts010g0Ts02T.The derivative ofs2 is
(77)ṡ2=f2+g2u+w2−ẋ2d.Now the expression of the control inputu can be obtained by substituting Eq. (77) into Eq. (63) as
(78)ṡ2=w2−k2s2−e^2s2−g1Ts1In summary, the entire closed-loop system with the control law can be composed of the dynamic surface vector, the filtered error vector, and the upper bound estimation error of the uncertainty as(79)ṡ0=g0s1∗+b1+w0−vrk0s0−1rε0sats0,d0,ṡ1=g1s2+b2+w1−k1s1−e^1s1−g0Ts010g0Ts02T,ṡ2=w2−k2s2−e^2s2−g1Ts1,ḃi=−τi−1bi−ẋiv,i=1,2,e~̇i=−υisiTsi−μiei+μie~i,i=1,2.The Lyapunov function is defined as(80)V=∑i=02Vsi+∑i=12Vbi+∑i=12Vei,where Vsi=1/2siTsi,i=0,1,2, Vbi=1/2biTbi,i=1,2, and Vei=1/2υie~i2,i=1,2.The derivatives of the Lyapunov function (80) with respect to time are
(81)V̇s0=s0Tṡ0=s0Tg0s1∗+b1+s0Tw0−vrs0Tk0s0−1rs0Tε0sats0,d0,(82)V̇s1=s1Tṡ1=s1Tg1s2+b2+s1Tw1−s1Tk1s1−e^1s1Ts1−s1∗Tg0Ts0,(83)V̇s2=s2Tṡ2=s2Tw2−s2Tk2s2−e^2s2Ts2−s2Tg1Ts1,(84)V̇bi=biTḃi,i=1,2,(85)V̇ei=1υie~ie~̇i,i=1,2.According to Young’s inequality, we can get [23]
(86)s0Tg0b1≤s0Tg0b1≤s0g0b1≤s02g022+12b12.Sinceki,i=0,1,2 is a diagonal matrix, the following inequality is established.
(87)−siTkisi≤−kiminsiTsi=−kiminsi2,i=0,1,2,where kimin is the minimum value of the diagonal elements in each gain matrix ki. Since the uncertainty w0 satisfies the following inequality,
(88)w0≤e0.e0 is a positive constant. So, the following inequality is established [7]:
(89)s0Tw0≤e0s01+s02,and
(90)s0Tε0sats0,d0=ε01s01+ε02s02s01>d01ands02>d02,ε01s012d01+ε02s02s01≤d01ands02>d02,ε01s012d01+ε02s022d02s01≤d01ands02≤d02,ε01s01+ε02s022d02s01>d01ands02≤d02.The follow inequality existed as(91)ε01s012d01≤ε01s01,ε02s022d02≤ε02s02,if
(92)ε01d01≥e0maxrs01,ε02d02≥e0maxrs02.Then,(93)ε01s012d01≥e0s012maxrs01≥e0s01r,ε02s022d02≥e0s022maxrs02≥e0s02r.Then(94)re0s01+e0s02≤ε01s012d01+ε02s022d02.Substituting Eq. (91) into Eq. (94) yields
(95)re0s01+e0s02≤ε01s012d01+ε02s022d02.Substituting Eq. (91) into Eq. (95) yields
(96)re0s01+e0s02≤ε01s01+ε02s02,ε01s012d01+ε02s02,ε01s012d01+ε02s022d02,ε01s01+ε02s022d02.Combining Eq. (90) with Eq. (96) yields
(97)re0s01+e0s02≤s0Tε0sats0,d0.Combining Eq. (89) and Eq. (97) yields
(98)rs0Tw0≤s0Tε0sats0,d0.Combining (88) and (90), the following inequalities are true:
(99)s0Tw0−1rs0Tε0sats0,d0≤0.Substituting Eq. (86), Eq. (87), and (99) into Eq. (81) yields
(100)V̇s0≤s0Tg0s1∗+b1−vrs0Tk0s0≤s0Tg0s1∗+s02g022+12b12−vrk0mins02.Then, using Young’s inequality again, we can get [7]
(101)s1Tg1b2≤s1Tg1b2≤s12g122+12b22,and
(102)s1Tw1≤s1Tw1≤s1w1≤e1s1≤e14+e1s12.Substituting Eq. (87), Eq. (101), and Eq. (102) into Eq. (82) yields
(103)V̇s1≤s1Tg1s2+s12g122+12b22+e14+e1s12−k1mins12−e^1s12−s1∗Tg0Ts0≤s1Tg1s2+s12g122+12b22+e~1s12−k1mins12+e14−s1∗Tg0Ts0.Considering that(104)s2Tw2≤e24+e2s22.Substituting Eq. (87) and Eq. (104) into Eq. (83) yields
(105)V̇s2≤−s2Tg1Ts1+e24+e~2s22−k2mins22.Substituting Eq. (67) into Eq. (84) yields
(106)V̇bi=−biTτi−1bi−biTẋiv≤−1τimaxbi2+biẋiv≤−1τimaxbi2+biMi≤Mi2−1τimaxbi2+14.Substituting Eq. (71) into Eq. (85) yields
(107)V̇ei=−e~isiTsi−μiei+μie~i≤−e~isi2−μie~i2+μi2ei2+e~i2=−e~isi2−μie~i22+μi2ei2.The derivatives of the Lyapunov function like Eq. (100), (103), (105), (106), and (107) are summed to obtain [5]
(108)V̇=∑i=02V̇si+∑i=12V̇bi+∑i=12V̇ei≤s0Tg0s1∗+s02g022+12b12−vrk0mins02+s1Tg1s2+s12g122+12b22+e~1s12−k1mins12+e14−s1∗Tg0Ts0−s2Tg1Ts1+e24+e~2s22−k2mins22+∑i=12Mi2−1τimaxbi2+12−∑i=12e~isi2−∑i=12μie~i22+∑i=12μi2ei2=g022−vrk0mins02+g122−k1mins12−k2mins22+∑i=1212+Mi2−1τimaxbi2−∑i=12μie~i22+e14+e24+12+∑i=12μi2ei2.Let(109)c=e14+e24+12+∑i=12μi2ei2.Then,c is a positive constant. Because g0 and g1 are a tight set, that is,
(110)g0≤G0,g1≤G1.Therefore, when parameter selection in the control law Eq. (65) and Eq. (92) satisfies the following inequalities:
(111)k0min≥maxrvγ+G022,k1min≥γ+G122,k2min≥γ2,τimax≤21+γ+2Mi2,i=1,2,μi≥γυi,i=1,2,where γ is a positive constant. The derivatives of the Lyapunov function (80) with respect to time satisfies
(112)V̇≤−γV+c.Then, the function satisfies the following inequality(113)Vt≤V0e−γt+cγ1−e−γt,when
(114)V≥cγ.V̇ is negative. We select state variables as
(115)x=s0Ts1Ts2Tb1Tb2Te~1υ1e~2υ2T.Then,(116)x=2V.That is, whenx≥2c/γ, V̇ was negative. For any given constant υi,i=1,2, as long as the selected element of ki,i=0,1,2 and μi,i=1,2 is large enough and the element of τi,i=1,2 is small enough, it can make γ large enough the same as 2c/γ which is small enough, so that the final bounds of s0, s1, and s2 are sufficiently small. So the boundaries of xF, λ̇D, λ̇T, and β can be small enough. It indicates that the missile can eventually hit the target point at the required terminal angular and the side-slip angle remains near zero during the flight.


## 3.2.1. Design of the Adaptive Dynamic Surface Back-Stepping IGC Method

To design an adaptive dynamic surface back-stepping control method, Eq. (37) is required to nonsingular [23] as
(39)detg1=αsecβ≠0.Therefore, let us assume to makeg1 nonsingularity as
(40)α,β∈Ω:α,β0<α<π2,β<π2,in the whole dive phase. And Ω is the bounded closed set in R2, so Ω is a tight set.The design steps of the adaptive dynamic surface back-stepping IGC method are as follows [6].Step 1.
Design the dynamic surface [5]
(41)s0=λ̇D+kFvrxFλ̇TT=s01s02T,where kF is the angular error coefficient, of which the magnitude determines the weight of the terminal angular error in the dynamic surface.
The rate of change of the dynamic surfaces0 is designed according to the exponential approach law.
(42)ṡ0=−vk01rs01−ε01rsats01,d01−vk02rs02−ε02rsats02,d02,where k01 and k02 are positive gain constants, ε01 andε02 are the gain-of-saturation function, and d01 and d02 are the boundary layer thickness. sats,d is a saturation function which defined as [7]
(43)sats,d=1,s>d,sd,s≤d,−1,s<−d.
When the missile is far from the target point andr is large, the rate of approaching law in (42) will be slow so that the missile has a small overload in the initial stage. When the missile approaches the target point and r is small, the rate of approaching law in (42) increases so that λ̇D does not diverge to improve the accuracy of hitting. The selection of the saturation function can effectively eliminate the chattering of the first virtual control variable.
Deriving the dynamic surfaces0 and then combining (42) can get the first virtual control input x1v as
(44)x1v∗=g0−1−vrk0s0−1rε0sats0,d0−f0′−f0.
In Eq. (44), the result of deriving the terminal angle is
(45)f0′=kFvrx01+kFv̇r−ṙvr2xF0.
In Eq. (44), k0=diagk01,k02, ε0=diagε01,ε02, and d0=d01d02T.
According tox1v∗=CLvcosγVvCLvsinγVvT, it can be translated as
(46)γVv=tan−1x1v∗2x1v∗1,CLv=x1v∗1cosγVv.
According to the aerodynamic model of the variable span missile, the structural form of the lift coefficient caused by the attack angle is(47)CLv=CL0+CLααv+CLα3αv3.
It can be solved as(48)αv=CLv−CL0−CLα3αv,t−13CLα.
In Eq. (48), take αv as the main factor affecting for lift, and αv,t−13 is the previous moment attack angle.
Then,(49)x1v=αvγVvT.
In order to exert the advantage of the variable span hypersonic missile, the morphing rate and attack angle of the wing are coordinated controls. The roadmap of the IGC method for a hypersonic morphing missile based on variable span auxiliary control is shown in Figure4, while the IGC method with invariable span is shown in Figure 5.
The principle of coordinated control shown in Figure4 is to translate the change of the lift coefficient into a change of the attack angle and the morphing rate, so that the span variant can withstand the change demand of the lift command and reduce the change demand of the attack angle control. While in the traditional IGC method with invariable span, the change demand of the lift command is only reached by the attack angle. It would improve the stability and velocity of missile IGC mission. That is,
(50)ĊLv=ĊLαξαv+CLαξα̇v.
Solve Eq. (50) to get
(51)ĊLαξ=ĊLv−CLαξt−1α̇v−ĊL0,t−1ξt−1−ĊL,t−1α3ξt−1αv3αv.
In Eq. (51), take αv as the main factor affecting for lift; ξt−1,ĊL0,t−1,ĊL,t−1α3 is the previous moment value of each variable.
According to the model of the aerodynamic coefficient of the missile,(52)CLαξ=1ξ⋅CLα,2×3ξ,where CLαξ is the coefficient matrix of CLα with respect to the morphing rate, and the expression is
(53)CLαξ=CLα,1,1ξCLα,1,2ξCLα,1,3ξCLα,2,1ξCLα2,2ξCLα,2,3ξ.
Let us substitute Eq. (53) into (3.35) and deriving on both sides to get
(54)ĊLαξ=ξ̇⋅CLα2,2ξ.
Combining Eq. (51) and Eq. (54), the variation law of the morphing rate can be obtained as
(55)ξ̇c=ĊLv−CLαξt−1α̇v−ĊL0,t−1ξt−1−ĊL,t−1α3ξt−1αv3αv⋅CLα2,2ξ.
In order to prevent differential explosion, it will pass through the first-order filter as(56)τξξ¨+ξ̇=ξ̇cξ̇0=ξ̇c0.Step 2.
Design the dynamic surface(57)s1=x1−x1d∗,where x1d∗=x1d10x1d2T, x1d is the output value of the following first-order filter as
(58)τ1ẋ1d+x1d=x1vx1d0=x1v0,where τ1=diagτ11,0,τ13 is the time constant of the filter.
Design the second virtual control input as(59)x2v=−g1−1f1+e^1s1+k1s1−ẋ1d∗+g0Ts010g0Ts02T,where ẋ1d∗=ẋ1d10ẋ1d2T, k1=diagk11,k12,k13, and e^1 is an estimate of the upper bound e1 of the uncertainty w1, which can be obtained by the following adaptive law
(60)e^̇1=υ1s1Ts1−μ1e^1,e^10=0,where υ1 and μ1 are positive constants.Step 3.
Design the dynamic surface(61)s2=x2−x2d.
Similarly,x2d is the output value of the following first-order filter as
(62)τ2ẋ2d+x2d=x2vx2d0=x2v0,where τ2=diagτ21,τ22,τ23 is the time constant of the filter. The control input of the system is designed as
(63)u=−g2−1f2+e^2s2+k2s2−ẋ2d+g1Ts1,where k2=diagk21,k22,k23 and e^2 is an estimate of the upper bound e2 of the uncertainty w2, which can be obtained by the following adaptive law as
(64)e^̇2=υ2s2Ts2−μ2e^2,e^20=0,where υ2 and μ2 are positive constants.
In summary, the adaptive dynamic surface back-stepping IGC method for variable span control of the system (3.20) can be expressed as(65)s0=x01+kFvRxFx02T,x1v=g0−1−vrk0s0−1rε0sats0,d0−f0′−f0,ξ̇c=ĊLv−CLαξα̇vαvTαvαvT−1CLα,2ξTCLα,2ξCLα,2ξT−1,τξξ¨+ξ̇=ξ̇cξ̇0=ξ̇c0,τ1ẋ1d+x1d=x1v,x1d0=x1v0,x1d∗=x1d10x1d2T,ẋ1d∗=ẋ1d10ẋ1d2T,s1=x1−x1d∗,x2v=−g1−1f1+e^1s1+k1s1−ẋ1d∗+g0Ts010g0Ts02T,e^̇1=υ1s1Ts1−μ1e^1,e^10=0,τ2ẋ2d+x2d=x2v,x2d0=x2v0,s2=x2−x2d,u=−g2−1f2+e^2s2+k2s2−ẋ2d+g1Ts1,e^̇2=υ2s2Ts2−μ2e^2,e^20=0.Figure 4
The roadmap of the IGC method with variable span auxiliary control.Figure 5
The roadmap of the IGC method with invariable span.


## 3.2.2. Stability Analysis

The stability of the closed-loop system obtained using the control law in Eq. (65) is analyzed below.Define the filter errorb1 and b2, respectively, as
(66)bi=xid−xiv,i=1,2.Then, the derivative of the filtering error is(67)ḃi=τi−1xiv−xid−ẋiv=−τi−1bi−ẋiv,i=1,2.AssumingBi=−ẋiv,i=1,2, according to the assumption in Section 3.2.1, Bi is a tight set, then it has an upper bound as [7]
(68)Bi≤Mi,i=1,2.We define the upper bound estimation error of the norm as [7]
(69)e~i=ei−e^i,i=1,2.Then, its derivative is(70)e~̇i=−e^̇i,i=1,2.Substituting the parameter adaptive law (60) and Eq. (64) into Eq. (70) yields
(71)e~̇i=−υisiTsi−μie^i=−υisiTsi−μiei+μie~i,i=1,2.s1∗ can be defined as
(72)s1∗=x1∗−x1d=x1∗−b1+x1v.The derivative ofs0 is
(73)ṡ0=f0+g0s1∗+b1+x1v+w0+f0′.Substituting Eq. (44) into Eq. (73) yields
(74)ṡ0=g0s1∗+b1+w0−vrk0s0−1rε0sats0,d0.The derivative ofs1 is
(75)ṡ1=f1+g1s2+x2v+b2+w1−ẋ1d∗.Substituting Eq. (59) into Eq. (75) yields
(76)ṡ1=g1s2+b2+w1−k1s1−e^1s1−g0Ts010g0Ts02T.The derivative ofs2 is
(77)ṡ2=f2+g2u+w2−ẋ2d.Now the expression of the control inputu can be obtained by substituting Eq. (77) into Eq. (63) as
(78)ṡ2=w2−k2s2−e^2s2−g1Ts1In summary, the entire closed-loop system with the control law can be composed of the dynamic surface vector, the filtered error vector, and the upper bound estimation error of the uncertainty as(79)ṡ0=g0s1∗+b1+w0−vrk0s0−1rε0sats0,d0,ṡ1=g1s2+b2+w1−k1s1−e^1s1−g0Ts010g0Ts02T,ṡ2=w2−k2s2−e^2s2−g1Ts1,ḃi=−τi−1bi−ẋiv,i=1,2,e~̇i=−υisiTsi−μiei+μie~i,i=1,2.The Lyapunov function is defined as(80)V=∑i=02Vsi+∑i=12Vbi+∑i=12Vei,where Vsi=1/2siTsi,i=0,1,2, Vbi=1/2biTbi,i=1,2, and Vei=1/2υie~i2,i=1,2.The derivatives of the Lyapunov function (80) with respect to time are
(81)V̇s0=s0Tṡ0=s0Tg0s1∗+b1+s0Tw0−vrs0Tk0s0−1rs0Tε0sats0,d0,(82)V̇s1=s1Tṡ1=s1Tg1s2+b2+s1Tw1−s1Tk1s1−e^1s1Ts1−s1∗Tg0Ts0,(83)V̇s2=s2Tṡ2=s2Tw2−s2Tk2s2−e^2s2Ts2−s2Tg1Ts1,(84)V̇bi=biTḃi,i=1,2,(85)V̇ei=1υie~ie~̇i,i=1,2.According to Young’s inequality, we can get [23]
(86)s0Tg0b1≤s0Tg0b1≤s0g0b1≤s02g022+12b12.Sinceki,i=0,1,2 is a diagonal matrix, the following inequality is established.
(87)−siTkisi≤−kiminsiTsi=−kiminsi2,i=0,1,2,where kimin is the minimum value of the diagonal elements in each gain matrix ki. Since the uncertainty w0 satisfies the following inequality,
(88)w0≤e0.e0 is a positive constant. So, the following inequality is established [7]:
(89)s0Tw0≤e0s01+s02,and
(90)s0Tε0sats0,d0=ε01s01+ε02s02s01>d01ands02>d02,ε01s012d01+ε02s02s01≤d01ands02>d02,ε01s012d01+ε02s022d02s01≤d01ands02≤d02,ε01s01+ε02s022d02s01>d01ands02≤d02.The follow inequality existed as(91)ε01s012d01≤ε01s01,ε02s022d02≤ε02s02,if
(92)ε01d01≥e0maxrs01,ε02d02≥e0maxrs02.Then,(93)ε01s012d01≥e0s012maxrs01≥e0s01r,ε02s022d02≥e0s022maxrs02≥e0s02r.Then(94)re0s01+e0s02≤ε01s012d01+ε02s022d02.Substituting Eq. (91) into Eq. (94) yields
(95)re0s01+e0s02≤ε01s012d01+ε02s022d02.Substituting Eq. (91) into Eq. (95) yields
(96)re0s01+e0s02≤ε01s01+ε02s02,ε01s012d01+ε02s02,ε01s012d01+ε02s022d02,ε01s01+ε02s022d02.Combining Eq. (90) with Eq. (96) yields
(97)re0s01+e0s02≤s0Tε0sats0,d0.Combining Eq. (89) and Eq. (97) yields
(98)rs0Tw0≤s0Tε0sats0,d0.Combining (88) and (90), the following inequalities are true:
(99)s0Tw0−1rs0Tε0sats0,d0≤0.Substituting Eq. (86), Eq. (87), and (99) into Eq. (81) yields
(100)V̇s0≤s0Tg0s1∗+b1−vrs0Tk0s0≤s0Tg0s1∗+s02g022+12b12−vrk0mins02.Then, using Young’s inequality again, we can get [7]
(101)s1Tg1b2≤s1Tg1b2≤s12g122+12b22,and
(102)s1Tw1≤s1Tw1≤s1w1≤e1s1≤e14+e1s12.Substituting Eq. (87), Eq. (101), and Eq. (102) into Eq. (82) yields
(103)V̇s1≤s1Tg1s2+s12g122+12b22+e14+e1s12−k1mins12−e^1s12−s1∗Tg0Ts0≤s1Tg1s2+s12g122+12b22+e~1s12−k1mins12+e14−s1∗Tg0Ts0.Considering that(104)s2Tw2≤e24+e2s22.Substituting Eq. (87) and Eq. (104) into Eq. (83) yields
(105)V̇s2≤−s2Tg1Ts1+e24+e~2s22−k2mins22.Substituting Eq. (67) into Eq. (84) yields
(106)V̇bi=−biTτi−1bi−biTẋiv≤−1τimaxbi2+biẋiv≤−1τimaxbi2+biMi≤Mi2−1τimaxbi2+14.Substituting Eq. (71) into Eq. (85) yields
(107)V̇ei=−e~isiTsi−μiei+μie~i≤−e~isi2−μie~i2+μi2ei2+e~i2=−e~isi2−μie~i22+μi2ei2.The derivatives of the Lyapunov function like Eq. (100), (103), (105), (106), and (107) are summed to obtain [5]
(108)V̇=∑i=02V̇si+∑i=12V̇bi+∑i=12V̇ei≤s0Tg0s1∗+s02g022+12b12−vrk0mins02+s1Tg1s2+s12g122+12b22+e~1s12−k1mins12+e14−s1∗Tg0Ts0−s2Tg1Ts1+e24+e~2s22−k2mins22+∑i=12Mi2−1τimaxbi2+12−∑i=12e~isi2−∑i=12μie~i22+∑i=12μi2ei2=g022−vrk0mins02+g122−k1mins12−k2mins22+∑i=1212+Mi2−1τimaxbi2−∑i=12μie~i22+e14+e24+12+∑i=12μi2ei2.Let(109)c=e14+e24+12+∑i=12μi2ei2.Then,c is a positive constant. Because g0 and g1 are a tight set, that is,
(110)g0≤G0,g1≤G1.Therefore, when parameter selection in the control law Eq. (65) and Eq. (92) satisfies the following inequalities:
(111)k0min≥maxrvγ+G022,k1min≥γ+G122,k2min≥γ2,τimax≤21+γ+2Mi2,i=1,2,μi≥γυi,i=1,2,where γ is a positive constant. The derivatives of the Lyapunov function (80) with respect to time satisfies
(112)V̇≤−γV+c.Then, the function satisfies the following inequality(113)Vt≤V0e−γt+cγ1−e−γt,when
(114)V≥cγ.V̇ is negative. We select state variables as
(115)x=s0Ts1Ts2Tb1Tb2Te~1υ1e~2υ2T.Then,(116)x=2V.That is, whenx≥2c/γ, V̇ was negative. For any given constant υi,i=1,2, as long as the selected element of ki,i=0,1,2 and μi,i=1,2 is large enough and the element of τi,i=1,2 is small enough, it can make γ large enough the same as 2c/γ which is small enough, so that the final bounds of s0, s1, and s2 are sufficiently small. So the boundaries of xF, λ̇D, λ̇T, and β can be small enough. It indicates that the missile can eventually hit the target point at the required terminal angular and the side-slip angle remains near zero during the flight.


## 4. Simulation Analysis

### 4.1. Performance Analysis of IGC Method for a Variable Span Missile

In order to verify the effectiveness of the IGC method with a variable span auxiliary control in the dive phase of the hypersonic missile, the model is numerically simulated and the control group is set up to observe the effect of the IGC method under the condition of variable span and invariable span. The initial state is shown in Table2.Table 2
The initial state of the missile and the latitude and longitude of the target point.
Name
Value
Name
Value
Name
Value
v0 (m/s)
3000
α0 (°)
2
ϕ0 (°)
0
θ0 (°)
0
β0 (°)
0
ψ0 (°)
0
σ0 (°)
0
γV0 (°)
0
ϕT (°)
1
ωx0 (o/s)
0
x (km)
0
ψT (°)
1
ωy0 (°/s)
0
y (km)
30
ξ10
0.40
ωz0 (∘/s)
0
z (km)
0
ξ20
0.55
λF (°)
60During the simulation process, it is necessary to control the amplitude and rate of change of the attack angle command and the rate of change of the bank angle command as follows.(117)0<αv≤15∘/s,α̇v≤5°/s,γ̇Vv≤20°/s.At the same time, the fin deflections of the missile are restricted as follows:(118)δx≤5∘,δy≤20∘,δz≤20∘.The rate of change of the morphing rate of the wing is as follows.(119)ξ̇≤30%30%/s.Firstly, the control parameters in Eq. (65) are optimized based on the genetic algorithm. The miss distance of the missile, the deviation of the terminal angular, the jitter of the attitude rates, and the jitter of fin deflections are optimized as fitness functions. The attitude rates and the jitter of the fin deflections can be expressed by the sum of the absolute values of their respective changes as
(120)fit=kddmiss+kθxF+kω∫0tdωx+dωy+dωz+kδ∫0tdδx+dδy+dδz,where dmiss is the missile miss distance and t is the flight time of the missile in the dive phase. kd, kθ, kω, and kδare the weight coefficients for measuring the miss distance, the error of the terminal angular, the jitter of rate of the angle, and the jitter of fin deflections. The initial population size is set to 60. After genetic optimization, the optimal controller parameters were obtained. The average fitness value changes with the genetic algebra as shown in Figures 6 and 7. The optimized parameter values are shown in Tables 2 and 3.Figure 6
The fitness value with genetic algebra of the variable span missile.Figure 7
The fitness value with the genetic algebra of the invariable span missile.Table 3
Parameter optimization value of the IGC method for variable span missile.
Name
Value
Name
Value
Name
Value
k01
2.56
k11
10.12
k21
2.25
k02
0.51
k12
3.12
k22
2.78
ε01
0.62
k13
2.41
k23
5.32
ε02
0.33
τ11
0.71
τ21
0.44
d01
0.44
τ12
0.50
τ22
0.22
d02
0.27
υ1
50.11
τ23
1.3
kF
0.33
μ1
2.36
υ2
47.29
τξ
0.21
μ2
1.2The simulation results are as follows.Compare the simulation results obtained by the variable span missile with the invariable span missile as follows.In Table5, γD is the final terminal angular. Vt is the final velocity. T is the total flight time. Qt is the final dynamic pressure. δ¯i,i=x,y,z are the average values of rolling, yawing, and pitching deflection angles, respectively.As shown in Figures8 and 9, the IGC method designed by the adaptive dynamic surface method can well complete the falling point attack in the dive phase in the case of variable span and invariable span. The miss distance of the variable span missile is less than that of the variable span missile. In Figure 9, the total flight time and the lateral distance error of the trace of the variable span are less than those in the invariable span.Figure 8
3D trajectory curve in the dive phase. Inva_sp: invariable span; Va-sp: variable span.Figure 9
Curves representing the missile’s position projected onto the launch coordinate frame. Inva_sp: invariable span; Va-sp: variable span.As shown in Figure10, the trajectory height of the variable span in the dive phase is lower than that of the invariable span. The velocity and Mach number are larger than those of the invariable span. This is because the variable extension characteristic can adjust the lift-to-drag ratio of the missile, so that the resistance and velocity loss is smaller. At the same time, the terminal angular of the variable span missile is closer to the predetermined value than that of the invariable span missile which indicates that it is easier to realize the terminal angular constraint by the rapid change of lift due to the variable span. As shown in Figure 11, the final dynamic pressure of the variable span missile is much larger than that of the invariable span, which is caused by the difference in the final velocity. It demonstrates that the variable span characteristics are beneficial to achieving a fast and high-velocity attack in the dive phase and reduce the flight time.Figure 10
Curves representing the height, velocity, and terminal angle. Inva_sp: invariable span; Va-sp: variable span.Figure 11
Curves representing the Mach number, dynamic pressure, and air density. Inva_sp: invariable span; Va-sp: variable span.It can be seen from Figure12 that the attack angle of the variable span in the dive phase is gentler than the invariant span and the total average value is smaller. Because the required change of the lift is performed by the change of span and attack angle, the change in the attack angle is gentler than the invariant span and the amplitude is also smaller. Because of the influence of the span variant on the aerodynamic coefficient, the changing amplitude of the side-slip angle under the variable span is smaller than the invariant span, while the changes in the bank angles of the two are more similar. It can be seen from Figure 13 that the rolling rate of the variable span has little difference from the invariable span while the yawing rate and the pitching rate have smaller amplitudes and smaller oscillation. As shown in Figure 14, the variation angle of the pitch fin deflections of the variable span is smaller than that of the invariant span which indicates that the variable span reduces the workload for the attitude control system. Figure 15 shows the morphing rate of the variable span missile in the dive phase. It can be seen that the variation of the span in the initial stage is more severe while in the middle and later stages the variation of the span is gentler, but the amplitude is increased. The combination of the change of span and attack angle achieves the change of lift request in the whole mission.Figure 12
Curves representing the attack, sideslip, and bank angles. Inva_sp: invariable span; Va-sp: variable span.Figure 13
Curves representing the rolling, yawing, and pitching rates. Inva_sp: invariable span; Va-sp: variable span.Figure 14
Curves representing the fin deflections. Inva_sp: invariable span; Va-sp: variable span.Figure 15
Curve representing the morphing rate. Inva_sp: invariable span; Va-sp: variable span.


### 4.2. Robustness Analysis of the IGC Method for the Variable Span Missile

In order to verify the robust performance of the integrated guidance control method with variation of span, the influence of the aerodynamic coefficient uncertainty and the atmospheric density uncertainty of the system model is considered in the 6DOF simulation verification model. The distributions of the above uncertainties are as follows.(121)CD=1+ν1C¯D,CL=1+ν2C¯L,CN=1+ν3C¯N,mx=1+ν4m¯x,my=1+ν5m¯y,mz=1+ν6m¯z,ρ=1+ν7ρ¯,ν=ν1ν2⋯ν7T,whereC¯D, C¯L, C¯N, m¯x, m¯y, m¯z, and ρ¯ are, respectively, the nominal aerodynamic coefficient and atmospheric density. ν is the uncertainty parameter vector and satisfies
(122)νi∼N0,0.152,i=1,2,…,6,ν7∼N0,0.22.The initial conditions and parameter settings are the same as those in Table3 and Table 4. After 100 times of Monte Carlo simulation experiments, the results of the variable span missile are as follows.Table 4
Parameter optimization value of the IGC method for the invariable span missile.
Name
Value
Name
Value
Name
Value
k01
3.22
k11
10.21
k21
2.14
k02
0.71
k12
2.14
k22
2.41
ε01
0.23
k13
1.75
k23
5.60
ε02
0.51
τ11
0.41
τ21
0.74
d01
0.13
τ12
0.56
τ22
0.44
d02
0.21
υ1
40.21
τ23
0.45
kF
0.35
μ1
0.75
υ2
75.30
μ2
0.12Table 5
Results of the adaptive dynamic surface IGC method.
Terms
Va_span
Inva_span
dmissm
30.86
43.21
γD∘
60.47
57.63
Vtm/s
2617
2518
Ts
45.05
46.33
QtkPa
3978
3660
δ¯z∘
6.1843
9.9400
δ¯y∘
0.1269
0.0880
δ¯x∘
0.0605
0.0579Table 6
Results of Monte Carlo simulation experiments.
Terms
Variable span
Invariable span
CEPm
51.31
77.20
dmissm
μ
50.61
74.76
σ
55.23
80.29
δ¯z∘
μ
6.2196
9.9341
σ
1.4145
1.4006
δ¯y∘
μ
0.0344
0.0478
σ
0.0056
0.0066
δ¯z∘
μ
0.0370
0.0476
σ
0.0076
0.0103In Table6, μ represents the mean value and σ represents the standard deviation. CEP in Table 6 is defined as the radius of the circle in which the center is the target point and the missile has a 50% probability of hitting the circle. The CEP could be calculated as follows.
(123)CEP=0.615σN+0.562σE,σN<σE,1.177σNσE,σN=σE,0.615σE+0.562σN,σN>σE,where σN and σE are the standard deviations of the north and east, respectively.From Figures16–26, it can be seen that the variable span missile can better meet the robust design requirements of the IGC method under the influence of the uncertainty. The side-slip angle is maintained in a small variation range in the dive phase and the attitude remains bounded which indicates that the designed IGC law can keep the state of the missile stable. In the case of variable span, the mean value of the miss distance is about 50.61 m, the standard deviation is about 55.23 m, and the CEP value is 51.31 m; these variables are smaller than the corresponding value of the invariable span. It indicates that the variable span missile has stronger robustness for the uncertainty aerodynamic caused by the aerodynamic coefficient and atmospheric density. The reason is that the change in the attack angle and span can better adjust the lift requested. In Figure 21, under the uncertainties, the change in morphing rate shows a large difference indicating that span variants can better offset the impact of uncertainties. The mean values of the pitching, yawing, and rolling fin deflections are 6.2196∘, 0.0344∘, and 0.0370∘ with variable span which are less than those of the invariable span, respectively. It shows that fin deflection usage of the variable span is smaller than that of the invariant span, but the standard deviation of the two cases is similar.Figure 16
3D trajectory curve in the dive phase.Figure 17
Curves representing the missile’s position projected onto the launch coordinate frame.Figure 18
Curves representing the height, velocity, and terminal angle.Figure 19
Curves representing the attack, sideslip, and bank angles.Figure 20
Curves representing the fin deflections.Figure 21
Curve representing the morphing rate.Figure 22
3D trajectory curve in the dive phase.Figure 23
Curves representing the missile’s position projected onto the launch coordinate frame.Figure 24
Curves representing the height, velocity, and terminal angle.Figure 25
Curves representing the attack, sideslip, and bank angles.Figure 26
Curves representing the fin deflections.


## 4.1. Performance Analysis of IGC Method for a Variable Span Missile

In order to verify the effectiveness of the IGC method with a variable span auxiliary control in the dive phase of the hypersonic missile, the model is numerically simulated and the control group is set up to observe the effect of the IGC method under the condition of variable span and invariable span. The initial state is shown in Table2.Table 2
The initial state of the missile and the latitude and longitude of the target point.
Name
Value
Name
Value
Name
Value
v0 (m/s)
3000
α0 (°)
2
ϕ0 (°)
0
θ0 (°)
0
β0 (°)
0
ψ0 (°)
0
σ0 (°)
0
γV0 (°)
0
ϕT (°)
1
ωx0 (o/s)
0
x (km)
0
ψT (°)
1
ωy0 (°/s)
0
y (km)
30
ξ10
0.40
ωz0 (∘/s)
0
z (km)
0
ξ20
0.55
λF (°)
60During the simulation process, it is necessary to control the amplitude and rate of change of the attack angle command and the rate of change of the bank angle command as follows.(117)0<αv≤15∘/s,α̇v≤5°/s,γ̇Vv≤20°/s.At the same time, the fin deflections of the missile are restricted as follows:(118)δx≤5∘,δy≤20∘,δz≤20∘.The rate of change of the morphing rate of the wing is as follows.(119)ξ̇≤30%30%/s.Firstly, the control parameters in Eq. (65) are optimized based on the genetic algorithm. The miss distance of the missile, the deviation of the terminal angular, the jitter of the attitude rates, and the jitter of fin deflections are optimized as fitness functions. The attitude rates and the jitter of the fin deflections can be expressed by the sum of the absolute values of their respective changes as
(120)fit=kddmiss+kθxF+kω∫0tdωx+dωy+dωz+kδ∫0tdδx+dδy+dδz,where dmiss is the missile miss distance and t is the flight time of the missile in the dive phase. kd, kθ, kω, and kδare the weight coefficients for measuring the miss distance, the error of the terminal angular, the jitter of rate of the angle, and the jitter of fin deflections. The initial population size is set to 60. After genetic optimization, the optimal controller parameters were obtained. The average fitness value changes with the genetic algebra as shown in Figures 6 and 7. The optimized parameter values are shown in Tables 2 and 3.Figure 6
The fitness value with genetic algebra of the variable span missile.Figure 7
The fitness value with the genetic algebra of the invariable span missile.Table 3
Parameter optimization value of the IGC method for variable span missile.
Name
Value
Name
Value
Name
Value
k01
2.56
k11
10.12
k21
2.25
k02
0.51
k12
3.12
k22
2.78
ε01
0.62
k13
2.41
k23
5.32
ε02
0.33
τ11
0.71
τ21
0.44
d01
0.44
τ12
0.50
τ22
0.22
d02
0.27
υ1
50.11
τ23
1.3
kF
0.33
μ1
2.36
υ2
47.29
τξ
0.21
μ2
1.2The simulation results are as follows.Compare the simulation results obtained by the variable span missile with the invariable span missile as follows.In Table5, γD is the final terminal angular. Vt is the final velocity. T is the total flight time. Qt is the final dynamic pressure. δ¯i,i=x,y,z are the average values of rolling, yawing, and pitching deflection angles, respectively.As shown in Figures8 and 9, the IGC method designed by the adaptive dynamic surface method can well complete the falling point attack in the dive phase in the case of variable span and invariable span. The miss distance of the variable span missile is less than that of the variable span missile. In Figure 9, the total flight time and the lateral distance error of the trace of the variable span are less than those in the invariable span.Figure 8
3D trajectory curve in the dive phase. Inva_sp: invariable span; Va-sp: variable span.Figure 9
Curves representing the missile’s position projected onto the launch coordinate frame. Inva_sp: invariable span; Va-sp: variable span.As shown in Figure10, the trajectory height of the variable span in the dive phase is lower than that of the invariable span. The velocity and Mach number are larger than those of the invariable span. This is because the variable extension characteristic can adjust the lift-to-drag ratio of the missile, so that the resistance and velocity loss is smaller. At the same time, the terminal angular of the variable span missile is closer to the predetermined value than that of the invariable span missile which indicates that it is easier to realize the terminal angular constraint by the rapid change of lift due to the variable span. As shown in Figure 11, the final dynamic pressure of the variable span missile is much larger than that of the invariable span, which is caused by the difference in the final velocity. It demonstrates that the variable span characteristics are beneficial to achieving a fast and high-velocity attack in the dive phase and reduce the flight time.Figure 10
Curves representing the height, velocity, and terminal angle. Inva_sp: invariable span; Va-sp: variable span.Figure 11
Curves representing the Mach number, dynamic pressure, and air density. Inva_sp: invariable span; Va-sp: variable span.It can be seen from Figure12 that the attack angle of the variable span in the dive phase is gentler than the invariant span and the total average value is smaller. Because the required change of the lift is performed by the change of span and attack angle, the change in the attack angle is gentler than the invariant span and the amplitude is also smaller. Because of the influence of the span variant on the aerodynamic coefficient, the changing amplitude of the side-slip angle under the variable span is smaller than the invariant span, while the changes in the bank angles of the two are more similar. It can be seen from Figure 13 that the rolling rate of the variable span has little difference from the invariable span while the yawing rate and the pitching rate have smaller amplitudes and smaller oscillation. As shown in Figure 14, the variation angle of the pitch fin deflections of the variable span is smaller than that of the invariant span which indicates that the variable span reduces the workload for the attitude control system. Figure 15 shows the morphing rate of the variable span missile in the dive phase. It can be seen that the variation of the span in the initial stage is more severe while in the middle and later stages the variation of the span is gentler, but the amplitude is increased. The combination of the change of span and attack angle achieves the change of lift request in the whole mission.Figure 12
Curves representing the attack, sideslip, and bank angles. Inva_sp: invariable span; Va-sp: variable span.Figure 13
Curves representing the rolling, yawing, and pitching rates. Inva_sp: invariable span; Va-sp: variable span.Figure 14
Curves representing the fin deflections. Inva_sp: invariable span; Va-sp: variable span.Figure 15
Curve representing the morphing rate. Inva_sp: invariable span; Va-sp: variable span.


## 4.2. Robustness Analysis of the IGC Method for the Variable Span Missile

In order to verify the robust performance of the integrated guidance control method with variation of span, the influence of the aerodynamic coefficient uncertainty and the atmospheric density uncertainty of the system model is considered in the 6DOF simulation verification model. The distributions of the above uncertainties are as follows.(121)CD=1+ν1C¯D,CL=1+ν2C¯L,CN=1+ν3C¯N,mx=1+ν4m¯x,my=1+ν5m¯y,mz=1+ν6m¯z,ρ=1+ν7ρ¯,ν=ν1ν2⋯ν7T,whereC¯D, C¯L, C¯N, m¯x, m¯y, m¯z, and ρ¯ are, respectively, the nominal aerodynamic coefficient and atmospheric density. ν is the uncertainty parameter vector and satisfies
(122)νi∼N0,0.152,i=1,2,…,6,ν7∼N0,0.22.The initial conditions and parameter settings are the same as those in Table3 and Table 4. After 100 times of Monte Carlo simulation experiments, the results of the variable span missile are as follows.Table 4
Parameter optimization value of the IGC method for the invariable span missile.
Name
Value
Name
Value
Name
Value
k01
3.22
k11
10.21
k21
2.14
k02
0.71
k12
2.14
k22
2.41
ε01
0.23
k13
1.75
k23
5.60
ε02
0.51
τ11
0.41
τ21
0.74
d01
0.13
τ12
0.56
τ22
0.44
d02
0.21
υ1
40.21
τ23
0.45
kF
0.35
μ1
0.75
υ2
75.30
μ2
0.12Table 5
Results of the adaptive dynamic surface IGC method.
Terms
Va_span
Inva_span
dmissm
30.86
43.21
γD∘
60.47
57.63
Vtm/s
2617
2518
Ts
45.05
46.33
QtkPa
3978
3660
δ¯z∘
6.1843
9.9400
δ¯y∘
0.1269
0.0880
δ¯x∘
0.0605
0.0579Table 6
Results of Monte Carlo simulation experiments.
Terms
Variable span
Invariable span
CEPm
51.31
77.20
dmissm
μ
50.61
74.76
σ
55.23
80.29
δ¯z∘
μ
6.2196
9.9341
σ
1.4145
1.4006
δ¯y∘
μ
0.0344
0.0478
σ
0.0056
0.0066
δ¯z∘
μ
0.0370
0.0476
σ
0.0076
0.0103In Table6, μ represents the mean value and σ represents the standard deviation. CEP in Table 6 is defined as the radius of the circle in which the center is the target point and the missile has a 50% probability of hitting the circle. The CEP could be calculated as follows.
(123)CEP=0.615σN+0.562σE,σN<σE,1.177σNσE,σN=σE,0.615σE+0.562σN,σN>σE,where σN and σE are the standard deviations of the north and east, respectively.From Figures16–26, it can be seen that the variable span missile can better meet the robust design requirements of the IGC method under the influence of the uncertainty. The side-slip angle is maintained in a small variation range in the dive phase and the attitude remains bounded which indicates that the designed IGC law can keep the state of the missile stable. In the case of variable span, the mean value of the miss distance is about 50.61 m, the standard deviation is about 55.23 m, and the CEP value is 51.31 m; these variables are smaller than the corresponding value of the invariable span. It indicates that the variable span missile has stronger robustness for the uncertainty aerodynamic caused by the aerodynamic coefficient and atmospheric density. The reason is that the change in the attack angle and span can better adjust the lift requested. In Figure 21, under the uncertainties, the change in morphing rate shows a large difference indicating that span variants can better offset the impact of uncertainties. The mean values of the pitching, yawing, and rolling fin deflections are 6.2196∘, 0.0344∘, and 0.0370∘ with variable span which are less than those of the invariable span, respectively. It shows that fin deflection usage of the variable span is smaller than that of the invariant span, but the standard deviation of the two cases is similar.Figure 16
3D trajectory curve in the dive phase.Figure 17
Curves representing the missile’s position projected onto the launch coordinate frame.Figure 18
Curves representing the height, velocity, and terminal angle.Figure 19
Curves representing the attack, sideslip, and bank angles.Figure 20
Curves representing the fin deflections.Figure 21
Curve representing the morphing rate.Figure 22
3D trajectory curve in the dive phase.Figure 23
Curves representing the missile’s position projected onto the launch coordinate frame.Figure 24
Curves representing the height, velocity, and terminal angle.Figure 25
Curves representing the attack, sideslip, and bank angles.Figure 26
Curves representing the fin deflections.


## 5. Conclusion

In this paper, the IGC method for hypersonic variable span missiles under variable auxiliary control is studied. The physical model and deformation mode of the variable span missile are studied and the aerodynamic model is established by software. The aerodynamic model contains the influence of the morphing process on aerodynamics. In the case of fully considering the span variant of the wing, the motion model of the variable span missile is established which includes the influence of additional force and additional moment caused by the span variant. An IGC design model in the dive phase with terminal angular constraint is established. Based on the adaptive dynamic surface method, the IGC method is designed. At the same time, the stability of the system is proved by the Lyapunov theorem. By utilizing the variable span characteristics of the wing, the synergistic change of the span and the attack angle is completed to realize the control of the lift and achieve rapid and stable control of the centroid motion.The parameters of the IGC method for span variant auxiliary control established above are optimized. The simulations of the variable span and invariable span are set for comparison. Results show that the adaptive dynamic surface method can well complete the motion mission of the dive phase of the hypersonic variable span missile. At the same time, the variable span missile has a smaller miss distance and flight time than the invariable span missile and the terminal angular is closer to the design value. The variations of the attack angle pitching fin deflection angle are smaller. The Monte Carlo simulation method is used to analyze the robustness of the IGC method for variable span missiles. Compared with the invariable span of the exhibition, the results show that the CEP of the variable span missile is less than that of the invariable span, which indicates that the robustness of the variable span missile is stronger than that of the invariable span missile.In addition, in the study of the IGC method of hypersonic variable span missile, the accuracy of the hypersonic aerodynamic model needs more attention in the subsequent research.



---

*Source: 6413410-2019-05-05.xml*