# Missile IGC Based on Improved Model Predictive Control and Sliding Mode Observer

Authors: Li Ma; Jiayuan Shan; Junhui Liu; Yan Ding
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

Considering recurrent optimization process in model predictive control (MPC), the model uncertainties and disturbances terms in the missile’s guidance and control model can degrade recursive feasibility, and there are control mutation problems in common MPC algorithm. This paper presents a disturbance rejection model predictive control algorithm for missile integrated guidance and control (IGC). Firstly, a sliding mode observer (SMDO) is designed to estimate the unknown disturbances caused by target maneuvering. Secondly, the method of optimizing control increment is adopted in MPC to avoid the phenomenon of control mutation in the model calculation. By limiting the control increment in each cycle, it ensures the continuity of the control input. Thirdly, by combining the SMDO and MPC, an IGC algorithm is presented, and the stability of the algorithm is proved by using Lyapunov stability theory. Finally, the simulation results with different impact angles verify the effectiveness of the proposed algorithm for intercepting maneuver target.


---

## Body

## 1. Introduction

In the traditional missile guidance and control system design, the guidance loop and the control loop are designed separately, which leads to the time delay problem. And, at the end of engagement, the two loops couple fiercely, which leads to instability of guidance loop and causes large miss distance. At present, the difficulties and key problems in the integrated design of guidance and control can be summarized in these two aspects: (1) The integrated system of guidance and control has a higher order and more state variables in the mathematical model. And it exists a certain coupling between these state variables. (2) Mathematical model of guidance and control system has a lot of uncertainties. The aerodynamic parameters and the aerodynamic torque parameters are complex function of Mach number, atmosphere density, angle of attack, and deflection of rudder. However, there are uncertainties in the actual atmospheric density, wind speed, and target maneuver. Therefore, how to deal with the uncertainties in the system is also the key problem that requires to be considered in the integrated design of guidance and control system.In the past decades, integrated guidance and control has been extensively studied in flight vehicle guidance and control fields. References [1–3] designed the integrated control law based on optimal control method. Reference [1] designed a three-dimensional guidance law using a stable numerical state-dependent coefficient (SDC) method. And the Ricatti equation was solved via the state-dependent Ricatti equation (SDRE) method. Compared with the SDRE, Hamilton-Jacobi-Bellman (HJB) equation was solved by using the θ‐D method [2], which can achieve an approximate suboptimal solution. In Ref. [3], polynomial expansion of the force and torque coefficients on their respective parameters (angle of attack, sideslip angle, and rudder deflection angle) was utilized; the equation of state can be transformed into an affine linear equation dependent on the state. However, the main problems of these optimal integrated control laws are as follows: in order to get the solvable form of optimal control, the model needs more approximation and simplification; how to solve the optimal control when there are uncertainties; how to select the optimal index weight and so on. In Refs. [4–8], sliding mode control (SMC) and active disturbance rejection control (ADRC) were used to cope with model uncertainties and external disturbances in the guidance and control system. References [4, 5] considered multiple layers SMC. In Ref. [6], an integral sliding mode surface which can shorten response time was designed without considering the constraints of variables. The bound of the uncertain dynamics was considered in Refs. [7, 8] against stationary targets. Although the uncertainties of flight control dynamics were considered in the design of guidance law in Refs. [4–8], but in fact, the guidance and control were designed separately. References [9–19] combined backstepping control (BC) and robust control to design IGC algorithms. In Ref. [9], adaptive neural networks were trained online to compensate for unknown uncertainties were proposed. In Ref. [10], an adaptive nonlinear IGC method was proposed, which made the line-of-sight angular rate stabilized under the uncertainties of target maneuvering and missile model. Reference [11] proposed an IGC algorithm based on adaptive fuzzy sliding mode control, which is robust to system uncertainties and external disturbances. In Ref. [12], a simple adaptive block dynamic surface control algorithm was proposed, which can reduce computation cost. Combined with the extended state observer, the three-dimensional control of dynamic surface control algorithm was proposed in Refs. [13–15] to obtain more stable control input. In Ref. [16], an IGC method based on nonsingular terminal dynamic surface control was proposed to obtain better control effect without singular phenomenon. In Ref. [17], the systematic backstepping technique is adopted in IGC law design, and a new smooth second-order sliding mode differentiator is developed to estimate the derivatives of the virtual control laws. In Ref. [18], IGC law is developed by combining backstepping sliding mode approach. An adaptive fuzzy approximator is employed to approximate the system uncertain functions, and a sliding mode is used to compensate the effects of approximate error and disturbances. In Ref. [19], the sliding mode control method and nonlinear extended disturbance observer technique were combined to derive the control accuracy of the impact angle under the nonlinear disturbances. In Refs. [9–19], backstepping control method was utilized to realize integrated guidance and control, but there is still intermediate delay in this method. References [20–24] studied model predictive control-based guidance and control of missile. Only ideal model of missile dynamics was used to design model predictive control algorithm. In Refs. [20–22], the target maneuver and other disturbances terms had not been fully considered. Reference [23] proposed a neural network to optimize model predictive control, but considered the acceleration of target as a disturbance. A robust model predictive control was proposed in Ref. [24] against a stationary target. In Refs. [20–24], the conditions of target and disturbance are not considered in the actual environment.Based on the above analysis, to solve the recursive feasibility and control continuity of MPC-based IGC considering the indeterminate and disturbance terms in system, this paper proposes a new IGC algorithm by combining SMDO and MPC. Firstly, considering the uncertainties and interferences problem in recurrent optimization, a new sliding mode observer (SMDO) is introduced and compensates the uncertainties and disturbances in the system to solve the recursive problem in MPC-based IGC. Secondly, this paper adopts the method that optimizes control increment in MPC instead of control input and limits the control increment in each cycle to ensure the continuity of the control command. In brief, based on the new continuous MPC and SMDO composite algorithm, the integrated missile guidance and control algorithm is proposed, and this method is verified under various simulation cases. Simulation results demonstrate that the proposed algorithm is superior to the standard model predictive control. The main contributions of this paper are as follows:(1)
A new sliding mode observer (SMDO) is designed to improve the anti-interference capability and algorithm’s stability without measuring the target’s motion information and then solves the recursive problem in optimization(2)
An improved model predictive control algorithm (MPC) that optimizes the control increments is proposed instead of the standard model predictive control that optimizes the control command. The proposed MPC algorithm limits the control increment in each cycle and ensures the continuity of the control command, thus avoiding the control mutation in system(3)
Combining the model predictive control and sliding mode observer, a new composite MPC algorithm is proposed. The stability of proposed algorithm is proved. And the simulation results show that the proposed MPC method improves the performance of the closed system compared with the standard MPCThe rest of this paper is orgnized as follows: in Section2, the model of IGC is established; in Section 3, the model predictive control and SMDO-based IGC method is proposed; in Section 4, a new disturbance observer is designed; in Section 5, the stability of the method is proved. Section 6 is the results of simulation. Finally, Section 7 is the conclusion.


## 2. Integrated Dynamic Model

### 2.1. Nonlinear Kinematics Model

Suppose that the missile and the target move in the pitch plane; then, the relative motion of missile and target is shown in Figure1, where M and T represent the missile and the target, respectively. R denotes the relative distance between missile and target, q denotes the line-of-sight angle between missile and target (LOS angle). Suppose that the missile and target have constant velocities VM and VT, and their accelerations are aM and aT, respectively. The flight path angles during flight are θM and θT.Figure 1
The relative motion of the missile and target.The relative motion of the missile and the target is modeled as follows:(1)Rq•=VMsinq−θM−VTsinq−θT,R•=−VMcosq−θM+VTcosq−θT.Differentiate Equation (1), and it yields the following:
(2)Rq••+R•q•=−VT•sinq−θT+VM•sinq−θM+VTθ•Tcosq−θT−VMθ•Mcosq−θM−q•VTcosq−θT+q•VMcosq−θM.Since theθ•M=aM/VM,θ•T=aT/VT and according to the assumption that the missile and the target have constant velocities, then we have V•M=0,V•T=0. Equation (2) can be simplified as follows:
(3)q••=−2R•q•+aTcosq−θT−aMcosq−θMR.


### 2.2. Missile Dynamics Model

The missile dynamics in the longitude plane is modeled as follows:(4)α•=ωz+1mVMmgcosθM−Y,ω•z=1JzM0+Mδzδz,υ•=ωz,θM=υ−α,where α is the attack angle, m is the mass of missile, g is gravity constant, Jz is the moment of inertia about the Z axis, υ is the pitch angle, ωz is the pitching rate, δz is the deflection angle of elevator for pitch control, Y is the lift force, Mδz is the control contribution to the angular acceleration, M0=M0α,ωz,VM,mzα,mzωz represents the angular acceleration contributions from all other sources, the mzα denotes pitch moment derivative with respect to α, and mzωz represents the pitch moment derivative with respect to ωz.The lift and moment terms are defined as follows:(5)Y=57.3QScyαα,M0=57.3QSLmzαα+QSL2mzωzωzVM,Mδz=57.3QSLmzδz,where Q is the dynamic pressure, which equals to 1/2ρVM2, Sis the missile reference area, L is the reference length, cyα is the lift force derivative with respect to α, and mzδz represents the pitch moment derivatives with respect to δz.According to Equation (4), we have θ•M=υ•−α•⟶aM=VMθ•M=VMY−mgcosθM/mVM; then, Equation (3) can be rewritten as follows:
(6)q••=1R−2R•q•+acTosq−θT−aMcosq−θM=−2R•Rq•−Y−mgcosθMcosq−θMmR+acTosq−θTR=−57.3QScyαcosq−θMmRα−2R•Rq•+gcosθMcosq−θM+acTosq−θTR.Definex1,x2,x3 as q•,α,ωz, respectively, and set u=δz. From (4)–(6), the integrated equation of guidance and control can be derived as
(7)x•1=a1x1+a2x2+d1,x•2=x3+k2x2+d2,x•3=k3x2+k4x3+bu,where d1,d2 and a1,a2,k2,k3,k4,b are defined as follows:
(8)d1=gcosθMcosq−θM+acTosq−θTR,d2=gcosθMVM+d¯2Δ,a1=−2R•R,a2=−57.3QScyαcosq−θMmR,k2=−57.3QScyαmVM,k3=57.3QSLmzαJz,k4=QSL2mzωzJzVM,b=57.3QSLmzδzJz.The disturbanced1 contains the unknown target acceleration and the time-varying term relative to system states, while the d2 represents the time-varying disturbance caused by the uncertainties of aerodynamic parameters and gravity, where d¯2Δ denotes the disturbance due to aerodynamic uncertainties.Based on the control problem formulated above, we will design a composite MPC scheme that consists of the nominal input and the estimated disturbances, i.e.,ut=κu~t,d⌢t to offset the effect caused by the disturbances, where d^t is the estimation of dt and u~t is obtained by minimizing the cost function of MPC. The idea is shown in Figure 2.Definition 1.
The impact angle is expressed asi, which means the angle between the velocity vector of missile and the velocity vector of target at impact time, so i=θT−θM. When the interception is realized, it can be obtained that q•=0, and the following equation can be obtained from Equation (1).
(9)VMsinq−θM=VTsinq−θT.Figure 2
Block diagram of composite MPC.According to Equation (9), there is always a unique expected terminal line-of-sight angle for a predetermined target with an expected impact angle. Therefore, the control problem of the impact angle can be transformed into the control of the final line-of-sight angle.SubstitutingθM=θT−i into Equation (9), then it can be given that
(10)VMsinq−θT+i=VTsinq−θT,(11)VMsinq−θTcosi+cosq−θTsini=VTsinq−θT.From (10) and (11), it can be obtained that
(12)q=θT−arctansinicosi−VT/VM.Finally, the control objective is to design a control law that enables the missile to impact the maneuvering target with the required impact anglei.


## 2.1. Nonlinear Kinematics Model

Suppose that the missile and the target move in the pitch plane; then, the relative motion of missile and target is shown in Figure1, where M and T represent the missile and the target, respectively. R denotes the relative distance between missile and target, q denotes the line-of-sight angle between missile and target (LOS angle). Suppose that the missile and target have constant velocities VM and VT, and their accelerations are aM and aT, respectively. The flight path angles during flight are θM and θT.Figure 1
The relative motion of the missile and target.The relative motion of the missile and the target is modeled as follows:(1)Rq•=VMsinq−θM−VTsinq−θT,R•=−VMcosq−θM+VTcosq−θT.Differentiate Equation (1), and it yields the following:
(2)Rq••+R•q•=−VT•sinq−θT+VM•sinq−θM+VTθ•Tcosq−θT−VMθ•Mcosq−θM−q•VTcosq−θT+q•VMcosq−θM.Since theθ•M=aM/VM,θ•T=aT/VT and according to the assumption that the missile and the target have constant velocities, then we have V•M=0,V•T=0. Equation (2) can be simplified as follows:
(3)q••=−2R•q•+aTcosq−θT−aMcosq−θMR.


## 2.2. Missile Dynamics Model

The missile dynamics in the longitude plane is modeled as follows:(4)α•=ωz+1mVMmgcosθM−Y,ω•z=1JzM0+Mδzδz,υ•=ωz,θM=υ−α,where α is the attack angle, m is the mass of missile, g is gravity constant, Jz is the moment of inertia about the Z axis, υ is the pitch angle, ωz is the pitching rate, δz is the deflection angle of elevator for pitch control, Y is the lift force, Mδz is the control contribution to the angular acceleration, M0=M0α,ωz,VM,mzα,mzωz represents the angular acceleration contributions from all other sources, the mzα denotes pitch moment derivative with respect to α, and mzωz represents the pitch moment derivative with respect to ωz.The lift and moment terms are defined as follows:(5)Y=57.3QScyαα,M0=57.3QSLmzαα+QSL2mzωzωzVM,Mδz=57.3QSLmzδz,where Q is the dynamic pressure, which equals to 1/2ρVM2, Sis the missile reference area, L is the reference length, cyα is the lift force derivative with respect to α, and mzδz represents the pitch moment derivatives with respect to δz.According to Equation (4), we have θ•M=υ•−α•⟶aM=VMθ•M=VMY−mgcosθM/mVM; then, Equation (3) can be rewritten as follows:
(6)q••=1R−2R•q•+acTosq−θT−aMcosq−θM=−2R•Rq•−Y−mgcosθMcosq−θMmR+acTosq−θTR=−57.3QScyαcosq−θMmRα−2R•Rq•+gcosθMcosq−θM+acTosq−θTR.Definex1,x2,x3 as q•,α,ωz, respectively, and set u=δz. From (4)–(6), the integrated equation of guidance and control can be derived as
(7)x•1=a1x1+a2x2+d1,x•2=x3+k2x2+d2,x•3=k3x2+k4x3+bu,where d1,d2 and a1,a2,k2,k3,k4,b are defined as follows:
(8)d1=gcosθMcosq−θM+acTosq−θTR,d2=gcosθMVM+d¯2Δ,a1=−2R•R,a2=−57.3QScyαcosq−θMmR,k2=−57.3QScyαmVM,k3=57.3QSLmzαJz,k4=QSL2mzωzJzVM,b=57.3QSLmzδzJz.The disturbanced1 contains the unknown target acceleration and the time-varying term relative to system states, while the d2 represents the time-varying disturbance caused by the uncertainties of aerodynamic parameters and gravity, where d¯2Δ denotes the disturbance due to aerodynamic uncertainties.Based on the control problem formulated above, we will design a composite MPC scheme that consists of the nominal input and the estimated disturbances, i.e.,ut=κu~t,d⌢t to offset the effect caused by the disturbances, where d^t is the estimation of dt and u~t is obtained by minimizing the cost function of MPC. The idea is shown in Figure 2.Definition 1.
The impact angle is expressed asi, which means the angle between the velocity vector of missile and the velocity vector of target at impact time, so i=θT−θM. When the interception is realized, it can be obtained that q•=0, and the following equation can be obtained from Equation (1).
(9)VMsinq−θM=VTsinq−θT.Figure 2
Block diagram of composite MPC.According to Equation (9), there is always a unique expected terminal line-of-sight angle for a predetermined target with an expected impact angle. Therefore, the control problem of the impact angle can be transformed into the control of the final line-of-sight angle.SubstitutingθM=θT−i into Equation (9), then it can be given that
(10)VMsinq−θT+i=VTsinq−θT,(11)VMsinq−θTcosi+cosq−θTsini=VTsinq−θT.From (10) and (11), it can be obtained that
(12)q=θT−arctansinicosi−VT/VM.Finally, the control objective is to design a control law that enables the missile to impact the maneuvering target with the required impact anglei.


## 3. IGC Design Based on Improved Model Predictive Control

### 3.1. Design SMDO for Unknown Disturbances

In this section, we consider the case that the disturbances are unknown. Compared with the first-order sliding mode observer, the second-order sliding mode observer (SOSM) [25] or higher-order sliding mode observer (HOSM) [26] significantly reduces the undesired chattering phenomenon, which is very valuable in practical applications. Based on the previous analysis, in this section, we propose a new SMDO based on SOSM and adaptive control method to estimate target maneuver.An enhanced nonlinear disturbance observer is proposed to estimate the disturbance of the IGC system. The nonlinear observer is designed as follows:(13)d^1=z1+λ1x1,z•1=−λ1z1−λ12x1−L1M1x1+M2x2,where d^1 denotes the estimated value of disturbance d1, z1∈Rl is the internal state of the nonlinear observer, λ1 and L1 are the gain coefficient of the nonlinear observer, and M1=−R•/R,M2=−57.3QScyαcosq−θM/mR.Set the error as follows:(14)e1=d1−d^1.Similarly, the estimated valued^2 is designed as follows:
(15)d^2=z2+λ2x2,z•2=−λ2z2−λ22x2−L2x3+k2x2,where d^2 denotes the estimated value of disturbance d2, z2∈Rl is the internal state of the nonlinear observer, and λ2 and L2 are the gain coefficient of the nonlinear observer, and set the error as follows:
(16)e2=d2−d^2.Design the sliding mode surface as follows:(17)S=c1e1+c2e2.And, the reaching law of sliding mode surface is designed as follows:(18)S•=−ks−εsgnsk>0,ε>0.In order to prove the convergence of the estimated error, firstly define a Lyapunov function:(19)V=12S2.Then, differentiate (19), and it can be obtained that
(20)V•=s⋅s•=s−ks−εsgns=−ks2−εs=−sks+ε.Due tok>0,ε>0, it can be given that
(21)ks+ε≥0,V•=−sks+ε≤0.SinceV•≤0, the sliding mode surface is stable based on Lyapunov’s theorem. And then, according to the finite time convergence theory [27], it can be concluded that the estimated error converges to zero in finite time. Finally, the nonlinear observers (13) and (15) are also finite-time stable, which means that the estimated error converges to zero in finite time.Remark 2.
The gain coefficients of the nonlinear observer should be positive constants. And, it needs to be selected carefully to limit the estimated error.


### 3.2. Linearization of Systems

Since both kinematics model and dynamics model are nonlinear, they should be linearized to utilize linear model predictive control method. Compared to nonlinear model predictive control, linear model predictive control has better real-time performance and is easier to analyze and calculate. The standard MPC algorithms in [22, 23] optimize control command directly. But in this section, the method that optimizes control increment is adopted.According to the integrated model, the state variables and control variables are as follows:(22)x=q•αωz,u=00δz.The corresponding Jacobian matrix is as follows:(23)A=a1a200k210k3k4,B=00b.The linear time-varying model is derived as follows:(24)x¯•=A¯x¯+B¯u¯,where
(25)x¯=q•−q•0α−α0ωz−ωz0,u¯=00δz−δz0,A¯=I+TA=1−Ta1Ta2001+Tk2T0Tk31+Tk4,B¯=TB=00Tb,where q•0,α0,ωz0,δz0 are the initial values of q•,α,ωz,δz.Then, the discrete linear model is obtained as follows:(26)x¯k+1=A¯x¯k+B¯u¯k,where k denotes the current step and k+1 denotes the next step. And the output of the system is as follows:
(27)yk=Cx¯k.In order to solve control mutation problem of MPC, the method that optimizes control increment in MPC instead of control command is adopted. It is assumed that(28)ηk∣t=x¯k∣tu¯k−1∣t,where t is the time. Then, a new state equation can be obtained as follows:
(29)ηk+1∣t=A~ηk∣t+B~Δuk∣t,ξk∣t=C~ηk∣t,where
(30)A~=A¯B¯0I,B~=B¯I,C~=C0.Through continuous iteration, it can be given that(31)ηk+2∣t=A~ηk+1∣t+B~Δuk+1∣t=A~2ηk+1∣t+A~B~Δut+B~Δuk+1∣t,ηk+N∣t=A~ηk+N−1∣t+B~Δuk+N−1∣t=A~Nηk+N∣t+A~N−1B~Δut+⋯+B~Δuk+N−1∣t,where N is the predictive time domain, and the larger the N is, the better the predictive effect produces, but more calculation cost will take. Finally, the state and output of the system in the future time are expressed as follows:
(32)ηk+1∣tηk+2∣t⋯ηk+N∣t=A~A~2⋯A~Nηt∣t+B~0⋯0A~B~B~⋯0⋯⋯⋯⋯A~N−1B~A~N−2B~B~Δuk∣tΔuk+1∣t⋯Δuk+N−1∣t.Then, it can be simplified as follows:(33)Xt=Atηt∣t+BtΔUt∣t,where
(34)Xt=ηk+1∣tηk+2∣t⋯ηk+N∣t,ΔUt∣t=Δuk∣tΔuk+1∣t⋯Δuk+N−1∣t,At=A~A~2⋯A~N,Bt=B~0⋯0A~B~B~⋯0⋯⋯⋯⋯A~N−1B~A~N−2B~B~.The state and output in the predictive time domain can be obtained from the current state quantityηt∣t and the control increment ΔUt∣t of the system, which is the predictive function in the model predictive control algorithm. The schematic is presented in Figure 3.Figure 3
MPC schematic.


## 3.1. Design SMDO for Unknown Disturbances

In this section, we consider the case that the disturbances are unknown. Compared with the first-order sliding mode observer, the second-order sliding mode observer (SOSM) [25] or higher-order sliding mode observer (HOSM) [26] significantly reduces the undesired chattering phenomenon, which is very valuable in practical applications. Based on the previous analysis, in this section, we propose a new SMDO based on SOSM and adaptive control method to estimate target maneuver.An enhanced nonlinear disturbance observer is proposed to estimate the disturbance of the IGC system. The nonlinear observer is designed as follows:(13)d^1=z1+λ1x1,z•1=−λ1z1−λ12x1−L1M1x1+M2x2,where d^1 denotes the estimated value of disturbance d1, z1∈Rl is the internal state of the nonlinear observer, λ1 and L1 are the gain coefficient of the nonlinear observer, and M1=−R•/R,M2=−57.3QScyαcosq−θM/mR.Set the error as follows:(14)e1=d1−d^1.Similarly, the estimated valued^2 is designed as follows:
(15)d^2=z2+λ2x2,z•2=−λ2z2−λ22x2−L2x3+k2x2,where d^2 denotes the estimated value of disturbance d2, z2∈Rl is the internal state of the nonlinear observer, and λ2 and L2 are the gain coefficient of the nonlinear observer, and set the error as follows:
(16)e2=d2−d^2.Design the sliding mode surface as follows:(17)S=c1e1+c2e2.And, the reaching law of sliding mode surface is designed as follows:(18)S•=−ks−εsgnsk>0,ε>0.In order to prove the convergence of the estimated error, firstly define a Lyapunov function:(19)V=12S2.Then, differentiate (19), and it can be obtained that
(20)V•=s⋅s•=s−ks−εsgns=−ks2−εs=−sks+ε.Due tok>0,ε>0, it can be given that
(21)ks+ε≥0,V•=−sks+ε≤0.SinceV•≤0, the sliding mode surface is stable based on Lyapunov’s theorem. And then, according to the finite time convergence theory [27], it can be concluded that the estimated error converges to zero in finite time. Finally, the nonlinear observers (13) and (15) are also finite-time stable, which means that the estimated error converges to zero in finite time.Remark 2.
The gain coefficients of the nonlinear observer should be positive constants. And, it needs to be selected carefully to limit the estimated error.


## 3.2. Linearization of Systems

Since both kinematics model and dynamics model are nonlinear, they should be linearized to utilize linear model predictive control method. Compared to nonlinear model predictive control, linear model predictive control has better real-time performance and is easier to analyze and calculate. The standard MPC algorithms in [22, 23] optimize control command directly. But in this section, the method that optimizes control increment is adopted.According to the integrated model, the state variables and control variables are as follows:(22)x=q•αωz,u=00δz.The corresponding Jacobian matrix is as follows:(23)A=a1a200k210k3k4,B=00b.The linear time-varying model is derived as follows:(24)x¯•=A¯x¯+B¯u¯,where
(25)x¯=q•−q•0α−α0ωz−ωz0,u¯=00δz−δz0,A¯=I+TA=1−Ta1Ta2001+Tk2T0Tk31+Tk4,B¯=TB=00Tb,where q•0,α0,ωz0,δz0 are the initial values of q•,α,ωz,δz.Then, the discrete linear model is obtained as follows:(26)x¯k+1=A¯x¯k+B¯u¯k,where k denotes the current step and k+1 denotes the next step. And the output of the system is as follows:
(27)yk=Cx¯k.In order to solve control mutation problem of MPC, the method that optimizes control increment in MPC instead of control command is adopted. It is assumed that(28)ηk∣t=x¯k∣tu¯k−1∣t,where t is the time. Then, a new state equation can be obtained as follows:
(29)ηk+1∣t=A~ηk∣t+B~Δuk∣t,ξk∣t=C~ηk∣t,where
(30)A~=A¯B¯0I,B~=B¯I,C~=C0.Through continuous iteration, it can be given that(31)ηk+2∣t=A~ηk+1∣t+B~Δuk+1∣t=A~2ηk+1∣t+A~B~Δut+B~Δuk+1∣t,ηk+N∣t=A~ηk+N−1∣t+B~Δuk+N−1∣t=A~Nηk+N∣t+A~N−1B~Δut+⋯+B~Δuk+N−1∣t,where N is the predictive time domain, and the larger the N is, the better the predictive effect produces, but more calculation cost will take. Finally, the state and output of the system in the future time are expressed as follows:
(32)ηk+1∣tηk+2∣t⋯ηk+N∣t=A~A~2⋯A~Nηt∣t+B~0⋯0A~B~B~⋯0⋯⋯⋯⋯A~N−1B~A~N−2B~B~Δuk∣tΔuk+1∣t⋯Δuk+N−1∣t.Then, it can be simplified as follows:(33)Xt=Atηt∣t+BtΔUt∣t,where
(34)Xt=ηk+1∣tηk+2∣t⋯ηk+N∣t,ΔUt∣t=Δuk∣tΔuk+1∣t⋯Δuk+N−1∣t,At=A~A~2⋯A~N,Bt=B~0⋯0A~B~B~⋯0⋯⋯⋯⋯A~N−1B~A~N−2B~B~.The state and output in the predictive time domain can be obtained from the current state quantityηt∣t and the control increment ΔUt∣t of the system, which is the predictive function in the model predictive control algorithm. The schematic is presented in Figure 3.Figure 3
MPC schematic.


## 4. MPC with Disturbance Compensation

The purpose of receding horizon optimization is to find the optimal control solution, which is an online optimization. It recalculates the current error at each time and finally achieves desired control performance by optimizing one or some cost functions. Therefore, receding horizon optimization may not obtain the global optimal solution, but it can respond to the state at each time to achieve local optimal. The idea of MPC is presented as Figure4.Figure 4
Receding horizon optimization.Therefore, it is necessary to design appropriate cost function to obtain the best possible control sequence. The cost function can be written as a quadratic function of state and control input.(35)J=∑i=1NXtTQ~Xt+ΔUTR~ΔU,where Q~ and R~ are the weight matrices, and they are positive definite matrices which satisfy Q~≥0,R~≥0. The first term in (35) is to tuning the tracking performance of the reference trajectory, and the second term is the constraint on control consumption.In MPC method, the missile needs to follow the reference trajectory and satisfies the constraints of states. Therefore, the following aspects should be considered when designing the cost function: (1) maximum limits on angle of attack, (2) maximum limits on pitching rate, (3) maximum limits on elevator deflection angle and the increment, and (4) terminal state of line-of-sight angular rate.(36)minΔutJ,(37)αmint+k≤αt+k≤αmaxt+k,(38)ωmint+k≤ωt+k≤ωmaxt+k,(39)umint+k≤ut+k≤umaxt+k,(40)Δumint+k≤Δut+k≤Δumaxt+k,(41)q•t+k=0.Remark 3.
Considering practical application limitations, such as strapdown seeker’s field of view, maneuverability of missile, and limits of actuator’s response, appropriate constraints and parameters should be selected carefully.After receding horizon optimization, control increment sequence is obtained in each control cycle.(42)ΔUt∣t=Δut∣t,Δut+1∣t,⋯,Δut+N−1∣t.The first control increment of the sequence is applied to the system as the actual control input increment.(43)u~t=ut−1+Δut∣t,where u~t is the control input to the system at current time. At the next control time, the system predicts the output in the next period according to the new state information and then obtains a new set of control increments through optimization. Next, repeat until the whole control process is completed. The MPC algorithm is summarized in the following Algorithm 1.Algorithm 1: MPC algorithm.
Algorithm: MPC methodology1: Obtain the initial variables and state matrices of the system through the initial state of the missile.2: Discrete the original system to obtain (26).3: Take the iteration processing to obtain (33), and the N is the number of iterations.4: The optimization objective function (35) is determined by (33), and then, the quadratic programming is performed under constraints (37)–(41).5: Finally, a set of control increments (42) is obtained, and the Δut∣t is imposed on the original system to obtain the minimizing control sequence u~t. And the estimate of the disturbance d^t is solved based on SMDO.6: Apply the composite control signalut=κu~t,d⌢t to the real system.7: Update the time oft⟶t+1 again and return to the first step to restart.Remark 4.
It can be seen from the above algorithm that the states of the system are updated by the actual states of each step, and this problem needs to be solved online. This strategy obtains an optimal control about the current state. Therefore, this method has the same advantages as the online optimization method.


## 5. Stability Proof

Assumption. Define F⋅ as a terminal cost, Xf denotes the terminal constraint set, and ψf⋅ denotes the time-varying control set. Suppose that F⋅,Xf,ψf⋅ satisfy the following conditions:A1: Xf⊂R. Xf is closed-loop, and 0∈Xf,∀x∈Xf (x is the state, and state constraint satisfied in Xf).A2: ψfx∈R,∀u∈ψfx (u is the control, and control constraint satisfied in ψf⋅).A3: fx,u,d∈Xf,∀x∈Xf,∀d∈Dx,u (d is the disturbance, D⋅ is a disturbance sequence, and f⋅ is a function about x,u,d).A4: F∗+lx,ψfx,d≤0,∀x∈Xf,∀d∈Dx,ψfx (F∗x,u represents the change due to increment of x, and l⋅ denotes a phase cost function).IfF is a Lyapunov function with robust control, based on these assumptions, one can obtain the following:
(44)VN0∗+lx,ψNx,d≤0.The above equation holds for all appropriate sets of x and all d∈Dx,ψNx. Therefore, according to [28], the above conditions can finally ensure the stability of the system.Notation. For any function F⋅, F∗x,u represents the change of state from x to x+=fx,u, and it can be defined as follows:
(45)F∗x,u=Ffx,u−Fx.Proof.
Since the system in this paper has complex characteristics such as nonlinear, time-varying, and constrained, the system needs to be discretized firstly and then prove the stability.
First, determine a phase cost functionl⋅, and then, the discrete system is obtained according to Equation (22).
(46)xk+1=fxk,uk,dk,(47)yk=hxk,where f⋅ is a function about xk. Then, the above equation can be rewritten as follows:
(48)x+=fx,u,d,(49)y=Hx,u,
When the policyπ=u is used, the following specific cost function form is defined as follows:
(50)Jx,π,d=∑i=0N−1lzi,di+FxN,where xi=xπ,di;x,0 and ly,d=1/2y2−γ2d2. According to (50), the function VNx,π is defined as
(51)VNx,π=maxJx,π,d∣d∈DNx,π,where DNx,u is a set of constrained values for disturbances. The optimal control is now defined as follows:
(52)ΦNx:VN0x=minVNx,π∣π∈ΓNx,where ΓNx is a set that satisfies the conditions of optimal control.
To simplify the stability argument,Dx,u is now assumed.
(53)Dx,u=d∣d≤δHx,u,where setting the Q=HTH>0, δ<1/γ, and it can be obtained that
(54)ly,d=12y2−γ2d21−δ2γ2y2≥12=σ2y2≥u2x2,where σ>0, μ>0, and the equation satisfies all x, u, andd∈Dx,u.
Next, in order to determineXf and F⋅, Equation (22) is rewritten as follows:
(55)x+=Ax+Bu+Gd,where A=fx0,0, B=fu0,0, G=fd0,0, and its realization of the stage cost is follows:
(56)∑i=0∞l′xi,ui,di.
Define the stage cost function as follows:(57)l′x,u,d=lHx,d,d+12ρx2,ρ>1.
When the data is appropriate, especially whenγ is large enough, and utilizing the optimal control strategy u=ψfx, then the following Lyapunov function can be selected as the cost function:
(58)Fx=12xTPx,P>0.
Set the cost functionF⋅ to satisfy the following:
(59)FAx+Bψfx+Gd−Fx+lx,ψfx,d+12ρx2≤0,where all the x and d satisfy the above equation. If the function f⋅ is smooth enough, then the F⋅ has the following set:
(60)Xf≔x∣Fx≤c.
And satisfying the following inequalities,(61)F∗+lx,ψfx,d≤0,where all the x∈Xf,d∈Dx,u are satisfied. From the above assumption, the F⋅,Xf,ψf⋅ satisfy the A1-A4 conditions, so it can be obtained that
(62)VN0∗+lx,ψNx,d≤0,where it satisfies all the x∈XN,x+∈Fx,ψNx,d∈Dx,u. Therefore, it is finally proved that the system is stable.☐Remark 5.
It should be emphasized that the stability proof conditions of MPC are limited, which means that the state variables should change within a certain range according to different cases. Therefore, the system is stable only when the set conditions are satisfied. Due to the particularity of the MPC model and the difference of state matrix at each time, this paper only needs to prove the stability of the initial time. And the proof method at the later time is the same, which only needs to change within the boundary.


## 6. Numerical Simulation

### 6.1. Simulation Parameters

The initial conditions of missile and target are as follows: the initial positions are set asxA0yA0=0m0m andxB0yB0=5000/3m5000m; the velocities are VA0=500m/s and VB0=250m/s; the initial flight path angles are set as θA0=30∘ and θB0=120∘; the initial pitch angle is set as υA0=30∘; the initial pitching rate is set as ωAz0=0; the initial distance and line of sight angle between missile and target are R0=10000m and q0=30∘, respectively. Coefficients of disturbance observer are selected as λ1=10, L1=20, λ2=10, and L2=30. The nominal aerodynamic parameters which are used to generate lift L and pitch torque M are set as follows:
(63)k2=−57.3QscyαmVM=−0.3487,k3=57.3QSLmzαJz=−17.801,k4=QSL2mzωzJzVM=−0.2741,b=57.3QSLmzδzJz=−31.267.


### 6.2. Simulation Results

For comparison, two other control algorithms are introduced in the simulation. The first one is the standard model predictive control in literature [20]. The standard MPC usually optimizes control input, not the control increment, which means the control was obtained from the receding horizon optimization of Equation (18). The second one is the reference value, which comes from an integral sliding mode control algorithm. The method is introduced briefly as follows:According to literature [6], an integral sliding mode surface with all state variables is designed as follows:
(64)S=c1x1+c2x2+x3,where c1>0,c2>0, and they are all constants. It can be obtained from Equations (7) and (64).
(65)S•=c1x2+d1+c2x3+k2x2+d2+k3x2+k4x3+bu.According to the design method of reaching law, it can be designed as follows:(66)S•=−εsgnS−kS+c1d~1+c2d~2.According to Equations (65) and (66), the control law is obtained.
(67)u=−1bεsgnS+ks+c1x2+d^1+c2x3+k2x2+d^2+k3x2+k4x3,where the coefficients are c1=12, c2=9, k=3, and ε=1.Three cases are presented to demonstrate the superiority of the proposed IGC algorithm as follows:Case 1.
Target interception with unknown target’s maneuvering.
In this case, assume that the target is maneuvering with accelerationaT=62.5sin0.25tm/s2. The LOS angular rate obtained by these two algorithms is given in Figure 5(b). It can be observed that under the standard MPC law, the LOS angular rate diverges at the end of the terminal guidance phase, which results in large LOS angle tracking errors. By contrast, due to introducing SMDO, it achieves high control accuracy under the same conditions. Meanwhile, it can also be seen from Figure 5(b) that the response trend of the two algorithms is closing when t<5 s, which means that the proposed MPC algorithm preserves the performance of the standard MPC algorithm. Figures 5(c) and 5(d) also indicate that the attack angle and pitch angle of the SMDO-MPC algorithm diverge slower than standard MPC. Furthermore, Figures 5(e) and 5(f) indicate that the proposed method needs smaller control input and has satisfactory disturbance estimation performance.Case 2.
Comparison of proposed algorithm with different predictive steps.
To evaluate the influence of predictive level on algorithm’s performance, three different predictive step parameters are used in the simulation, that is,Tsim=20s, 35 s, and 50s. In this case, the parameters of SMDO are designed the same with Case 1. The LOS angle, the attack angle, and the pitch angle are given in Figures 6(b)–6(d). From these figures, we can see that those with larger predictive steps have faster convergence speed and more accurate tracking performance. But it may lead to larger guidance and control errors. So it needs to select the appropriate predictive steps. It can be concluded that the predictive step level can significantly affect the performance of integrated guidance and control. Therefore, under certain conditions, this parameter can be designed according to the performance requirement.Case 3.
Validation of SMDO-MPC under different impact angle constraints.
Assume that the target is maneuvering with accelerationaT=20sin0.25tm/s2, and the impact angles i=30∘ and 60∘. Under the proposed IGC law, the simulation results are shown as presented in Figures 7(a) and 7(b). The results indicate that the proposed IGC law has satisfactory performance, and both the LOS angle and trajectory control are within an acceptable range. It means that the proposed IGC law is within a reasonable range of impact angles, and the impact angle constraint is achieved. Figure 7(e) indicates that there is an acceptable field-of-view angle in simulation. Therefore, the proposed IGC law can intercept maneuvering targets with impact angle constraints.Figure 5
Simulation results under Case1: (a) intercept trajectory, (b) the LOS angular rate, (c) the attack angle, (d) the pitch angle, (e) the deflection angle, and (f) real value and estimate value of disturbance.
(a)(b)(c)(d)(e)(f)Figure 6
Simulation results under Case2: (a) the LOS angular rate, (b) the LOS angle, (c) the attack angle, and (d) the pitch angle.
(a)(b)(c)(d)Figure 7
Simulation results under Case3: (a) intercept trajectory, (b) the LOS angle, (c) the attack angle, (d) the pitch angle, and (e) the field-of-view angle.
(a)(b)(c)(d)(e)


## 6.1. Simulation Parameters

The initial conditions of missile and target are as follows: the initial positions are set asxA0yA0=0m0m andxB0yB0=5000/3m5000m; the velocities are VA0=500m/s and VB0=250m/s; the initial flight path angles are set as θA0=30∘ and θB0=120∘; the initial pitch angle is set as υA0=30∘; the initial pitching rate is set as ωAz0=0; the initial distance and line of sight angle between missile and target are R0=10000m and q0=30∘, respectively. Coefficients of disturbance observer are selected as λ1=10, L1=20, λ2=10, and L2=30. The nominal aerodynamic parameters which are used to generate lift L and pitch torque M are set as follows:
(63)k2=−57.3QscyαmVM=−0.3487,k3=57.3QSLmzαJz=−17.801,k4=QSL2mzωzJzVM=−0.2741,b=57.3QSLmzδzJz=−31.267.


## 6.2. Simulation Results

For comparison, two other control algorithms are introduced in the simulation. The first one is the standard model predictive control in literature [20]. The standard MPC usually optimizes control input, not the control increment, which means the control was obtained from the receding horizon optimization of Equation (18). The second one is the reference value, which comes from an integral sliding mode control algorithm. The method is introduced briefly as follows:According to literature [6], an integral sliding mode surface with all state variables is designed as follows:
(64)S=c1x1+c2x2+x3,where c1>0,c2>0, and they are all constants. It can be obtained from Equations (7) and (64).
(65)S•=c1x2+d1+c2x3+k2x2+d2+k3x2+k4x3+bu.According to the design method of reaching law, it can be designed as follows:(66)S•=−εsgnS−kS+c1d~1+c2d~2.According to Equations (65) and (66), the control law is obtained.
(67)u=−1bεsgnS+ks+c1x2+d^1+c2x3+k2x2+d^2+k3x2+k4x3,where the coefficients are c1=12, c2=9, k=3, and ε=1.Three cases are presented to demonstrate the superiority of the proposed IGC algorithm as follows:Case 1.
Target interception with unknown target’s maneuvering.
In this case, assume that the target is maneuvering with accelerationaT=62.5sin0.25tm/s2. The LOS angular rate obtained by these two algorithms is given in Figure 5(b). It can be observed that under the standard MPC law, the LOS angular rate diverges at the end of the terminal guidance phase, which results in large LOS angle tracking errors. By contrast, due to introducing SMDO, it achieves high control accuracy under the same conditions. Meanwhile, it can also be seen from Figure 5(b) that the response trend of the two algorithms is closing when t<5 s, which means that the proposed MPC algorithm preserves the performance of the standard MPC algorithm. Figures 5(c) and 5(d) also indicate that the attack angle and pitch angle of the SMDO-MPC algorithm diverge slower than standard MPC. Furthermore, Figures 5(e) and 5(f) indicate that the proposed method needs smaller control input and has satisfactory disturbance estimation performance.Case 2.
Comparison of proposed algorithm with different predictive steps.
To evaluate the influence of predictive level on algorithm’s performance, three different predictive step parameters are used in the simulation, that is,Tsim=20s, 35 s, and 50s. In this case, the parameters of SMDO are designed the same with Case 1. The LOS angle, the attack angle, and the pitch angle are given in Figures 6(b)–6(d). From these figures, we can see that those with larger predictive steps have faster convergence speed and more accurate tracking performance. But it may lead to larger guidance and control errors. So it needs to select the appropriate predictive steps. It can be concluded that the predictive step level can significantly affect the performance of integrated guidance and control. Therefore, under certain conditions, this parameter can be designed according to the performance requirement.Case 3.
Validation of SMDO-MPC under different impact angle constraints.
Assume that the target is maneuvering with accelerationaT=20sin0.25tm/s2, and the impact angles i=30∘ and 60∘. Under the proposed IGC law, the simulation results are shown as presented in Figures 7(a) and 7(b). The results indicate that the proposed IGC law has satisfactory performance, and both the LOS angle and trajectory control are within an acceptable range. It means that the proposed IGC law is within a reasonable range of impact angles, and the impact angle constraint is achieved. Figure 7(e) indicates that there is an acceptable field-of-view angle in simulation. Therefore, the proposed IGC law can intercept maneuvering targets with impact angle constraints.Figure 5
Simulation results under Case1: (a) intercept trajectory, (b) the LOS angular rate, (c) the attack angle, (d) the pitch angle, (e) the deflection angle, and (f) real value and estimate value of disturbance.
(a)(b)(c)(d)(e)(f)Figure 6
Simulation results under Case2: (a) the LOS angular rate, (b) the LOS angle, (c) the attack angle, and (d) the pitch angle.
(a)(b)(c)(d)Figure 7
Simulation results under Case3: (a) intercept trajectory, (b) the LOS angle, (c) the attack angle, (d) the pitch angle, and (e) the field-of-view angle.
(a)(b)(c)(d)(e)


## 7. Conclusion

To enhance the recursive feasibility of standard MPC-based IGC, this paper proposes an improved MPC-based IGC algorithm by combining the SMDO with MPC. Firstly, an SMDO is designed to estimate disturbance caused by target maneuvering and model uncertainties. And in order to avoid the phenomenon of control mutation in the model calculation, this paper adopts the method of optimizing control increment in MPC and limits the control increment in each cycle, which ensures the continuity of the control command. Then, the new designs of SMDO and MPC are combined into a new composite controller. And the stability of the system is proved by utilizing Lyapunov stability theory. Simulation results indicate that the proposed IGC method has satisfactory strike performance, which ensures interception precision under unknown target maneuvering, and achieves the desired terminal impact angle. Compared with the standard MPC, the proposed method needs smaller rudder deflection input. In the future, the composite algorithm needs to be further analyzed in satisfying more constraints and the stability needs to be proved strictly.



---

*Source: 9680346-2021-11-16.xml*