Cubature Kalman Filters Model Predictive Static Programming Guidance Method with Impact Time and Angle Constraints Considering Modeling Errors

Abstract

This paper proposes a CKF-MPSP guidance method for hitting stationary targets with impact time and angle constraints for missiles in the presence of modeling errors. This innovative guidance scheme is composed of three parts: First, the model predictive static programming (MPSP) algorithm is used to design a nominal guidance method that simultaneously satisfies impact time and angle constraints. Second, the cubature Kalman filter (CKF) is introduced to estimate values of the influence of the inevitable modeling errors. Finally, a one-step compensation scheme is proposed to eliminate the modeling errors’ influence. The proposed method uses a real missile dynamics model, instead of a simplified one with a constant-velocity assumption, and eliminates the effects of modeling errors with the compensation scheme; thus, it is more practical. Simulations in the presence of modeling errors are conducted, and the results illustrate that the CKF-MPSP guidance method can reach the target with a high accuracy of impact time and angles, which demonstrates the high precision and strong robustness of the method.

1. Introduction

Guidance methods have always been a hot research topic in the field of missiles. The early guidance-law design considers only the minimized miss distance requirement, such as proportional navigation (PN). With the development of military science and technology, classical guidance methods no longer satisfy combat requirements [1]. There is an urge to research advanced guidance methods with multiple constraints. Impact angle and impact time constraints are essential for advanced terminal guidance methods. Impact angles are divided into path angle and azimuth angle. Attacking the target with a missile with proper impact angles may improve the destructive effect and hit weak parts of the target. The impact time is vital for attacking time-sensitive targets. Moreover, the guidance method with impact angle and time constraints gives the multi-missile cooperative guidance capability. Thus, investigating the impact time and angle-constrained guidance method (ITACG) is very important.

In this paper, a CKF-MPSP guidance method with impact time and angle constraints for a stationary target is proposed considering modeling errors. A baseline guidance method with impact time and angle constraints is designed based on the MPSP algorithm. The modeling errors are estimated by the CKF and compensate the baseline guidance method to eliminate their effects. The main contributions of this paper are shown below.

(1)

An ITACG is designed for a stationary target based on the MPSP algorithm, which can simultaneously achieve impact time and angle constraints. This guidance method considers the missile’s dynamic model instead of a constant-velocity model. Therefore, the proposed method is more suitable for practical missiles.

(2)

The proposed guidance method takes the desired time as a terminal condition for static planning. Time-to-go information is not required during the guidance process, which avoids the influence of time-to-go estimation errors on time control accuracy.

(3)

A CKF-based modeling error compensation scheme is proposed to solve the problem of the MPSP algorithm being unable to be used for error conditions. This improvement enhances the feasibility of the guidance method in practical applications since the modeling error is inevitable.

It is worth noting that, although the CKF-MPSP guidance method is proposed in the scenario of a missile attacking a stationary target on the ground, in this paper, it can also be used for position control, such as aircraft landing. Furthermore, by predicting and introducing the target’s motion model, the CKF-MPSP method can be used for moving-target interception. However, these applications are not the focus of this paper and require further research in the future.

This article is organized as follows. Section 2 provides a literature review of existing achievements in related research fields. Section 3 formulates the problem researched in this paper. Section 4 proposes the CKF-MPSP guidance method to implement time- and angle-constrained guidance. Simulation results are given in Section 5. Section 6 gives the conclusion.

2. Literature Review

Many scholars have conducted much research on ITACG problems. The mainstream methods can be divided into non-predictive guidance and predictive guidance methods.

The guidance laws that utilize the current relative motion information between missiles and targets are called non-predictive guidance. These kinds of methods only adopt the relative motion model instead of the actual model of the missile, making the design of such methods relatively simple. The existing non-predictive guidance methods basically follow two design paradigms.

The first design paradigm is to design guidance laws both in the line-of-sight (LOS) and normal LOS direction for time control and angle constraints, respectively. These kinds of methods are widely used in cooperative guidance scenarios. Zhang [2] proposed an ITACG with finite time convergence. Yu [1], Chen [3], and Lin [4] proposed fixed-time ITACGs based on the sliding mode theory. Ma [5] designed a disturbance-observer-based ITACG to enable the interception of maneuvering targets. Wang [6] proposed a decoupled three-dimensional sliding mode guidance law achieving simultaneous arrival at the target for multiple missiles with angle constraints. Jing [7] proposed a predefined-time convergence ITACG method for a multi-missile cooperative guidance scenario. Because of the design paradigm, the methods presented in Refs. [1,2,3,4,5,6,7] all need a control force both in the LOS and normal LOS direction for time control and angle constraints, respectively. However, most existing missiles are thrust-free and controlled by aerodynamic force in the terminal guidance period. In reality, missiles cannot provide the guidance command in the LOS direction, limiting the practical application of such guidance methods.

The second design paradigm is to design guidance laws only in the normal direction of LOS/velocity, which is more practical but more challenging compared to the above methods. Chen [8] simplified missile dynamics under a small heading error approximation and derived an optimal guidance law with impact time and angle constraints against a stationary target. Kim [9] introduced a polynomial guidance method considering impact time and angle constraints. Zhao [10] designed the trajectory as a function with two undetermined parameters and adjusted them to control the impact time and angle. Kang [11] derived a look-angle shaping scheme for ITACG. Hou [12,13] proposed a time-to-go estimation scheme for terminal sliding mode guidance with an impact angle constraint and further designed a nonsingular terminal sliding mode guidance law considering impact time and angle simultaneously. Chen [14] designed a two-stage guidance that satisfies time and angle constraints through a proper guidance-switching strategy. Zhang [15] proposed an ITACG by introducing an impact time feedback control term based on biased PN. Yan [16] proposed a computational geometry guidance against stationary targets, satisfying the constraints by iterating parameters of the geometry curve. Majumder [17] proposed a sliding-mode-control-based nonlinear guidance scheme for controlling both impact angle and impact time simultaneously. Liu [18] designed an adaptive sliding mode ITACG method, increasing its adaptability and robustness. Wang [19] designed a two-stage guidance method, achieving ITACG through reasonable switching between two guidance rules. In Refs. [8,9,10,11,12,13,14,15,16,17,18,19], the guidance methods rely on time-to-go or range-to-go estimation. However, the estimation is difficult because of the uncontrollable varying velocity. To simplify the design process, the missile’s velocity is assumed as constant, which may cause significant time estimation errors, especially when the trajectory is winding due to impact angle constraints. Thus, the guidance effect is not satisfactory in reality. Overall, there are difficulties in the practical application of the non-predictive ITACG methods in missiles because of the model mismatch.

Unlike non-predictive guidance, the predictive guidance methods predict the terminal states using a real model, which can avoid the model mismatch and thus can derive a better guidance performance. As one of the predictive guidance methods, the MPSP-based guidance methods have received widespread attention in recent years. The MPSP algorithm was first introduced in Ref [20]. Combining the philosophy of approximate dynamic programming and model predictive control, the MPSP algorithm obtains the terminal estimation of the output vector by integral prediction. Then, it efficiently solves the optimization problem by turning the dynamic programming problem into a static one. It has been applied in guidance problems due to its ability to deal with varying velocities. The existing MPSP-based guidance methods are mainly focused on the angle constraint only. Oza [21] designed an angle-constrained guidance method for an air-to-ground missile and verified the feasibility of MPSP guidance. Maity [22] further introduced the static Lagrange multiplier in MPSP guidance, improving the computational efficiency. Refs. [23,24,25,26,27] improved the computational efficiency of the MPSP algorithm through further improvements and designed angle-constrained guidance methods, respectively. Refs. [21,22,23,24,25,26,27] verified that the MPSP algorithm is feasible for solving guidance problems online. However, during integral prediction, modeling errors will cause the accumulation of estimation errors, affecting the algorithm’s performance and stability. Although the receding horizon strategy can correct some previous errors, it cannot eliminate the influence of modeling errors. Thus, the MPSP algorithm is highly dependent on the model’s accuracy. To the best of the authors’ knowledge, no paper has studied MPSP guidance in the presence of modeling errors so far.

Through the analysis of the ITACG method literature, we can draw the following conclusions:

(1)

Compared to non-predictive guidance, predictive guidance may present better performance for unpowered missile reality applications.

(2)

As one of the predictive guidance methods, MPSP-based guidance can avoid the model mismatch and thus can derive a better guidance performance, which has been verified in Refs. [21,22,23,24,25,26,27].

(3)

The existing MPSP-based guidance methods are mainly focused on the angle constraint only. The MPSP-based ITACG methods still need further research.

(4)

As an inevitable but significant factor affecting the MPSP algorithm, the modeling errors have not been considered so far.

Based on the above analysis, this paper focuses on the ITACG problems for unpowered missiles in the presence of modeling errors, trying to fill existing research gaps.

3. Problem Description

A 3-D terminal guidance scenario of a missile attacking a stationary ground target is considered in this paper, as shown in Figure 1. $M$ is the missile whose position is $\left(x_0,y_0,z_0\right)$. $T$ is a stationary ground target whose position is $\left(x_t,y_t,z_t\right)$. The missile is supposed to arrive at the target point with the desired impact path angle ${\theta }_f$, impact azimuth angle ${\psi }_{vf}$, and impact time $t_f$.

The dynamic model of the missile with an unknown modeling error is written as

$$\dot{X}=f\left(X,U,d\right)=\left[\begin{array}{l}V\cos \theta \cos {\psi }_v\hfill \\ V\sin \theta \hfill \\ -V\cos \theta \sin {\psi }_v\hfill \\ -D\left(\alpha ,\beta \right)-g\sin \theta +d_x\hfill \\ \frac{L\left(\alpha \right)-g\cos \theta }{V}+d_y\hfill \\ -\frac{Z\left(\beta \right)}{V\cos \theta }+d_z\hfill \end{array}\right],$$

(1)

where $X={\left[x,y,z,V,\theta ,{\psi }_v\right]}^{\mathrm{T}}$ is the state vector. x, y, and z are 3-D position coordinates of the missile. V, $\theta$, and ${\psi }_v$ are the missile’s velocity, path angle, and azimuth angle, respectively. The angle of attack (AOA) $\alpha$ and the side slip angle (SSA) $\beta$ compose the control vector $U={\left[\alpha ,\beta \right]}^{\mathrm{T}}$. The drag acceleration D, lift acceleration L, and lateral acceleration Z are written as

$$\left\{\begin{array}{l}D=C_x\left(\alpha ,\beta \right)q{S}_{ref}/m\\ L=C_y\left(\alpha \right)q{S}_{ref}/m\\ Z=C_z\left(\beta \right)q{S}_{ref}/m\end{array}\right..$$

(2)

In Equation (2), aerodynamic coefficients C_x, C_y, and C_z are related to AOA and SSA. q is dynamic pressure, S_ref is reference area, m is the missile’s mass. $d={\left[d_x,d_y,d_z\right]}^{\mathrm{T}}$ is the unknown modeling error, which may be caused by unmodeled dynamics, uncertain parameters, external disturbances, etc.

Selecting the output vector as $Y={\left[\begin{array}{lllll}x\hfill & y\hfill & z\hfill & \theta \hfill & {\psi }_v\hfill \end{array}\right]}^{\mathrm{T}}$, the purpose of our guidance method is to determine proper control commands $U$, making sure $Y\to Y_d$ when $t\to t_f$, where $Y_d={\left[\begin{array}{lllll}{x}_t\hfill & y_t\hfill & z_t\hfill & {\theta }_f\hfill & {\psi }_{vf}\hfill \end{array}\right]}^{\mathrm{T}}$.

4. CKF-MPSP Terminal Guidance Method

A CKF-MPSP guidance method is proposed and applied in the terminal guidance scenario to achieve offset-free control in the presence of modeling errors. The CKF-MPSP method comprises three parts: nominal MPSP guidance, CKF modeling error estimation, and one-step modeling-error compensation. The nominal MPSP guidance generates a baseline guidance command ignoring modeling errors. The CKF modeling error estimation generates the guessed modeling error, and the one-step modeling-error compensation introduces an error compensation term in the baseline guidance command to maintain precision and stability.

The schematic of the CKF-MPSP guidance method is shown in Figure 2. The nominal MPSP guidance method is firstly used to obtain the nominal control command $U_n$ according to initial states $X_0$ and a guess value of control command $U_0$. Then, regarding the disturbance as an initial guess value ${\widehat{d}}_0$, the one-step modeling-error compensation method is utilized to eliminate the modeling error’s influence and generate the control command $U$, which is substituted into the dynamic model to update the states $X$. The CKF algorithm is utilized to generate the estimation of the states and the modeling error, which are denoted as $\widehat{X}$ and $\widehat{d}$, according to the measurement $\widehat{y}$ of the global navigation satellite system (GNSS). $\widehat{X}$, $\widehat{d}$, and $U$ are used as initial values to calculate the control vector in the next guidance period. Repeating the process until hit, the ITACG is achieved.

4.1. Nominal MPSP Guidance

Ignoring the unknown modeling error, the MPSP method [20] is used to generate the baseline control commands. The nominal dynamic model is represented as

$$\dot{X}=f\left(X,U\right)=\left[\begin{array}{l}V\cos \theta \cos {\psi }_v\hfill \\ V\sin \theta \hfill \\ -V\cos \theta \sin {\psi }_v\hfill \\ -D-g\sin \theta \hfill \\ \frac{L-g\cos \theta }{V}\hfill \\ -\frac{Z}{V\cos \theta }\hfill \end{array}\right].$$

(3)

Discretizing Equation (3), the discretized dynamic model can be described as

$$\left\{\begin{array}{l}{X}_{k+1}=F_k\left(X_k,U_k\right)=X_k+h{f}_k\left(X_k,U_k\right)\\ Y_k=C{X}_k\end{array}\right.,$$

(4)

where $k=1,2,3,\cdots ,N$ represents the time grids and h is the simulation step. The output matrix $C$ is shown as

$$C=\left[\begin{array}{llllll}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right].$$

(5)

After settling the desired terminal time $N$ and the desired terminal output $Y_d={\left[\begin{array}{lllll}{x}_t\hfill & y_t\hfill & z_t\hfill & {\theta }_f\hfill & {\psi }_{vf}\hfill \end{array}\right]}^{\mathrm{T}}$, the prediction output vector $Y_N$ can be obtained through Runge–Kutta integration using initial states and previous control commands. Then, the updated control commands can be obtained according to the deviation between $Y_N$ and $Y_d$.

Denote the terminal output deviation as $\Delta Y_N=Y_N-Y_d$. Expanding $\Delta Y_N$ at $Y_d$ and ignoring the high-order terms, we can obtain

$$\Delta Y_N\cong \mathrm{d}{Y}_N=\left(\frac{\partial Y_N}{\partial X_N}\right)\mathrm{d}{X}_N.$$

(6)

According to Equation (4), we can determine $\partial Y_N/\partial X_N=C$, and the following formula holds.

$$\mathrm{d}{X}_{k+1}=\left(\frac{\partial F_k}{\partial X_k}\right)\mathrm{d}{X}_k+\left(\frac{\partial F_k}{\partial U_k}\right)\mathrm{d}{U}_k,$$

(7)

where $\mathrm{d}{X}_k$ and $\mathrm{d}{U}_k$ are the change in state and control at the k-th step, respectively.

The partial derivative of $F_k$ with respect to $X_k$ and $U_k$ are shown as below:

$$\begin{array}{cc}\hfill \frac{\partial F_k}{\partial X_k}& =I_{6\times 6}+h\frac{\partial f_k}{\partial X_k}\hfill \\ & =I_{6\times 6}+h\left[\begin{array}{llllll}0\hfill & 0\hfill & 0\hfill & \cos \theta \cos {\psi }_v\hfill & -V\sin \theta \cos {\psi }_v\hfill & -V\cos \theta \sin {\psi }_v\hfill \\ 0\hfill & 0\hfill & 0\hfill & \sin \theta \hfill & V\cos \theta \hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & -\cos \theta \sin {\psi }_v\hfill & V\sin \theta \sin {\psi }_v\hfill & -V\cos \theta \cos {\psi }_v\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & -g\cos \theta \hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & -\frac{L-g\cos \theta }{V^2}\hfill & \frac{g\sin \theta }{V}\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & \frac{Z}{V^2\cos \theta }\hfill & -\frac{Z\sin \theta }{V{\cos }^2\theta }\hfill & 0\hfill \end{array}\right]\hfill \end{array},$$

(8)

$$\frac{\partial F_k}{\partial U_k}=h\frac{\partial f_k}{\partial U_k}=h\left[\begin{array}{ll}0\hfill & 0\hfill \\ 0\hfill & 0\hfill \\ 0\hfill & 0\hfill \\ \frac{-C_x^{\alpha }q{S}_{ref}}{m}\hfill & \frac{-C_x^{\beta }q{S}_{ref}}{m}\hfill \\ \frac{C_y^{\alpha }q{S}_{ref}}{mV}\hfill & 0\hfill \\ 0\hfill & -\frac{C_z^{\beta }q{S}_{ref}}{mV\cos \theta }\hfill \end{array}\right],$$

(9)

where $C_x^{\alpha }$ is the partial derivative of the drag coefficient with respect to AOA, $C_y^{\alpha }$ is the derivative of the lift coefficient with respect to AOA, $C_x^{\beta }$ is the partial derivative of the drag coefficient with respect to SSA, and $C_z^{\beta }$ is the derivative of the lateral coefficient with respect to SSA. They can be obtained from the aerodynamic data.

Substituting Equation (7) into Equation (6), we can obtain

$$\mathrm{d}{Y}_N=\left(\frac{\partial Y_N}{\partial X_N}\right)\left[\left(\frac{\partial F_{N-1}}{\partial X_{N-1}}\right)\mathrm{d}{X}_{N-1}+\left(\frac{\partial F_{N-1}}{\partial U_{N-1}}\right)\mathrm{d}{U}_{N-1}\right].$$

(10)

In Equation (10), $\mathrm{d}{X}_{N-1}$ can be expressed as

$$\mathrm{d}{X}_{N-1}=\left(\frac{\partial F_{N-2}}{\partial X_{N-2}}\right)\mathrm{d}{X}_{N-2}+\left(\frac{\partial F_{N-2}}{\partial U_{N-2}}\right)\mathrm{d}{U}_{N-2}$$

(11)

And $\mathrm{d}{X}_{N-2}$ can be further expressed by $\mathrm{d}{X}_{N-3}$ and $\mathrm{d}{U}_{N-3}$. Repeating the above process until $\mathrm{d}{X}_1$ and $\mathrm{d}{U}_1$, it is clear that Equation (10) can be rewritten as

$$\mathrm{d}{Y}_N=A\mathrm{d}{X}_1+B_1\mathrm{d}{U}_1+B_2\mathrm{d}{U}_2+\cdots +B_{N-1}\mathrm{d}{U}_{N-1},$$

(12)

where

$$\begin{array}{l}A\triangleq \frac{\partial Y_N}{\partial X_N}\frac{\partial F_{N-1}}{\partial X_{N-1}}\frac{\partial F_{N-2}}{\partial X_{N-2}}\cdots \frac{\partial F_1}{\partial X_1}\\ \left\{\begin{array}{ll}{B}_k\triangleq \frac{\partial Y_N}{\partial X_N}\frac{\partial F_{N-1}}{\partial X_{N-1}}\frac{\partial F_{N-2}}{\partial X_{N-2}}\cdots \frac{\partial F_{k+1}}{\partial X_{k+1}}\frac{\partial F_k}{\partial U_k},\hfill & k=1,2,\cdots ,N-2\hfill \\ B_k\triangleq \frac{\partial Y_N}{\partial X_N}\frac{\partial F_{N-1}}{\partial U_{N-1}},\hfill & k=N-1\hfill \end{array}\right.\end{array}.$$

(13)

Since the initial states are with no errors ($\mathrm{d}{X}_1=0$), the final output error is only decided by control commands as

$$\mathrm{d}{Y}_N={\displaystyle \sum _{k=1}^{N-1}{B}_k\mathrm{d}{U}_k}.$$

(14)

The purpose of guidance is to find a series of control commands $U_k=U_k^0-\mathrm{d}{U}_k\left(k=1,2,\cdots ,N\right)$ to make $\mathrm{d}{Y}_N\to 0$, where $U_k^0$ is the previous control history solution. It is worth noting that Equation (14) has $2\times \left(N-1\right)$ unknowns and 5 equations. Usually, $2\times \left(N-1\right)>5$; thus, the solutions are not unique. To maximize guidance performance, set the following energy-optimal performance index and aim to minimize it.

$$J=\frac{1}{2}{\displaystyle \sum _{k=1}^{N-1}{\left(U_k^0-\mathrm{d}{U}_k\right)}^{\mathrm{T}}{R}_k\left(U_k^0-\mathrm{d}{U}_k\right)},$$

(15)

where $R_k$ is a positive definite weight coefficient matrix.

Equations (14) and (15) constitute a static optimization problem, whose solution at every time step $k=1,2,\cdots ,N$, according to static optimization theory, is

$$U_k^{*}=U_k^0-\mathrm{d}{U}_k=R_k^{-1}{B}_k^{\mathrm{T}}{A}_{\lambda }^{-1}\left(\mathrm{d}{Y}_N-b_{\lambda }\right),$$

(16)

where $A_{\lambda }\triangleq -{\displaystyle \sum _{k=1}^{N-1}{B}_k{R}_k^{-1}{B}_k^T}$,$b_{\lambda }\triangleq {\displaystyle \sum _{k=1}^{N-1}{B}_k{U}_k^0}$.

4.2. Modeling Error Estimation Based on CKF

The MPSP guidance method highly relies on modeling accuracy because of the integral prediction. However, the realistic model inevitably has unknown modeling errors or external disturbances. It has been pointed out in the literature [28] that, in the presence of model mismatch, the MPSP method cannot realize the desired terminal states. Estimating and compensating for modeling errors are common ways to achieve offset-free terminal state control. This section uses the CKF algorithm to estimate states and modeling errors simultaneously for subsequent compensation.

To estimate the modeling errors, consider them as constants and extend them to states. The dynamic model (1) can be rewritten as

$${\dot{X}}^E=f_E\left(X^E\right)=\left[\begin{array}{l}V\cos \theta \cos {\psi }_v\hfill \\ V\sin \theta \hfill \\ -V\cos \theta \sin {\psi }_v\hfill \\ -D-g\sin \theta +d_x\hfill \\ \frac{L-g\cos \theta }{V}+d_y\hfill \\ -\frac{Z}{V\cos \theta }+d_z\hfill \\ 0\hfill \\ 0\hfill \\ 0\hfill \end{array}\right]+w,$$

(17)

where $X^E={\left[x,y,z,V,\theta ,{\psi }_v,d_x,d_y,d_z\right]}^{\mathrm{T}}$ is the expansion state vector, $w$ is Gaussian-distributed process noise, and $E\left[w{w}^T\right]=Q$.

During flight, GNSS measures the missile’s motion in real time. So, the measurement equations can be denoted as

$${\widehat{y}}^E=X^E+v,$$

(18)

where $v$ is Gaussian-distributed measurement noise and $E\left[v{v}^T\right]=R$. ${\widehat{y}}^E={\left[\widehat{x},\widehat{y},\widehat{z},\widehat{V},\widehat{\theta },{\widehat{\psi }}_v,{\widehat{d}}_x,{\widehat{d}}_y,{\widehat{d}}_z\right]}^{\mathrm{T}}$ is the expansion output vector.

Discretizing Equations (17) and (18), we can obtain the nonlinear filter model:

$$\left\{\begin{array}{l}{X}_k^E=F_E\left(X_{k-1}^E\right)+w_{k-1}\hfill \\ Y_k^E=X_k^E+v_{k-1}\hfill \end{array}\right.,$$

(19)

where $F_E\left(X_{k-1}^E\right)=X_{k-1}^E+h{f}_E\left(X_{k-1}^E\right)$.

The CKF algorithm consists of two procedures: Time Update and Measurement Update [29]. Combined with the filter model, the CKF algorithm process is shown below.

4.2.1. Time Update

Assume at time k that the posterior probability density function $p\left(\left.X_{k-1}^E\right|y_{k-1}\right)=N\left({\widehat{X}}_{k-1}^E,P_{k-1}\right)$ is known. Denote the Cholesky factorization of the error covariance $P_{k-1}$ as $S_{k-1}$.

$$P_{k-1}=S_{k-1}{S}_{k-1}^{\mathrm{T}}.$$

(20)

Calculate the cubature points ${\chi }_{k-1}^{\left(i\right)}$ based on the third-degree cubature rule:

$$\begin{array}{cc}{\chi }_{k-1}^{\left(i\right)}={\widehat{X}}_{k-1}^E+S_{k-1}{\xi }_i,& i=1,2,\cdots ,2n\end{array},$$

(21)

where n is the dimension of states. ${\xi }_i=\sqrt{n}{\left[1\right]}_i$ is the basic cubature point set. The point set $\left[1\right]$ is defined as

$$\left[1\right]=\underset{2n}{\underbrace{\left\{{\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right]}_{n\times 1},\left[\begin{array}{c}0\\ 1\\ ⋮\\ 0\end{array}\right],\cdots ,\left[\begin{array}{c}0\\ 0\\ ⋮\\ 1\end{array}\right],\left[\begin{array}{c}-1\\ 0\\ ⋮\\ 0\end{array}\right],\left[\begin{array}{c}0\\ -1\\ ⋮\\ 0\end{array}\right],\cdots ,\left[\begin{array}{c}0\\ 0\\ ⋮\\ -1\end{array}\right]\right\}}},$$

(22)

and ${\left[1\right]}_i$ represents the i-th column vector in $\left[1\right]$.

Calculate the one-step prediction at time k and its error covariance:

$${\chi }_{\left.k\right|k-1}^{*\left(i\right)}=F_E\left({\chi }_{k-1}^{\left(i\right)}\right),$$

(23)

$${\widehat{X}}_{\left.k\right|k-1}^E=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}{\chi }_{\left.k\right|k-1}^{*\left(i\right)}},$$

(24)

$$P_{\left.k\right|k-1}=\frac{1}{2n}{{\displaystyle \sum _{i=1}^{2n}\left[{\chi }_{\left.k\right|k-1}^{*\left(i\right)}-{\widehat{X}}_{\left.k\right|k-1}^E\right]\left[{\chi }_{\left.k\right|k-1}^{*\left(i\right)}-{\widehat{X}}_{\left.k\right|k-1}^E\right]}}^{\mathrm{T}}+Q_{k-1}.$$

(25)

4.2.2. Measurement Update

Calculate the cubature points for Measurement Update, and then calculate measurement prediction:

$$P_{\left.k\right|k-1}=S_{\left.k\right|k-1}{S}_{\left.k\right|k-1}^{\mathrm{T}},$$

(26)

$$\begin{array}{cc}{\chi }_{\left.k\right|k-1}^{\left(i\right)}={\widehat{X}}_{\left.k\right|k-1}^E+S_{\left.k\right|k-1}{\xi }_i,& i=1,2,\cdots ,2n\end{array},$$

(27)

$${\widehat{y}}_{\left.k\right|k-1}^E=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}{\chi }_{\left.k\right|k-1}^{\left(i\right)}}.$$

(28)

Calculate the innovation covariance matrix $P_{xy}$ and the cross-covariance matrix $P_{yy}$:

$$P_{xy}=\frac{1}{2n}{{\displaystyle \sum _{i=1}^{2n}\left[{\chi }_{\left.k\right|k-1}^{\left(i\right)}-{\widehat{X}}_{\left.k\right|k-1}^E\right]\left[{\chi }_{\left.k\right|k-1}^{\left(i\right)}-{\widehat{y}}_{\left.k\right|k-1}^E\right]}}^{\mathrm{T}},$$

(29)

$$P_{yy}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}\left[{\chi }_{\left.k\right|k-1}^{\left(i\right)}-{\widehat{y}}_{\left.k\right|k-1}^E\right]{\left[{\chi }_{\left.k\right|k-1}^{\left(i\right)}-{\widehat{y}}_{\left.k\right|k-1}^E\right]}^{\mathrm{T}}}+R_k.$$

(30)

Calculate the Kalman gain:

$$K_k=P_{xy}{P}_{yy}^{-1}.$$

(31)

Estimate the updated state and the corresponding error covariance:

$${\widehat{X}}_k^E={\widehat{X}}_{\left.k\right|k-1}^E+K_k\left(y_k^E-{\widehat{y}}_{\left.k\right|k-1}^E\right),$$

(32)

$$P_k=P_{\left.k\right|k-1}-K_k{P}_{yy}{K}_k^{\mathrm{T}}.$$

(33)

4.3. One-Step Modeling-Error Compensation

After estimating the modeling errors, compensate for the effect of modeling errors on the system by attaching an additional control term $\Delta U_k$ to the MPSP optimal command $U_k^{*}$ at time k. Denoting the estimation modeling error vector as ${\widehat{d}}_k={\left[{\widehat{d}}_{xk},{\widehat{d}}_{yk},{\widehat{d}}_{zk}\right]}^{\mathrm{T}}$, the disturbed system model can be described as

$$\left\{\begin{array}{l}{\tilde{X}}_{k+1}={\tilde{X}}_k+h{f}_k\left({\tilde{X}}_k,{\tilde{U}}_k,{\widehat{d}}_k\right)\\ {\tilde{Y}}_{k+1}=C{\tilde{X}}_{k+1}\end{array}\right..$$

(34)

At time k, denote $Y_{k+1}$ as the output vector obtained by a one-step calculation with current state $X_k$ and the MPSP control command $U_k^{*}$. The objective of modeling-error compensation is to generate a modified control command ${\tilde{U}}_k=U_k^{*}+\Delta U_k$ making ${\tilde{Y}}_{k+1}\to Y_{k+1}$.

Denote the output error as

$$\Delta Y_{k+1}={\tilde{Y}}_{k+1}-Y_{k+1}.$$

(35)

Substituting Equations (4) and (34) into (35) and because of ${\tilde{X}}_k=X_k$, we can obtain

$$\Delta Y_{k+1}=hC\left[f_k\left(X_k,U_k^{*}+\Delta U_k,{\widehat{d}}_k\right)-f_k\left(X_k,U_k^{*}\right)\right].$$

(36)

Expanding $f_k\left(X_k,U_k^{*}+\Delta U_k,{\widehat{d}}_k\right)$ and ignoring high-order terms, we can obtain

$$f_k\left(X_k,U_k^{*}+\Delta U_k,{\widehat{d}}_k\right)\cong f_k\left(X_k,U_k^{*}\right)+\frac{\partial f_k}{\partial U_k^{}}\Delta U_k+\frac{\partial f_k}{\partial d_k^{}}{\widehat{d}}_k.$$

(37)

Substituting Equation (37) into (36), we can obtain

$$\Delta Y_{k+1}=hC\left[\frac{\partial f_k}{\partial U_k^{}}\Delta U_k+\frac{\partial f_k}{\partial d_k^{}}{\widehat{d}}_k\right].$$

(38)

The compensation term is desired to make the output error zero. According to Equation (38), we can obtain the desired additional control term as

$$\Delta U_k=-{\left(C\frac{\partial f_k}{\partial U_k^{}}\right)}^{-1}\left(C\frac{\partial f_k}{\partial d_k^{}}{\widehat{d}}_k\right).$$

(39)

The modified control command at time k is

$${\tilde{U}}_k=U_k^{*}+\Delta U_k.$$

(40)

4.4. CKF-MPSP Terminal Guidance Process

Considering the above, the CKF-MPSP guidance process is summarized below (Algorithm 1).

Algorithm 1: CKF-MPSP Terminal Guidance Scheme.

INPUT: current time k, desired terminal time N, desired terminal output $Y_d$, current states $X_k$, guessed control command $\left[U_k^0,U_{k+1}^0,\cdots ,U_{N-1}^0\right]$, estimated modeling error ${\widehat{d}}_k^{}$, CKF initial values ${\widehat{X}}_{k-1}^E$, ${\widehat{P}}_{k-1}$. 1: while current time k is no larger than N, do 2:  while terminal output deviation $\Delta Y_N$ is larger than tolerance value $\epsilon$, do 3:   predict the terminal output vector $Y_N$ through Runge-Kutta integration with the nominal dynamic model (4), the current states $X_k$ and the guessed control command $\left[U_k^0,U_{k+1}^0,\cdots ,U_{N-1}^0\right]$. 4:   calculate terminal output deviation $\Delta Y_N=Y_N-Y_d$. 5:   calculate matrices $\left[B_k^{},B_{k+1}^{},\cdots ,B_{N-1}^{}\right]$ according to Equation (13). 6:   calculate the optimal control command $\left[U_k^{*},U_{k+1}^{*},\cdots ,U_{N-1}^{*}\right]$ according to Equation (16). 7:   take $\left[U_{k+1}^{*},\cdots ,U_{N-1}^{*}\right]$ as the new guessed control command. 8:  end while. 9:  calculate the one-step output ${\tilde{Y}}_{k+1}$ with the disturbed system model (34) in presence of the estimated modeling error ${\widehat{d}}_k^{}$ 10: calculate the modified control command ${\tilde{U}}_k$ at time k, according to Equations (39) and (40). 11:  substitute ${\tilde{U}}_k$ into the realistic dynamic model (1) and obtain the updated state $X_{k+1}$. 12:  estimate the filter state ${\widehat{X}}_k^E$, error covariance ${\widehat{P}}_k$, and modeling error ${\widehat{d}}_{k+1}^{}$, using CKF algorithm (20)~(33). 13: time update, k = k + 1. 14: end while

Remark 1.

 The MPSP guidance method takes the desired impact time as the terminal time N of static planning. The target’s position and impact angles are regarded as desired terminal output $Y_d$. Making the terminal output deviation $\Delta Y_N$ no larger than tolerance value $\epsilon$ by iterating $U^{*}$, the impact time and angle constraints can be satisfied simultaneously.

Remark 2.

 The MPSP guidance method relies on initial guessed control commands $\left[U_k^0,U_{k+1}^0,\cdots ,U_{N-1}^0\right]$. The guessed control commands are quickly generated through some simple guidance laws in common. In this paper, the traditional PN guidance law is used to obtain the initial guessed control commands. Usually, the impact time of PN, represented by the symbol $N_P$, is different from the desired time N. For the MPSP algorithm, $N_P$ must be no less than N, so the outputs at time N can be predicted. In this paper, a protection mechanism is introduced to make sure the MPSP algorithm normally runs even if $N_P$ is less than N: Let $U_k^0=U_{N_P-1}^0$, for $N_P-1<k\le N-1$.

5. Simulations and Results

In this section, several numerical simulations are carried out to evaluate the performance of the proposed terminal guidance method in the presence of modeling errors. A three-dimensional guidance scenario of a missile attacking a stationary target is constructed. The initial simulation conditions are listed below (

Table 1. Simulation initial conditions.

Parameters	Values
Missile’s initial velocity $V$	200 m/s
Missile’s initial path angle $\theta$	0°
Missile’s initial azimuth angle ${\psi }_v$	0°
Missile’s initial position $\left(x,y,z\right)$	(0 m, 4000 m, 0 m)
Target’s position $\left(x_t,y_t,z_t\right)$	(5000 m, 0 m, 1000 m)

):

Besides achieving precision arrival, the impact time and impact angles are also required. The desired impact time is settled as 35 s, and the desired terminal path angle and azimuth angle are settled as −80° and −50°, respectively. The dynamic modeling errors in Equation (1) are settled as $d_x=-0.4\sin \left(t/400\right)\cos \left(t/400\right)$, $d_y=1$, and $d_z=2.5\cos \left(\pi t/200\right)$. The measurement errors from the GNSS system are assumed to be normally distributed. The position error, velocity error, and acceleration error are settled to be 10 m (3σ), 1 m/s (3σ), and 0.1 m/s² (3σ), respectively.

For the MPSP algorithm, the guessed control commands $\left[U_k^0,U_{k+1}^0,\cdots ,U_{N-1}^0\right]$ are needed. The traditional PN [30] is used to produce the initial values of $\left[U_k^0,U_{k+1}^0,\cdots ,U_{N-1}^0\right]$ in this paper, and the navigation ratio is settled as 6. In addition, a comparison with the MPSP guidance method presented in [21] is provided to validate the superiority of the method. The end condition for PN simulation is that the missile reaches the target. And the end conditions for the other two simulations are that the simulation times reach the desired impact time, which is 35 s on this occasion. For MPSP and CKF-MPSP methods, the tolerance value vector is settled as $\epsilon ={\left[1\mathrm{m},1\mathrm{m},1\mathrm{m},{0.1}^{\circ },{0.1}^{\circ }\right]}^{\mathrm{T}}$. The simulation results are shown below.

To evaluate the guidance accuracy, some crucial parameters of the three methods shown in Figure 3 are provided in

Table 2. Simulation results of PN method.

Method	Miss Distance (m)	Velocity (m/s)	Path Angle (°)	Azimuth Angle (°)	Impact Time (s)
PN	0.26	242.56	−46.42	−15.23	30.5

and

Table 3. Simulation results of MPSP and CKF-MPSP methods at the desired impact time (35 s).

Method	Miss Distance (m)	Velocity (m/s)	Path Angle (°)	Azimuth Angle (°)
MPSP	289.73	161.18	−90.24	218.31
CKF-MPSP	0.29	249.90	−79.99	−50.00

, which include the terminal miss distance, terminal velocity, terminal path angle, terminal azimuth angle, and impact time.

It is obvious that the CKF-MPSP method has good accuracy for miss distance while strictly constraining the terminal path angle, azimuth angle, and impact time. The PN method leads to the minimum miss distance. However, it cannot consider impact time and angle constraints. Because of the dynamic modeling errors, the MPSP method cannot find a feasible solution, which leads to a significant guidance error. Affected by modeling errors, the MPSP method’s velocity is lower than the CKF-MPSP’s at every identical moment, which is also the main reason for the MPSP method not reaching the destination. The CKF-MPSP method estimates the modeling errors with high accuracy. Referring to Figure 4, the estimation error of modeling errors is no larger than 0.04 m/s². The influence of modeling errors can be reduced by compensating for the nominal command acceleration. The simulation results illustrate the effectiveness and superiority of the CKF-MPSP method in the presence of modeling errors.

Based on the above simulation, the influence of atmospheric density deviation, aerodynamic parameter deviation, and the random variation in the dynamic modeling errors are further considered. The deviations are assumed to be normally distributed, and their values are 10% (3σ). The results of 200 Monte Carlo simulations are as follows.

The key indexes in Figure 5 is summarized in

Table 4. Monte Carlo simulation results.

Method	Miss Distance (m)	Velocity (m/s)	Path Angle (°)	Azimuth Angle (°)
Average	1.21	248.92	−79.95	−49.93
Maximum value	5.56	250.11	−79.21	−48.49
Minimum value	0.03	240.40	−80.12	−50.26
Standard error	0.95	1.84	0.16	0.34

. According to Table 4, the miss distances are no larger than 5.56 m and the average impact angle errors are 0.055° and 0.077°, respectively. The missile maintains high guidance accuracy in the presence of disturbances. The simulation results illustrate that the proposed guidance method has strong robustness.

6. Conclusions

ITACG is a vital field for missiles because it can improve destructive effects and hit weak parts of time-sensitive targets, and make it possible for multiple missiles to attack a target simultaneously. Studying ITACG can effectively enhance the combat effectiveness of missiles.

In this paper, an ITACG is proposed based on the MPSP algorithm. By taking the desired impact time and angles as terminal conditions, the guidance method can satisfy these constraints simultaneously. Furthermore, to eliminate the influence of modeling errors on prediction, the CKF algorithm is used for error estimation, and a compensation scheme is designed. The proposed guidance method considers the missile’s dynamic model instead of a constant-velocity model. Meanwhile, the modeling errors are estimated and compensated. Thus, this method is more practically significant. A terminal guidance scenario is settled, and the PN method, MPSP method, and CKF-MPSP method are used for simulation in the presence of modeling errors. According to the simulation results, the CKF-MPSP method can achieve impact time and angle constraint guidance, and maintain high accuracy within the influence of modeling errors. Furthermore, the Monte Carlo simulation is conducted, considering the influence of atmospheric density deviations, aerodynamic parameter deviations, and random variations in the dynamic modeling errors. According to the simulation results, the miss distances are no larger than 5.56 m, and the average impact angle errors are 0.055° and 0.077°, respectively. The missile maintains high guidance accuracy in the presence of disturbances. Comprehensively, the simulation results illustrate that the proposed CKF-MPSP guidance method has high precision and strong robustness.

It should be acknowledged that this article still has some limitations. The research of this paper is mainly focused on the ITACG against stationary targets. For moving targets, it is also necessary to introduce their motion models into the guidance method. However, the accurate estimation of the moving targets’ motion is still a difficult problem, because of insufficient target information and potential maneuvering. Motion model estimation and MPSP-based guidance for moving targets remain to be researched in the future.