Command filtered back-stepping missile integrated guidance and autopilot based on extended state observer

Abstract

This article focuses on the missile integrated guidance and autopilot design in the end-game phase with control input saturation. First, the nonlinear integrated guidance and autopilot model is developed with third actuator dynamics, where the control surface deflection has magnitude constraint. Second, three nonlinear extended state observers are used to estimate the aerodynamics coefficients’ uncertainties and unmatched time-varying disturbances. Third, a command filtered controller is designed step by step with linear and nonlinear sliding surfaces to improve the terminal performance. In the process, different command filters are implemented to avoid the influences of disturbances and repetitive derivation and solve the problem of unknown control direction. The stability of closed-loop system is proved by Lyapunov theory, and the principles abided by the controller parameters are concluded through the proof. Finally, 6-degree-of-freedom numerical simulations are presented to show the feasibility and validity of the proposed controller.

Introduction

Since 1980s, there have been lots of studies in the area of missile integrated guidance and control (IGC). In the conventional design, the missile guidance and control system is treated as two separated processes based on different operation frequencies. The outer guidance system creates acceleration or angle of attack (AOA) command, and the inner control system tracks it. When the two control loops are combined, the original performance objectives are lost and must be recovered. Thus, the resulting iterative design may not produce an optimal overall system. Additionally, under the condition of high speed, imprecisely known aerodynamics, complicated uncertainties, and external disturbances, the missile dynamics is obviously characterized by nonlinearities. The approaches which involve linearization about a set of equilibrium conditions or trim points within the flight envelop suffer from several disadvantages. Comparatively, the IGC design has provided better solutions in low-cost, modular growth, design flexibility, simple logistics and attracted great interest of researchers.

In Williams et al., ¹ the scheme of tactical missile IGC system is designed, where the information produced by inertial sensors of guidance system are used for the attitude control system. Because the accuracy of sensors in guidance system is much greater than that of autopilot, the IGC can reduce cost and improve the entire guidance and control system performance. In Evers and Cloutier, ² the discussion of IGC scheme goes a step further with three optimal control laws for a tactical missile. In Xin et al. ³ and Vaddi and Menon, ⁴ new optimal control methods are proposed to effectively design IGC system for missiles. However, the weighting matrices have too many elements, while the actuator dynamics is assumed to be sufficiently fast and is not modeled in the IGC development. In comparison, sliding mode control (SMC) is a robust method with not complicated structure for nonlinear control issues. It has been extensively used in IGC design. In Shima et al., ⁵ a brief review of SMC is made at first. A simple first-order actuator dynamics is considered in model development. Then, the missile IGC design is compared with two different separated designs. It is indicated that the inherent instability of decoupled guidance and control loops in terminal phase is postponed by the integrated design, so the interception accuracy is greatly improved. Besides, the control input chasing is depressed using zero-error-miss (ZEM) sliding surface. In Harl et al., ⁶ predicted-impact-point (PIP) heading error is used to make sliding surface, and terminal second-order SMC is designed to achieve the convergence in a finite time without chattering. In Shima et al. ⁷ and Tournes and Shtessel, ⁸ the IGC dynamics is built with two first-order actuators for dual-control missiles, and the controller complexity becomes staggering. Nevertheless, the true missile dynamics has high order for second- or third-order actuator dynamics. For high-order nonlinear system model in literatures,^9–12 back-stepping technique is a useful tool, which is always combined with SMC and disturbance observer. In Huang et al., ⁹ the nonlinear disturbance observer is implemented to estimate the nonlinearities, uncertainties, and disturbances in the system dynamics, thus the decrease in undesired chattering in control is achieved and the robustness of closed-loop system is enhanced. In Guo et al., ¹⁰ the extended state observers (ESOs) are used to estimate indirectly measured states and various parametric uncertainties. The nonlinear ESO in Xia et al. ¹¹ has better performance under complex uncertainties and measurement noises than linear ESO in Shao and Wang;^12,13 however, it has more parameters to tune. In Sun et al., ¹⁴ velocity tracking error is used to design a composite-error-based ESO in a feedback form, which makes estimation and tracking errors smaller without high gains.

To simplify the implementation and cancel out the system noises, command filter is introduced. In Farrell et al., ¹⁵ command filtered back-stepping is proposed to offer a means to get the time derivatives of the pseudo control signals. In Huo and colleagues,^16,17 low-pass filter is used to construct the derivative of pseudo control input. It solves the problem of “explosion of complexity” caused by the repeated differentiations of the pseudo control signals in dynamic surface control (DSC). In Pan et al., ¹⁸ the second-order command filters instead of the first-order filters are applied in DSC. With the help of command filters, the performance of back-stepping control scheme significantly improves in stability and steady-state tracking accuracy, while the analysis is made in detail in Pan and Yu. ¹⁹ In Pan and Yu, ²⁰ the stability of commander filtered back-stepping control is further improved by composite learning. In Farrell et al. ²¹ and Xu et al., ²² a second-order command filter is designed to impose magnitude and rate limitation on the control input. In Chwa, ²³ directly differentiating the pseudo control command with respect to time is avoided and the global uniform ultimate boundlessness of the tracking errors is guaranteed in the presence of the input constraints. In Erdos et al., ²⁴ low-pass filters of the adaptive control scheme guarantee the fast adaptive performance and robustness of the missile integrated guidance and autopilot (IGA) system. Deep discussion about control input constraint, finite time convergence in IGC design arises rapidly in Wang and Wang^25,26 and Shi et al. ²⁷ In Wen et al., ²⁸ first-order auxiliary dynamics is developed in addition to the system model to deal with the input constraints. In Chen et al.,^29,30 Nussbaum function is introduced to compensate for the nonlinear term arising from input saturation and solve the problem of unknown control direction. In Chen et al., ³⁰ a novel Nussbaum gain is proposed for multiple unknown actuator directions and time-varying nonlinearities.

In summary, IGC study should be more focused on the model development with actuator dynamics, coefficients’ uncertainties, and control input constraints. Then, the application is able to be studied and discussed more practically and effectively. Motivated by (1) the IGC design, (2) ESO-based disturbance estimate, and (3) command filtered back-stepping control in the presence of control input constraint, this work proposes a novel missile IGA scheme.

The article is organized as follows: in section “Introduction,” the integrated design and related control techniques are introduced. In section “Missile IGA model,” the missile flight dynamics in end-game phase with third-order actuator model is developed, and the control surface deflection has constraint of magnitude saturation. In section “ESO,” ESO is designed to synthetically estimate the uncertainty and time-varying disturbance, and three ESOs are used in different channels of the IGA model. In section “Command filtered controller,” the command filtered back-stepping technique is implemented based on the estimations of ESOs. Command filters are used to get the derivate of pseudo inputs produced by the sliding mode controllers, and Nussbaum function is introduced to solve the problem of unknown control direction when differentiating the saturated control input. In section “Stability analysis,” the stability proof of closed-loop system is given by Lyapunov theory, and the principles followed by the controller parameters are concluded. In section “Numerical simulation,” 6-degree-of-freedom (DOF) numerical simulations verify the proposed design with two scenarios, and Monte Carlo simulations are made to test the performance under initial flight states bias and measurement noises. Finally, section “Conclusion” concludes the article.

Missile IGA model

In this section, the missile IGA model in the end-game phase is developed. In Figure 1, the coordinate systems are built without considering the earth rotation. The aerodynamics is built in the velocity coordinate system. V is the velocity of mass center, R is the distance between the missile and the target, q is the angle between line of sight (LOS) and local horizon, $\gamma$ is the flight path angle, $\alpha$ is the AOA, and ${\delta }_z$ is the control surface deflection. The angle directions shown in Figure 1 are positive.

Figure 1.

The planar geometry in end-game phase.

The state vector is chosen as $\mathit{x}=[x_1,x_2,x_3,x_4{]}^T=[\stackrel{·}{q},\alpha ,{\omega }_z,{\delta }_z{]}^T$ , and the control surface deflection is saturated as $y_{\mathrm{act}}=\mathrm{sat}(x_4)$ . Then the state-space form of IGC system can be depicted as follows

$${\stackrel{·}{x}}_1=f_1(x_1)+a{x}_2+d_1$$

(1a)

$${\stackrel{·}{x}}_2=f_2(x_2)+x_3+d_2$$

(1b)

$${\stackrel{·}{x}}_3=f_3(x_2,x_3)+b{y}_{\mathrm{act}}+d_3$$

(1c)

where $d_i(t),i=1,2,3$ are the unmatched time-varying disturbances in different channels. The nonlinear functions $f_i(x),i=1,2,3$ and $a,b$ , which contain uncertain aerodynamics coefficients ${\overline{L}}_{\alpha },{\overline{L}}_o,{\overline{M}}_{\alpha },{\overline{M}}_{\omega },{\overline{M}}_{{\delta }_z}$ , are depicted as follows

$$\begin{array}{c}{f}_1(x)=-\frac{x_1\stackrel{·}{R}}{R}+\frac{V}{R}\cos (q-\gamma )(x_1+\frac{g}{V}\cos \gamma -{\overline{L}}_o),\hfill \\ a=-\frac{V}{R}\cos (q-\gamma ){\overline{L}}_{\alpha }\hfill \end{array}$$

(2a)

$$f_2(x)=\frac{g}{V}\cos \gamma -{\overline{L}}_{\alpha }{x}_2-{\overline{L}}_o$$

(2b)

$$f_3(x)={\overline{M}}_{\alpha }{x}_2+{\overline{M}}_{\omega }{x}_3,b={\overline{M}}_{{\delta }_z}$$

(2c)

The lift force can be simplified as two parts: $L=L_{\alpha }\alpha +L_o$ when the velocity is high. Furthermore, $L_{\alpha }=Q{S}_{\mathrm{ref}}{C}_L^{\alpha }$ , where $C_L^{\alpha }$ is the lift force derivative with respect to AOA, $Q=(1/2)\rho V^2$ is the dynamic pressure, and $S_{\mathrm{ref}}$ is the reference size; $L_o=Q{S}_{\mathrm{ref}}{C}_{\mathrm{Lo}}$ is the other factors’ contribution to lift. Then ${\overline{L}}_o=L_o/\mathrm{mV},{\overline{L}}_{\alpha }=L_{\alpha }/\mathrm{mV}$ in equations (2a) and (2b) are the unified parts of lift, and $m$ is the mass. The aerodynamics moment is also described as two parts $M_z=M_o+M_{{\delta }_z}{\delta }_z$ , where $m_z^{{\delta }_z}$ is the moment derivative with respect to control surface deflection and $l_{\mathrm{ref}}$ is the reference length in $M_{{\delta }_z}=Q{S}_{\mathrm{ref}}{l}_{\mathrm{ref}}{m}_z^{{\delta }_z}$ , and $M_o$ is the other factors contribution to moment. Both are unified by the moment of inertia with respect to body z-axis $I_z:{\overline{M}}_o=M_o/I_z,{\overline{M}}_{{\delta }_z}=M_{{\delta }_z}/I_z$ in equation (2c).

The actuator dynamics is considered as a first-order inertial element with time constant ${\tau }_{\mathrm{act}}$ and a second-order element with damping rate ${\xi }_{\mathrm{act}}$ , natural frequency ${\omega }_{n,\mathrm{act}}$ , and constant gain $K_{\mathrm{act}}$ . Then the transfer function is expressed as follows

$$\frac{{\delta }_z}{{\delta }_{\mathrm{zc}}}=\frac{1}{{\tau }_{\mathrm{act}}s+1}\frac{K_{\mathrm{act}}}{s^2+2{\xi }_{\mathrm{act}}{\omega }_{n,\mathrm{act}}s+{\omega }_{n,\mathrm{act}}^2}$$

(3)

The state-space model of actuator can be depicted as follows

$${\stackrel{·}{x}}_4=-\frac{1}{{\tau }_{\mathrm{act}}}{x}_4+\frac{1}{{\tau }_{\mathrm{act}}}{x}_{4\tau }$$

(4a)

$${\stackrel{·}{x}}_{4\tau }=x_5$$

(4b)

$${\stackrel{·}{x}}_5=f_{\mathrm{act}}(x_{4\tau },x_5)+b_{\mathrm{act}}u+d_{\mathrm{act}}$$

(4c)

where $f_{\mathrm{act}}(x_{4\tau },x_5)$ is a linear function and $d_{\mathrm{act}}(t)$ is the time-varying disturbance in actuator dynamics. The IGA model is developed by equations (1) and (4) with state vector $\mathit{X}=[x_1,x_2,x_3,x_4,x_{4\tau },x_5{]}^T$ , and the control input is $u={\delta }_{\mathrm{zc}}$ .

The nonlinear saturation characteristic of the actuator in equation (1c) can be modeled as

$$\mathrm{sat}(x_4)=\{\begin{array}{c}{{\delta }_{z,}}_{\max }\mathrm{sgn}({\delta }_z),\left|{\delta }_z\right|\ge {{\delta }_{z,}}_{\max }\hfill \\ \multicolumn{1}{c}{{\delta }_z,\left|{\delta }_z\right|<{{\delta }_{z,}}_{\max }}\end{array}$$

(5)

where ${\delta }_{z,\max }$ is the magnitude constraint of control surface deflection. The saturated signal can be approximated by a smooth function defined as

$$g(x_4)={\delta }_{z,\max }\tanh \left(\frac{{\delta }_z}{{\delta }_{z,\max }}\right)={\delta }_{z,\max }\frac{e^{\frac{{\delta }_z}{{\delta }_{z,\max }}}-e^{-\frac{{\delta }_z}{{\delta }_{z,\max }}}}{e^{\frac{{\delta }_z}{{\delta }_{z,\max }}}+e^{-\frac{{\delta }_z}{{\delta }_{z,\max }}}}$$

(6)

then $\mathrm{sat}(x_4)=g(x_4)+o$ , where $o$ is the small approximation error.

ESO

The disturbance estimation technique is particularly crucial for disturbance attenuation. To compensate the IGA system under uncertainties and unmatched disturbances, the technique of ESO is implemented. It is developed as follows:

First, considering the following system

$${\stackrel{·}{x}}_{\mathrm{sys}}(t)=f_{\mathrm{sys}}(x)+b_{\mathrm{sys}}{u}_{\mathrm{sys}}(t)+d_{\mathrm{sys}}(t)$$

(7)

The uncertain nonlinearity and time-varying disturbance are considered as a whole ${\mathrm{\Delta }}_{\mathrm{sys}}=f_{\mathrm{sys}}(x)+d_{\mathrm{sys}}(t)$ to be estimated. If defining ${\mathrm{\Delta }}_{e\mathrm{s}o}$ and $x_{\mathrm{eso}}$ as the states of ESO, which estimate, respectively, ${\mathrm{\Delta }}_{\mathrm{sys}}$ and $x_{\mathrm{sys}}$ , then nonlinear ESO is designed as follows

$$\{\begin{array}{c}{\stackrel{·}{x}}_{\mathrm{eso}}={\mathrm{\Delta }}_{\mathrm{eso}}-{\beta }_1(x_{\mathrm{eso}}-x_{\mathrm{sys}})+b_{\mathrm{sys}}{u}_{\mathrm{sys}}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathrm{\Delta }}}_{\mathrm{eso}}=-{\beta }_2\mathrm{fal}(x_{\mathrm{eso}}-x_{\mathrm{sys}},\mu ,{\delta }_{\mathrm{fal}})}\end{array}$$

(8)

where $\mathrm{fal}(·)$ is a nonlinear function with good performance in dealing with measurement noises

$$\begin{array}{c}\mathrm{fal}(x_{\mathrm{eso}}-x_{\mathrm{sys}},\mu ,{\delta }_{\mathrm{fal}})\hfill \\ =\{\begin{array}{c}{|x_{\mathrm{eso}}-x_{\mathrm{sys}}|}^{\mu }\mathrm{sgn}(x_{\mathrm{eso}}-x_{\mathrm{sys}}),\mathrm{if}|x_{\mathrm{eso}}-x_{\mathrm{sys}}|>{\delta }_{\mathrm{fal}}\hfill \\ \multicolumn{1}{c}{\frac{x_{\mathrm{eso}}-x_{\mathrm{sys}}}{{\delta }^{1-\mu }},\mathrm{else}}\end{array}\hfill \end{array}$$

(9)

The estimator error of the whole disturbance $e_{\mathrm{\Delta }}={\mathrm{\Delta }}_{\mathrm{eso}}-{\mathrm{\Delta }}_{\mathrm{sys}}$ is determined by ${\beta }_1,{\beta }_2,\mu$ , and ${\delta }_{\mathrm{fal}}$ . The basic principle is to choose ${\beta }_1,{\beta }_2>0$ , $0<\mu <1$ , ${\delta }_{\mathrm{fal}}>0$ , and ${\beta }_1>>{\beta }_2$ . A smaller $\mu$ can make the estimate signal more smooth.

In the missile IGA design, three ESOs (ESO1, ESO2, and ESO3) are used to estimate the unmatched uncertain nonlinearities and time-varying disturbances ${\mathrm{\Delta }}_i=f_i(x)+d_i(t),i=1,2,3$ in equation (1); the results are respectively ${\mathrm{\Delta }}_{\mathrm{esoi}},i=1,2,3$

Command filtered controller

In this section, based on the estimate states by ESOs, the command filtered back-stepping sliding mode controller is designed to achieve LOS rate convergence:

• Step 1: The first sliding surface is chosen as $s_1=x_1$ . The pseudo control input ${\overline{x}}_2$ is designed according to proportional reaching law with a positive constant $k_1$

$${\overline{x}}_2=-\frac{1}{a}(k_1{s}_1+{\mathrm{\Delta }}_{\mathrm{eso}1})$$

(10)

If directly differentiating ${\overline{x}}_2$ with respect to time, the result is complex. So a low-pass command filter is used to get its derivative

$${\tau }_1{\stackrel{·}{x}}_{2r}+x_{2r}={\overline{x}}_2,x_{2r}(0)={\overline{x}}_2(0)$$

(11)

In equation (11), ${\tau }_1$ is a positive constant. The error produced by the low-pass filter is defined as ${\lambda }_1=x_{2r}-{\overline{x}}_2$ , and ${\lambda }_1$ can be sufficiently small when ${\tau }_1$ is appropriately set.

• Step 2: Defining the second sliding surface $s_2=x_2-x_{2r}$ , the pseudo control input ${\overline{x}}_3$ is obtained by the following control law

$${\overline{x}}_3=-k_2{s}_2-{\mathrm{\Delta }}_{\mathrm{eso}2}+{\stackrel{·}{x}}_{2r}$$

(12)

Another low-pass filter in the same form of equation (11) is used to avoid directly differentiating ${\overline{x}}_3$

$${\tau }_2{\stackrel{·}{x}}_{3r}+x_{3r}={\overline{x}}_3,x_{3r}(0)={\overline{x}}_3(0)$$

(13)

Then the error produced by the low-pass filter can be expressed as ${\lambda }_2=x_{3r}-{\overline{x}}_3$ under appropriate time constant ${\tau }_2$ .

• Step 3: In this step, there are two sub-steps to finish the design. First, the pseudo control input ${\overline{x}}_4$ is designed to achieve the convergence of $s_3=x_3-x_{3r}$ under the following control law

$${\overline{x}}_4=-k_3{s}_3-{\mathrm{\Delta }}_{\mathrm{eso}3}+{\stackrel{·}{x}}_{3r}$$

(14)

The third low-pass filter is used to get the first-order derivative of ${\overline{x}}_4$

$${\tau }_3{\stackrel{·}{x}}_{4r}+x_{4r}={\overline{x}}_4,x_{4r}(0)={\overline{x}}_4(0)$$

(15)

The error between ${\overline{x}}_4$ and $x_{4r}$ is given as ${\lambda }_3=x_{4r}-{\overline{x}}_4$ under time constant ${\tau }_3$ .

Second, if defining $s_4=g(x_4)-x_{4r}$ , then differentiating it with respect to time

$${\stackrel{·}{s}}_4=\xi {\stackrel{·}{x}}_4-{\stackrel{·}{x}}_{4r}$$

(16)

where $x_4$ is also given by a low-pass filter with time constant ${\tau }_{\mathrm{act}}$

$$\xi =\frac{\partial g(x_4)}{\partial x_4},{\tau }_{\mathrm{act}}{\stackrel{·}{x}}_4+x_4=x_{4c}$$

(17)

Then the proportional reaching law with positive constant $k_4$ is chosen to achieve the convergence of $s_4$ . The control direction is unknown because of $\xi$ , and a Nussbaum function $N(\chi )$ is used to design pseudo control input $x_{4u}$

$$x_{4u}={\tau }_{\mathrm{act}}N(\chi )x_{4\mathrm{uc}}$$

(18)

$$x_{4\mathrm{uc}}=(\frac{\xi }{{\tau }_{\mathrm{act}}}{x}_4+{\stackrel{·}{x}}_{4r}-k_4{s}_4)$$

(19)

The Nussbaum function is defined with the following properties

$$\underset{s\to \infty }{lim}sup\frac{1}{s}\underset{0}{\overset{s}{\int }}N(\zeta )d\zeta =+\infty$$

(20a)

$$\underset{s\to \infty }{lim}inf\frac{1}{s}\underset{0}{\overset{s}{\int }}N(\zeta )d\zeta =-\infty$$

(20b)

According to the properties, the following Nussbaum function is implemented

$$N(\chi )=e^{{\chi }^2}\cos (\chi )$$

(21)

The parameter $\chi$ is an adaptive parameter according to the following principle

$$\stackrel{·}{\chi }={\gamma }_{\chi }(\frac{\xi }{{\tau }_{\mathrm{act}}}{x}_4+{\stackrel{·}{x}}_{4r}-k_4{s}_4)s_4$$

(22)

At the end of Step 3, a third-order differentiator is also used as command filter to get first- and second-order derivatives of the pseudo control input $x_{4u}$

$$\begin{array}{c}{\stackrel{·}{x}}_{\mathrm{third}1}=x_{\mathrm{third}2}\hfill \\ {\stackrel{·}{x}}_{\mathrm{third}2}=x_{\mathrm{third}3}\hfill \\ {\stackrel{·}{x}}_{\mathrm{third}3}=\frac{4}{{\epsilon }^3}(-2^{\frac{5}{3}}({\left|v\right|}^{\frac{1}{3}}\mathrm{sgn}(v))-{\left|{\epsilon }^2{x}_{\mathrm{third}3}\right|}^{\frac{5}{3}}\mathrm{sgn}({\epsilon }^2{x}_{\mathrm{third}3})),\hfill \\ x_{\mathrm{third}1}(0)=x_{4u}(0)\hfill \end{array}$$

(23)

where $v=x_{\mathrm{third}1}-x_{4u}+{|(1/\epsilon )x_{\mathrm{third}2}|}^{9/7}\mathrm{sgn}(\epsilon x_{\mathrm{third}2})$ , and the error between $x_{\mathrm{third}3}$ and $x_{4u}$ can be small enough through choosing suitable $\epsilon$ .

• Step 4: The task of tracking $x_{4u}$ is completed by terminal sliding mode controller. Considering the following second-order dynamics in equation (4)

$${\stackrel{·}{x}}_{4\tau }=x_5$$

(24a)

$${\stackrel{·}{x}}_5=f_{\mathrm{act}}(x_{4\tau },x_5)+\mathrm{\Delta }{f}_{\mathrm{act}}+b_{\mathrm{act}}u(t)+d_{\mathrm{act}}(t)$$

(24b)

In equation (24b), $\mathrm{\Delta }{f}_{\mathrm{act}}$ is the uncertainty of $f_{\mathrm{act}}(x_{4\tau },x_5)$ and $d_{\mathrm{act}}(t)$ is the time-varying disturbance, which are bounded by positive constants $F$ and $D$

$$\left|\mathrm{\Delta }{f}_{\mathrm{act}}\right|\le F,\left|d_{\mathrm{act}}\right|\le D$$

(25)

Defining vector $E=[e_{\mathrm{act}},{\stackrel{·}{e}}_{\mathrm{act}}{]}^T$ with $e_{\mathrm{act}}=x_{4\tau }-x_{4u}$ , its derivative can be written as ${\stackrel{·}{e}}_{\mathrm{act}}={\stackrel{·}{x}}_{4\tau }-x_{\mathrm{third}2}$ . Then the following nonlinear sliding mode surface is designed

$$\sigma =\mathrm{CE}-\mathrm{CP}$$

(26)

In equation (26), $C=[c_1,c_2]$ is a constant vector, and $P=[p,\stackrel{·}{p}{]}^T$ is determined by the following nonlinear function

$$p(t)=\{\begin{array}{c}\begin{array}{c}{e}_{\mathrm{act}}+{\stackrel{·}{e}}_{\mathrm{act}}t+\frac{1}{2}{\stackrel{··}{e}}_{\mathrm{act}}{t}^2-(\frac{10}{T^3}{e}_{\mathrm{act}}+\frac{6}{T^2}{\stackrel{·}{e}}_{\mathrm{act}}+\frac{3}{2T}{\stackrel{··}{e}}_{\mathrm{act}})t^3\hfill \\ \multicolumn{1}{c}{+(\frac{15}{T^4}{e}_{\mathrm{act}}+\frac{8}{T^3}{\stackrel{·}{e}}_{\mathrm{act}}+\frac{3}{2T^2}{\stackrel{··}{e}}_{\mathrm{act}})t^4-(\frac{6}{T^5}{e}_{\mathrm{act}}+\frac{3}{T^4}{\stackrel{·}{e}}_{\mathrm{act}}+\frac{1}{2T^3}{\stackrel{··}{e}}_{\mathrm{act}})t^5,0\le t\le T}\end{array}\hfill \\ \multicolumn{1}{c}{0,t>T}\end{array}$$

(27)

where $T$ is the convergence time of the terminal sliding mode controller.

Finally, the control input is designed as follows

$$\begin{array}{c}\multicolumn{1}{c}{u=-\frac{1}{b_{\mathrm{act}}}(f_{\mathrm{act}}(x)-x_{\mathrm{third}3}-\stackrel{··}{p}+\frac{c_1}{c_2}({\stackrel{·}{e}}_{\mathrm{act}}-\stackrel{·}{p}))}\\ -\frac{1}{b_{\mathrm{act}}}\mathrm{sgn}(c_2\sigma )(F+D)\hfill \end{array}$$

(28)

Stability analysis

Defining $e_1=x_1,e_2=x_2-{\overline{x}}_2,e_3=x_3-{\overline{x}}_3,e_4=g(x_4)-{\overline{x}}_4$ , and differentiating them with respect to time

$${\stackrel{·}{e}}_1={\stackrel{·}{x}}_1=-k_1{e}_1+e_2+{\epsilon }_1$$

(29)

$${\stackrel{·}{e}}_2={\stackrel{·}{x}}_2-{\stackrel{·}{\overline{x}}}_2=-k_2{e}_2+e_3+{\epsilon }_2+k_2{\lambda }_1+{\stackrel{·}{\lambda }}_1$$

(30)

$${\stackrel{·}{e}}_3={\stackrel{·}{x}}_3-{\stackrel{·}{\overline{x}}}_3=-k_3{e}_3+e_4+{\epsilon }_3+k_3{\lambda }_2+{\stackrel{·}{\lambda }}_2$$

(31)

$${\stackrel{·}{e}}_4=x_{4\mathrm{uc}}(\xi N(\chi )-1)-k_4{e}_4+{\stackrel{·}{\lambda }}_3+k_4{\lambda }_3$$

(32)

Besides, the following Lyapunov function is chosen

$$V=\frac{1}{2}{e}_1^2+\frac{1}{2}{e}_2^2+\frac{1}{2}{e}_3^2+\frac{1}{2}{e}_4^2+\frac{1}{2}{\sigma }^2$$

(33)

Differentiating part of equation (33) with respect to time

$$\begin{array}{c}\multicolumn{1}{c}{\sigma \stackrel{·}{\sigma }=\sigma \{c_2(f_{\mathrm{act}}-{\stackrel{··}{x}}_{4u}-\stackrel{··}{p}+\frac{c_1}{c_2}(\stackrel{·}{e}-\stackrel{·}{p}))+c_2\mathrm{bu}+c_2(\mathrm{\Delta }{f}_{\mathrm{act}}+d_{\mathrm{act}})\}}\\ \le \sigma (c_2(x_{\mathrm{third}3}-{\stackrel{··}{x}}_{4u})-c_2\mathrm{sgn}(c_2\sigma )(F+D)+c_2(\mathrm{\Delta }{f}_{\mathrm{act}}+d_{\mathrm{act}}))\hfill \\ \le c_2\sigma (x_{\mathrm{third}3}-{\stackrel{··}{x}}_{4u})-\left|c_2\sigma \right|(F+D)+c_2\sigma (\mathrm{\Delta }{f}_{\mathrm{act}}+d_{\mathrm{act}})\hfill \end{array}$$

(34)

Because the differentiator is set appropriately such that $|x_{\mathrm{third}3}-{\stackrel{··}{x}}_{4u}|\le o_{\mathrm{third}}$ , $o_{\mathrm{third}}$ is small positive constant, and $\mathrm{\Delta }{f}_{\mathrm{act}},d_{\mathrm{act}}$ is bounded. We can find a positive constant $\mathrm{K}$ to yield to

$$\sigma \stackrel{·}{\sigma }\le -\mathrm{K}\left|c_2\sigma \right|$$

(35)

Through differentiating equation (33) and combining equation (34), we have

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{V}\le e_1{\stackrel{·}{e}}_1+e_2{\stackrel{·}{e}}_2+e_3{\stackrel{·}{e}}_3+e_4{\stackrel{·}{e}}_4+\sigma \stackrel{·}{\sigma }}\\ \le -k_1{e}_1^2-k_2{e}_2^2-k_3{e}_3^2-k_4{e}_4^2+e_1{e}_2+e_2{e}_3+e_1{\epsilon }_1\hfill \\ +e_2({\epsilon }_2+{\stackrel{·}{\eta }}_1+k_2{\eta }_1)+e_3({\epsilon }_3+{\stackrel{·}{\eta }}_2+k_3{\eta }_2)\hfill \\ +e_4({\stackrel{·}{\eta }}_3+k_4{\eta }_3)+x_{4\mathrm{uc}}{e}_4[\xi N(\chi )-1]\hfill \end{array}$$

(36)

If defining the following bounds

$$\kappa ={\epsilon }_{\max }+k_{\max }{\eta }_{\max }+{\stackrel{·}{\eta }}_{\max }$$

(37a)

$${\epsilon }_{\max }=\max [\left|{\epsilon }_1\right|,\left|{\epsilon }_2\right|,\left|{\epsilon }_3\right|]$$

(37b)

$$k_{\max }=\max [k_1,k_2,k_3]$$

(37c)

$${\eta }_{\max }=\max [\left|{\eta }_1\right|,\left|{\eta }_2\right|,\left|{\eta }_3\right|]$$

(37d)

$${\stackrel{·}{\eta }}_{\max }=\max [\left|{\stackrel{·}{\eta }}_1\right|,\left|{\stackrel{·}{\eta }}_2\right|,\left|{\stackrel{·}{\eta }}_3\right|]$$

(37e)

According to Young’s inequality, equation (36) can be rewritten as follows

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{V}\le -(k_1-1)e_1^2-(k_2-\frac{3}{2})e_2^2-(k_3-1)e_3^2}\\ -(k_4-1)e_4^2+\frac{1}{2}{\kappa }^2+x_{4\mathrm{uc}}{e}_4[\xi N(\chi )-1]\hfill \end{array}$$

(38)

Remark

If choosing $C=min[(k_1-1),(k_2-(3/2)),(k_3-1),(k_4-1)]$ , then the following inequality is yielded to

$$\stackrel{·}{V}\le -\mathrm{CV}+\frac{1}{2}{\kappa }^2+\frac{\stackrel{·}{\chi }}{{\gamma }_{\chi }}[\xi N(\chi )-1]$$

(39)

Integrating equation (39) directly, then we have

$$\begin{array}{c}V(t)\le V(0)e^{-\mathrm{Ct}}+\frac{{\kappa }^2}{2C}(1-e^{-\mathrm{Ct}})\hfill \\ +\frac{e^{-\mathrm{Ct}}}{{\gamma }_{\chi }}\underset{0}{\overset{t}{\int }}\stackrel{·}{\chi }[\xi N(\chi )-1]e^{C\tau }d\tau \hfill \end{array}$$

(40)

According to the proof in Wang and Wang, ²⁶ ${\int }_0^t\stackrel{·}{\chi }[\xi N(\chi )-1]e^{C\tau }d\tau$ is bounded, and $V(t)$ is bounded, which implied all the error is bounded.

Numerical simulation

In this section, 6-DOF nonlinear numerical simulations are provided to illustrate the control schemes proposed in the previous sections. First, the initial flight condition in end-game phase is set as follows

$$\begin{array}{c}{R}_0=20\mathrm{km},V=1200\mathrm{m}/\mathrm{s},g=9.8\mathrm{m}/{\mathrm{s}}^2,\hfill \\ {\gamma }_0=-10\circ ,{\alpha }_0=0\circ ,q_0=-60\circ \hfill \end{array}$$

(41)

The IGA model is built in two scenarios with the following aerodynamics coefficients based on Table 1 in Appendix 1.

• Scenario 1

$${\overline{L}}_{\alpha }=1,{\overline{L}}_o=-0.1,{\overline{M}}_{\delta }=1,{\overline{M}}_{\alpha }=0.1,{\overline{M}}_{\omega }=-1$$

(42a)

• Scenario 2

$$\begin{array}{c}{\overline{L}}_{\alpha }=1-0.01\alpha ,{\overline{L}}_o=-0.1,{\overline{M}}_{\delta }=1,\hfill \\ {\overline{M}}_{\alpha }=0.1-0.001\alpha ,{\overline{M}}_{\omega }=-1-0.001{\omega }_z\hfill \end{array}$$

(42b)

The third-order actuator dynamics is given by the following transfer function

$$G_{{\delta }_z}(s)=\frac{1}{0.0198s+1}\frac{4621}{s^2+66.34s+4848}$$

(43)

The time-varying disturbances in different channels are given as follows

$$\begin{array}{c}{d}_1=0.001\sin (2t),d_2=0.05\cos (2t),\hfill \\ d_3=0.05\sin (t)\cos (2t),d_{\mathrm{act}}=0.01\sin (t)\hfill \end{array}$$

(44)

Second, the three ESOs are set with parameters as follows

$$\begin{array}{c}{\mathit{\beta }}_1=[30,30,30],{\mathit{\beta }}_2=[5,5,5],\hfill \\ \mathit{\mu }=[0.15,0.15,0.15],{\mathit{\delta }}_{\mathrm{fal}}=[0.01,0.01,0.01]\hfill \end{array}$$

(45)

Third, the first three one-order command filters have time constants ${\tau }_1={\tau }_2={\tau }_3=0.0198\mathrm{s}$ . The third-order differentiator in equation (23) is set with $\epsilon =0.04$ . The control surface deflection is constraint with ${\delta }_{z,\max }=10\circ ,30\circ$ , and the Nussbaum function in equation (22) is set with ${\gamma }_{\chi }={10}^{-9}$ . The proportional reaching law of back-stepping sliding mode controller is designed with the following positive constants

$$k_1=0.75,k_2=0.75,k_3=0.75$$

(46)

Finally, the terminal sliding surface in equation (26) is set with $c_1=5$ and $c_2=1$ , and the bound of uncertainty and time-varying disturbance in equation (28) is given as $F=10$ and $D=3$ .

Comparison simulating

The simulation results in Scenario 1 (S1) under saturated ${\delta }_{\max }=10\circ$ are compared with those by the method of IGA without ESOs in Figure 2. As Figure 2(a) shows, $\stackrel{·}{q}$ finally converges to zero under the saturated control surface deflection. In Figure 2(d), the control surface deflection of the proposed method is well-constrained under saturation. The results indicate that the commander filtered controller in this article handles the nonlinear saturation phenomenon very well. The flight states $\alpha ,{\omega }_z$ in S1 with ESOs change in smaller range than without ESOs.

Figure 2.

Comparison with and without ESOs: (a) curves of $\stackrel{·}{q}$ , (b) curves of $\alpha$ , (c) curves of ${\omega }_z$ , and (d) curves of saturated ${\delta }_z$ .

The performance of three nonlinear ESOs in Scenario 2 (S2) is shown in Figures 3–5 through comparison between estimate states and actual states, and estimate synthetical uncertainties and preset synthetical uncertainties. First of all, the tracking errors and estimation errors in the closed-loop control system are uniformly ultimately bounded. It can also be seen that all the ESOs estimate the nonlinearities and disturbances with good accuracy, so the system dynamics is well-compensated. As a result, the controller performs better with ESOs.

Figure 3.

States’ comparison of ESO1: (a) comparison between $x_{\mathrm{eso}1}$ and $\stackrel{·}{q}$ and (b) comparison between ${\mathrm{\Delta }}_{\mathrm{eso}1}$ and ${\mathrm{\Delta }}_1$ .

Figure 4.

States’ comparison of ESO2: (a) comparison between $x_{\mathrm{eso}2}$ and $\alpha$ and (b) comparison between ${\mathrm{\Delta }}_{\mathrm{eso}2}$ and ${\mathrm{\Delta }}_2$ .

Figure 5.

States’ comparison of ESO3: (a) comparison between $x_{\mathrm{eso}3}$ and ${\omega }_z$ and (b) comparison between ${\mathrm{\Delta }}_{\mathrm{eso}3}$ and ${\mathrm{\Delta }}_3$ .

Two scenarios of S1 with equation 42(a) and S2 with equation 42(b) are set to test the performance under large aerodynamics coefficient uncertainties. In equation (42b), ${\overline{L}}_{\alpha }$ in the unified lift force is set with nonlinearity item of AOA, while ${\overline{M}}_{\alpha },{\overline{M}}_{\omega }$ in the unified moment are also nonlinear functions of AOA and pitch angle rate. The comparison results are shown in Figure 6, which indicate that the proposed method can cancel out the uncertainties and disturbances very well. From Figure 6(d), the control surface deflection is saturated in the beginning 4 s, and then works well in the range of $[-10\circ ,10\circ ]$ .

Figure 6.

Comparison between Scenario 1 and Scenario 2: (a) curves of $\stackrel{·}{q}$ , (b) curves of $\alpha$ , (c) curves of ${\omega }_z$ , and (d) curves of saturated ${\delta }_z$ .

The results shown in Figure 7 are comparison simulation with different magnitude constraints of control surface deflection ${\delta }_{\max }=10\circ$ and ${\delta }_{\max }=30\circ$ using aerodynamics coefficients in S1. It can be seen that the proposed controller performs well under both situations. The convergence of $\stackrel{·}{q}$ is faster with large control surface deflection, while $\alpha ,{\omega }_z$ change in larger range and vibration.

Figure 7.

Comparison with different ${\delta }_{zmax}$ in Scenario 1: (a) curves $\stackrel{·}{q}$ of under different ${\delta }_{z,max}$ , (b) curves of $\alpha$ under different ${\delta }_{z,max}$ , (c) curves of ${\omega }_z$ under different ${\delta }_{z,max}$ , and (d) curves of saturated ${\delta }_z$ .

Simulating with initial distance, velocity, measurement noise

With the help of Monte Carlo theory, large number of numerical simulations under a wide dispersion range of velocity, initial distance, and measurement noise are made. The results are shown in Figures 8–10. Figure 8 presents the Monte Carlo simulation results under initial distance dispersion between 18 and 22 km. The mean value is 20 km. It can be seen that the corresponding flight states converge well; the main difference lies in the convergence time under initial distance dispersion.

Figure 8.

Monte Carlo simulation under $R_0$ dispersion: (a) curves of $\stackrel{·}{q}$ under $R_0$ dispersion, (b) curves of $\alpha$ under $R_0$ dispersion, (c) curves of ${\omega }_z$ under $R_0$ dispersion, and (d) curves of saturated ${\delta }_z$ .

Figure 9.

Monte Carlo simulation under $V$ dispersion: (a) curves of $\stackrel{·}{q}$ under $V$ dispersion, (b) curves of $\alpha$ under $V$ dispersion, (c) curves of ${\omega }_z$ under $V$ dispersion, and (d) curves of saturated ${\delta }_z$ .

Figure 10.

Monte Carlo simulation under measurement noises: (a) curves of $\stackrel{·}{q}$ , (b) curves of $\alpha$ , (c) curves of ${\omega }_z$ , and (d) curves of saturated ${\delta }_z$ .

In Figure 9, the simulation results under velocity dispersion between 1000 and 1400 m/s are presented. The mean value is 1200 m/s. It is seen that the parameters’ setting of the proposed controller can satisfy large variations in flight velocity. From Figures 8 and 9, the proposed scheme is able to achieve the control object under large dispersion of flight states.

The measurement noises are modeled as a first-order inertia transfer function with a time constant ${\tau }_{\mathrm{senor}}$ and uncertain gain $K_{\mathrm{sensor}}$ in feedback

$$G_{\mathrm{sensor}}=\frac{K_{\mathrm{sensor}}}{{\tau }_{\mathrm{sensor}}s+1}$$

(47)

In Figure 10, the Monte Carlo simulation results under measurement noises are presented. It can be concluded that the uncertain gain of sensor dynamics has a very limited effect on the proposed controller.

Conclusion

In this article, a novel composite IGC scheme combined with third-order actuator dynamics under control input saturation and ESO is developed to address hypersonic missile flight control with multiple uncertainties and control constraint. The nonlinearities and unmatched time-varying disturbances are well estimated by three ESOs. Four differentiators including a third-order hybrid nonlinear differentiator are used to calculate the derivatives of pseudo control input. Thus, the peaking phenomenon and chasing of back-stepping sliding mode controller are greatly depressed. Besides, the noteworthy feature of the proposed IGC approach is that the control surface deflection is constraint under saturation with Nussbaum gain. It is very important in practical application.