Study of Model Uncertainties Influence on the Impact Point Dispersion for a Gasodynamicaly Controlled Projectile

Abstract

The article presents the analysis of the impact point dispersion reduction using lateral correction thrusters. Two types of control algorithms are used and four sources of uncertainties are taken into account: aerodynamic parameters, thrust curve, initial conditions and IMU errors. The Monte Carlo approach was used for simulations and Circular Error Probable was used as a measure of dispersion. Generic rocket mathematical and simulation model was created in MATLAB/Simulink 2020b environment. Results show that the use of control algorithms greatly reduces the impact point dispersion.

1. Introduction

Determining and lowering the impact point dispersion of an artillery projectile is an important factor when determining its usefulness and effectiveness. Due to the imperfections of modelling of the flight of such objects, unknowns, simplifications, and uncertainties in model parameters, it is important to take into account the various possibilities of scenarios and determine and quantify the dispersion of possible landing points. Key factors influencing the flight of the object, and the hardest to obtain accurately for modelling, are the aerodynamic characteristics, thrust curve, and initial conditions. To maximize the accuracy of the projectile, various types of control methods, algorithms, and target detection methods are utilized, such as H${}_{\infty }$ guidance law [1], proportional navigation guidance and its modifications [2,3], various optimal control methods [4,5,6,7] or model predictive approach [8].

In [9,10] the impact point dispersion of a lateral pulse jet controlled rocket following a reference trajectory and its robustness to the effects of the measurement noise were studied. Launch conditions’ uncertainties were studied in [11,12]. In [13,14] the impact point dispersion due to manufacturing errors using Monte Carlo method was studied. In [15] the effects of launch condition variability, atmospheric factors, and IMU errors on the guidance accuracy were investigated. The influence of missile initial conditions’ uncertainties and IMU errors on the impact point dispersion using lateral thrusters for control were investigated in [16] and the influence of uncertainties in rocket parameters on the performance of a cold launch were analyzed in [17]. IMU errors and noise impact on the guidance were also investigated in [18,19]. The analysis of the robustness of the control algorithm with respect to uncertainty regarding the launch environment and rocket conditions was presented in [20,21]. Monte Carlo analysis of the impact point dispersion due to the missile parameters and atmospheric conditions’ uncertainties was performed in [22] and the impact point dispersion reduction due to high spin motion was analyzed in [23].

In this paper, the analysis of the impact point dispersion caused by model uncertainties is presented. The uncertainties in the aerodynamic parameters, thrust curve, and initial conditions, as well as the uncertainties caused by the Inertial Measurement Unit model and its errors are analyzed. Two types of control algorithms were tested, Multi-Condition Control Algorithm (MCCA) and modified Proportional Navigation Guidance (mPNG). The analysis was performed using the Monte Carlo approach.

The prepared article is organised as follows: in Section 2 the mathematical model of the rocket is presented in the Section 2.1 with used assumptions and coordinate systems. The dynamic equations of motion, followed by additional kinematic equations are shown in that section. The external loads comprised of aerodynamics, gravity, propulsion, and correction thrusters are presented, and this is followed by the description of inertial parameters, atmosphere, and Inertial Measurement Unit modelling approach. At last, the used control algorithms are shown. In Section 2.2 the simulational model created in MATLAB/Simulink R2020b is presented. Section 3 presents the simulational study and its results, followed by the discussion and interpretation. The paper finishes with a short summary of the conclusions shown in Section 4.

The novelty in this paper is the comparison of the IMU errors’ influence on the accuracy of the two control algorithms utilizing lateral pulse jet control.

The developed numerical simulation might be used in design process of new sensors intended for projectile navigation. This tool might be used for fast prototyping of control schemes that reduces the overall system design process and costs.

2. Materials and Methods

2.1. Mathematical Model

2.1.1. Assumptions

For creating a mathematical and simulational model of a generic rocket, several assumptions were made. The rocket is modelled as a rigid body with six degrees of freedom and variable inertial parameters. It is controlled by a set of solid rocket motor thrusters, which use does not change the inertial and aerodynamic parameters of the rocket. Atmosphere model is taken from the International Standard Atmosphere [24]. Earth rotation and eccentricity are not modelled, the gravitational acceleration is constant and consistent with WGS-84 model [25]. The mathematical model includes the Inertial Measurement Unit (IMU) comprised of accelerometers and gyroscopes triades, modelled as second-order dynamical systems, with their disturbances that include noise, bias, scale factor and cross coupling as well as g-dependent factor for gyroscopes.

2.1.2. Coordinate Systems

The coordinate systems used are presented in Figure 1:

• The navigational coordinate system $O_n{x}_n{y}_n{z}_n$ is a right-handed, Cartesian coordinate system fixed to earth. Its origin is located at any point and the $O_n{x}_n{y}_n$ plane is tangent to the surface of the earth. The $O_n{x}_n$ axis points in the direction of a launch rail, $O_n{z}_n$ axis in the direction of the gravitational acceleration and $O_n{y}_n$ axis completes the right-handed coordinate system.

• The gravitational coordinate system $O_g{x}_g{y}_g{z}_g$ is a right-handed, Cartesian moving coordinate system fixed with the rocket. Its origin is located at the center of mass of the rocket and the whole system remains parallel to the navigational system during the whole flight of the rocket.

• The body coordinate system $O_b{x}_b{y}_b{z}_b$ is a right-handed Cartesian coordinate system. Its origin is located at any point of the rocket. The $O_b{x}_b$ axis is parallel to the rockets longitudinal axis and points forward, $O_b{y}_b$ axis points at the right wing and $O_b{z}_b$ axis completes the right-handed coordinate system. Orientation of the body coordinate system with respect to the gravitational coordinate system is described by the Euler angles of yaw $\mathsf{\Psi }$, pitch $\mathsf{\Theta }$ and roll $\mathsf{\Phi }$.

• The measuring coordinate system $O_e{x}_e{y}_e{z}_e$ is a right-handed, Cartesian coordinate system fixed with the airflow. Its origin is located at any point of the $O_b{y}_b{z}_b$ symmetry plane and the position of that point with respect to the $O_b{x}_b{y}_b{z}_b$ coordinate system is defined by the vector ${\mathbf{r}}_{\mathbf{e}}$. The $O_e{x}_e$ axis lies in the direction of the airflow, $O_e{z}_e$ axis points upwards and the $O_e{y}_e$ axis points to the right.

• The aerodynamic coordinate system $O_a{x}_a{y}_a{z}_a$ is the right-handed Cartesian coordinate system fixed with the airflow. The $O_a{x}_a$ axis points in the opposite direction that the airflow, $O_a{z}_a$ axis points downwards and the $O_a{y}_a$ completes the right-handed coordinate system.

2.1.3. Dynamic Equations of Motion

For developing the dynamic equations of motion, the linear and angular momentum change theorems for the rigid body were used. In the non-inertial frame $O_b{x}_b{y}_b{z}_b$ with the origin not located at the center of mass, they are given as [26,27]:

(1)$\frac{\tilde{\delta }\mathsf{\Pi }}{\tilde{\delta }t}+\mathsf{\Omega }\times \mathsf{\Pi }={\mathbf{F}}_{\mathbf{b}}$

(2)$\frac{\tilde{\delta }{\mathbf{K}}_{\mathbf{0}}}{\tilde{\delta }t}+\mathsf{\Omega }\times {\mathbf{K}}_{\mathbf{0}}+{\mathbf{V}}_{\mathbf{b}}\times \mathsf{\Pi }={\mathbf{M}}_{\mathbf{b}}$

where ${\mathbf{V}}_{\mathbf{b}}={\left[\begin{array}{ccc}U& V& W\end{array}\right]}^T$ is the velocity vector, $\mathsf{\Omega }={\left[\begin{array}{ccc}P& Q& R\end{array}\right]}^T$ is the angular velocity vector, ${\mathbf{F}}_{\mathbf{b}}={\left[\begin{array}{ccc}{X}_b& Y_b& Z_b\end{array}\right]}^T$ is the vector of external forces acting on the object, ${\mathbf{M}}_{\mathbf{b}}={\left[\begin{array}{ccc}{L}_b& M_b& N_b\end{array}\right]}^T$ is the vector of external torques with respect to point $O_b$ and $\frac{\tilde{\delta }}{\tilde{\delta }t}$ is the local derivative. Linear and angular momentum for a rigid body are [26]:

(3)$\mathsf{\Pi }=m\left({\mathbf{V}}_{\mathbf{b}}+\mathsf{\Omega }\times {\mathbf{r}}_{\mathbf{C}}\right)$

(4)${\mathbf{K}}_{\mathbf{0}}=\mathbf{I}\mathsf{\Omega }+{\mathbf{r}}_{\mathbf{C}}\times m{\mathbf{V}}_{\mathbf{b}}$

where m is the instantaneous mass of the rocket, $\mathbf{I}$ is the instantaneous moment of inertia tensor and ${\mathbf{r}}_{\mathbf{C}}$ is the center of mass position with respect to $O_b$. The dynamic equations of motion can be written as:

(5)$m{\dot{\mathbf{V}}}_{\mathbf{b}}+\dot{\mathsf{\Omega }}\times \mathbf{S}+\mathsf{\Omega }\times m{\mathbf{V}}_{\mathbf{b}}+\mathsf{\Omega }\times \left(\mathsf{\Omega }\times \mathbf{S}\right)={\mathbf{F}}_{\mathbf{b}}$

(6)$\mathbf{I}\dot{\mathsf{\Omega }}+\mathbf{S}\times {\dot{\mathbf{V}}}_{\mathbf{b}}+\dot{\mathbf{I}}\mathsf{\Omega }+\mathsf{\Omega }\times \mathbf{I}\mathsf{\Omega }+\mathsf{\Omega }\times \left(\mathbf{S}\times {\mathbf{V}}_{\mathbf{b}}\right)+{\mathbf{V}}_{\mathbf{b}}\times \left(\mathsf{\Omega }\times \mathbf{S}\right)={\mathbf{M}}_{\mathbf{b}}$

where $\mathbf{S}=m{\mathbf{r}}_{\mathbf{C}}$ is the first moment of mass. It must be noted that the propulsion terms resulting from expelling of the propellant by the main motor is included on the right side of above-mentioned equations. In the moments equations the jet damping effect was also omitted because for rocket artillery projectiles it is rather small when compared to aerodynamic damping. After changing cross-products to matrix multiplication, using the skew-symmetric matrix notation, ${\left[\right]}_x$, this can be written as:

(7)$\left[\begin{array}{cc}m\mathbf{1}& -{\left[\mathbf{S}\right]}_x\\ {\left[\mathbf{S}\right]}_x& \mathbf{I}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{V}}}_{\mathbf{b}}\\ \dot{\mathsf{\Omega }}\end{array}\right]+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& \dot{\mathbf{I}}\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_{\mathbf{b}}\\ \mathsf{\Omega }\end{array}\right]+\left[\begin{array}{cc}{\left[\mathsf{\Omega }\right]}_x& \mathbf{0}\\ [{\mathbf{V}}_{\mathbf{b}}{]}_x& {\left[\mathsf{\Omega }\right]}_x\end{array}\right]\left[\begin{array}{cc}m\mathbf{1}& -{\left[\mathbf{S}\right]}_x\\ {\left[\mathbf{S}\right]}_x& \mathbf{I}\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_{\mathbf{b}}\\ \mathsf{\Omega }\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_{\mathbf{b}}\\ {\mathbf{M}}_{\mathbf{b}}\end{array}\right]$

where $\mathbf{0}$ is the zero matrix and $\mathbf{1}$ is the unit matrix. In the short form this is:

(8)$\mathbf{A}\dot{\mathbf{x}}+\dot{\mathbf{A}}\mathbf{x}+\omega \mathbf{Ax}={\mathbf{F}}_{\mathbf{B}}$

where the state vector has the form $\mathbf{x}={\left[\begin{array}{cccccc}U& V& W& P& Q& R\end{array}\right]}^T$. This equation can be then numerically integrated to obtain the state vector.

2.1.4. Orientation

For determining object orientation, the quaternion algebra was used. Quaternion describes the orientation of the object in terms of rotation around a specific axis and is written as [28,29]:

(9)$\mathbf{e}=e_0+e_1\mathbf{i}+e_2\mathbf{j}+e_3\mathbf{k}$

where $e_0,e_1,e_2,e_3$ are the real numbers and $\mathbf{i},\mathbf{j},\mathbf{k}$ are the axes versors. The real parts of the quaternion can be written in terms of rotation axes direction cosines $E_x,E_y,E_z$ and rotation angle ${\delta }_E$ and are presented in Equations (10)–().

(10)$\begin{array}{ccc}\hfill e_0& =& \cos \frac{{\delta }_E}{2}\hfill \end{array}$

(11)$\begin{array}{ccc}\hfill e_1& =& E_x\sin \frac{{\delta }_E}{2}\hfill \end{array}$

(12)$\begin{array}{ccc}\hfill e_2& =& E_y\sin \frac{{\delta }_E}{2}\hfill \end{array}$

(13)$\begin{array}{ccc}\hfill e_3& =& E_z\sin \frac{{\delta }_E}{2}\hfill \end{array}$

The kinematic equation for the rate of change of the quaternion is given as [29,30]:

(14)$\left[\begin{array}{c}{\dot{e}}_0\\ {\dot{e}}_1\\ {\dot{e}}_2\\ {\dot{e}}_3\end{array}\right]=-\frac{1}{2}\left[\begin{array}{cccc}0& P& Q& R\\ -P& 0& -R& Q\\ -Q& R& 0& -P\\ -R& -Q& P& 0\end{array}\right]\left[\begin{array}{c}{e}_0\\ e_1\\ e_2\\ e_3\end{array}\right]-kE\left[\begin{array}{c}{e}_0\\ e_1\\ e_2\\ e_3\end{array}\right]$

where k is the feedback coefficient and E is the bounding equation violation coefficient $E={\left|\mathbf{e}\right|}^2-1$. It was assumed that $k=1$. Quaternions can be used to calculate the transformation matrix from the body to the navigation coordinate system as [29]:

(15)$\mathsf{\Lambda }=\left[\begin{array}{ccc}{e}_0^2+e_1^2-e_2^2-e_3^2& 2(e_1{e}_2-e_0{e}_3)& 2(e_0{e}_2+e_1{e}_3)\\ 2(e_0{e}_3+e_1{e}_2)& e_0^2-e_1^2+e_2^2-e_3^2& 2(e_2{e}_3-e_0{e}_1)\\ 2(e_1{e}_3-e_0{e}_2)& 2(e_0{e}_1+e_2{e}_3)& e_0^2-e_1^2-e_2^2+e_3^2\end{array}\right]$

and orientation angles of roll, pitch, and yaw, given as [29,30]:

(16)$\begin{array}{ccc}\hfill \mathsf{\Phi }& =& \arctan \frac{2(e_0{e}_1+e_2{e}_3)}{e_0^2-e_1^2-e_2^2+e_3^2}\hfill \end{array}$

(17)$\begin{array}{ccc}\hfill \mathsf{\Theta }& =& \arcsin 2(e_0{e}_2-e_1{e}_3)\hfill \end{array}$

(18)$\begin{array}{ccc}\hfill \mathsf{\Psi }& =& \arctan \frac{2(e_0{e}_3+e_1{e}_2)}{e_0^2+e_1^2-e_2^2-e_3^2}\hfill \end{array}$

Using the transformation matrix from (15) can be used in the second kinematic equation bounding linear velocities in $O_n$ coordinate system with velocities in $O_b$ coordinate system.

(19)$\left[\begin{array}{c}{\dot{x}}_n\\ {\dot{y}}_n\\ {\dot{z}}_n\end{array}\right]=\mathsf{\Lambda }\left[\begin{array}{c}U\\ V\\ W\end{array}\right]$

The initial quaternion can be determined from the initial orientation angles by means of equations [29,30,31]:

(20)$\begin{array}{ccc}\hfill e_0& =& \cos \frac{\mathsf{\Phi }}{2}\cos \frac{\mathsf{\Theta }}{2}\cos \frac{\mathsf{\Psi }}{2}+\sin \frac{\mathsf{\Phi }}{2}\sin \frac{\mathsf{\Theta }}{2}\sin \frac{\mathsf{\Psi }}{2}\hfill \end{array}$

(21)$\begin{array}{ccc}\hfill e_1& =& \sin \frac{\mathsf{\Phi }}{2}\cos \frac{\mathsf{\Theta }}{2}\cos \frac{\mathsf{\Psi }}{2}-\cos \frac{\mathsf{\Phi }}{2}\sin \frac{\mathsf{\Theta }}{2}\sin \frac{\mathsf{\Psi }}{2}\hfill \end{array}$

(22)$\begin{array}{ccc}\hfill e_2& =& \cos \frac{\mathsf{\Phi }}{2}\sin \frac{\mathsf{\Theta }}{2}\cos \frac{\mathsf{\Psi }}{2}+\sin \frac{\mathsf{\Phi }}{2}\cos \frac{\mathsf{\Theta }}{2}\sin \frac{\mathsf{\Psi }}{2}\hfill \end{array}$

(23)$\begin{array}{ccc}\hfill e_3& =& \cos \frac{\mathsf{\Phi }}{2}\cos \frac{\mathsf{\Theta }}{2}\sin \frac{\mathsf{\Psi }}{2}-\sin \frac{\mathsf{\Phi }}{2}\sin \frac{\mathsf{\Theta }}{2}\cos \frac{\mathsf{\Psi }}{2}\hfill \end{array}$

2.1.5. External Loads

The motion of the object is caused by the external forces and torques from aerodynamics ${\mathbf{F}}_{\mathbf{a}}$ and ${\mathbf{M}}_{\mathbf{a}}$, gravity ${\mathbf{F}}_{\mathbf{g}}$ and ${\mathbf{M}}_{\mathbf{g}}$, thrust ${\mathbf{F}}_{\mathbf{s}}$ and ${\mathbf{M}}_{\mathbf{s}}$ and reaction thrusters ${\mathbf{F}}_{\mathbf{sk}}$ and ${\mathbf{M}}_{\mathbf{sk}}$ [16]:

(24)$\begin{array}{ccc}\hfill {\mathbf{F}}_{\mathbf{b}}& =& {\mathbf{F}}_{\mathbf{a}}+{\mathbf{F}}_{\mathbf{g}}+{\mathbf{F}}_{\mathbf{s}}+{\mathbf{F}}_{\mathbf{sk}}\hfill \end{array}$

(25)$\begin{array}{ccc}\hfill {\mathbf{M}}_{\mathbf{b}}& =& {\mathbf{M}}_{\mathbf{a}}+{\mathbf{M}}_{\mathbf{g}}+{\mathbf{M}}_{\mathbf{s}}+{\mathbf{M}}_{\mathbf{sk}}\hfill \end{array}$

2.1.6. Aerodynamics

Aerodynamic force and moment vectors are given with respect to point $O_e$, so with respect to point $O_b$ they are given as:

(26)${\mathbf{F}}_{\mathbf{a}}=\left[\begin{array}{c}{X}_a\\ Y_a\\ Z_a\end{array}\right]$

(27)${\mathbf{M}}_{\mathbf{a}}=\left[\begin{array}{c}{L}_a\\ M_a\\ N_a\end{array}\right]={\mathbf{M}}_{\mathbf{a},O_e}+{\mathbf{r}}_{\mathbf{e}}\times {\mathbf{F}}_{\mathbf{a}}$

where ${\mathbf{r}}_{\mathbf{e}}={\mathbf{r}}_{\mathbf{we}}-{\mathbf{r}}_{\mathbf{wC}}+{\mathbf{r}}_{\mathbf{C}}$ is the vector describing the location of point $O_e$ with respect to point $O_b$, ${\mathbf{r}}_{\mathbf{w}}\mathbf{e}$ is the position of point $O_e$ with respect to rocket’s base and ${\mathbf{r}}_{\mathbf{wC}}$ is the position of point $O_b$ with respect to the rocket’s base. Aerodynamic force and moment are given as:

(28)${\mathbf{F}}_{\mathbf{a}}=\frac{1}{2}\rho {\left|{\mathbf{V}}_{\mathbf{b}}\right|}^{2}S\left[\begin{array}{c}{C}_X(\alpha ,\beta ,Ma)\\ C_Y(\alpha ,\beta ,Ma)\\ C_Z(\alpha ,\beta ,Ma)\end{array}\right]$

(29)${\mathbf{M}}_{\mathbf{a},O_e}=\frac{1}{2}\rho {\left|{\mathbf{V}}_{\mathbf{b}}\right|}^{2}Sd\left[\begin{array}{c}{C}_l(\alpha ,\beta ,Ma)\\ C_m(\alpha ,\beta ,Ma)\\ C_n(\alpha ,\beta ,Ma)\end{array}\right]$

where $\rho$ is the air density, S is the rocket’s cross section area and d is the rocket’s diameter. Aerodynamic incidence angles are given as [29,30]:

(30)$\begin{array}{ccc}\hfill \alpha & =& \arctan \frac{W}{U}\hfill \end{array}$

(31)$\begin{array}{ccc}\hfill \beta & =& \arcsin \frac{V}{|{\mathbf{V}}_{\mathbf{b}}|}\hfill \end{array}$

(32)$\begin{array}{ccc}\hfill Ma& =& \frac{|{\mathbf{V}}_{\mathbf{b}}|}{a}\hfill \end{array}$

where $\alpha$ is the angle of attack, $\beta$ is the sideslip angle, and $Ma$ is the Mach number given in (), where a is the local speed of sound. The aerodynamic coefficients are [9,32]:

(33)$\begin{array}{c}{C}_X=(C_{X_{base0}}+C_{X_{base{\alpha }^2}}{\alpha }^2+C_{X_{base{\beta }^2}}{\beta }^2)+(C_{X_{eng0}}+C_{X_{eng{\alpha }^2}}{\alpha }^2+C_{X_{eng{\beta }^2}}{\beta }^2){\delta }_{eng}\\ C_Y=C_{Y_0}+C_{Y_{\beta }}\beta \\ C_Z=C_{Z_0}+C_{Z_{\alpha }}\alpha \\ C_l=C_{l_0}+(C_{l_{p_0}}+C_{l_{p_{{\alpha }^2}}}{\alpha }^2+C_{l_{p_{{\beta }^2}}}{\beta }^2)\frac{Pd}{2|{\mathbf{V}}_{\mathbf{b}}|}\\ C_m=C_{m_0}+C_{m_{\alpha }}\alpha \\ C_n=C_{n_0}+C_{n_{\beta }}\beta \end{array}$

where $C_{X_0}$ is zero-longitudinal axial force coefficient, $C_{Y_{\beta }}$ is side force with angle of sideslip derivative, $C_{Z_{\alpha }}$ is normal force with respect to angle of attack derivative, $C_{l_0}$ is spin driving rolling moment coefficient and $C_{l_p}$ is spin damping derivative. $C_{m_{\alpha }}$ is pitching moment with respect to angle of attack derivative, $C_{n_{\beta }}$ is yawing moment derivative with respect to sideslip angle. $C_{m_q}$ is pitching moment coefficient derivative with pitch rate and $C_{n_r}$ is yawing moment coefficient derivative with yaw rate. ${\delta }_e$ is the parameter that describes the main motor state (${\delta }_e=0$ for active phase of flight and ${\delta }_e=1$ after main motor burnout, for gliding flight). When the main motor operates the projectile base drag is lower than after main motor burnout. $C_{X_0}$ was obtained for two system configurations (main motor on/off) and ${\delta }_e$ is used in a simulation to switch between aerodynamic data tables. Aerodynamic coefficients were obtained using commercially available software PRODAS (Projectile Rocket Ordnance Design & Analysis System). These coefficients were implemented into the Simulink model using Lookup-Table methodology. The aerodynamic coefficients are presented in Figure 2.

2.1.7. Gravity

Gravitational acceleration vector in the gravitational coordinate system is given as $\mathbf{g}={\left[\begin{array}{ccc}0& 0& g_0\end{array}\right]}^T$. Gravitational acceleration is assumed constant and consistent with WGS-84 reference model [25], i.e., $g_0=9.80665$ m/s${}^2$. Gravitational force and torques are given as:

(34)${\mathbf{F}}_{\mathbf{g}}={\mathbf{T}}_{\mathbf{b}}^{\mathbf{g}}m\left[\begin{array}{c}0\\ 0\\ g_0\end{array}\right]$

(35)${\mathbf{M}}_{\mathbf{g}}={\mathbf{r}}_{\mathbf{C}}\times {\mathbf{F}}_{\mathbf{g}}$

where $T_b^g$ is the transformation matrix from gravitational to body coordinate system, given as:

(36)${\mathbf{T}}_{\mathbf{g}}^{\mathbf{b}}=\left[\begin{array}{ccc}\cos \mathsf{\Theta }\cos \mathsf{\Psi }& \cos \mathsf{\Theta }\sin \mathsf{\Psi }& -\sin \mathsf{\Theta }\\ \sin \mathsf{\Phi }\sin \mathsf{\Theta }\cos \mathsf{\Psi }-\cos \mathsf{\Phi }\sin \mathsf{\Psi }& \sin \mathsf{\Phi }\sin \mathsf{\Theta }\sin \mathsf{\Psi }+\cos \mathsf{\Phi }\cos \mathsf{\Psi }& \sin \mathsf{\Phi }\cos \mathsf{\Theta }\\ \cos \mathsf{\Phi }\sin \mathsf{\Theta }\cos \mathsf{\Psi }+\sin \mathsf{\Phi }\sin \mathsf{\Psi }& \cos \mathsf{\Phi }\sin \mathsf{\Theta }\sin \mathsf{\Psi }-\sin \mathsf{\Phi }\cos \mathsf{\Psi }& \cos \mathsf{\Phi }\cos \mathsf{\Theta }\end{array}\right]$

2.1.8. Thrust

Thrust vector, which can deviate from the rocket’s longitudinal axis by angle ${\mathsf{\Theta }}_T$ in pitch plane and ${\mathsf{\Psi }}_T$ in the yaw plane, is given as [16]:

(37)${\mathbf{F}}_{\mathbf{s}}=F_p\left(t\right)\left[\begin{array}{c}\cos {\mathsf{\Theta }}_T\cos {\mathsf{\Psi }}_T\\ \cos {\mathsf{\Theta }}_T\sin {\mathsf{\Psi }}_T\\ -\sin {\mathsf{\Psi }}_T\end{array}\right]$

where $F_p\left(t\right)$ is the instantaneous value of the thrust force. Torque from the thrust force with respect to $O_b$ is given as:

(38)${\mathbf{M}}_{\mathbf{s}}=(-{\mathbf{r}}_{\mathbf{wC}}+{\mathbf{r}}_{\mathbf{C}})\times {\mathbf{F}}_{\mathbf{s}}$

2.1.9. Correction Thrusters

For the gasodynamic control system comprised of a set of identical correction thrusters placed radially in a set of parallel layers, the thrust and torque from the thrusters are given as:

(39)${\mathbf{F}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{j}}}=F_{p_{sk}}\left(t\right)\left[\begin{array}{c}0\\ \sin {\mathsf{\Phi }}_{i,j}\\ -\cos {\mathsf{\Phi }}_{i,j}\end{array}\right]$

(40)${\mathbf{M}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{j}}}=\left({\mathbf{r}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{j}}}-{\mathbf{r}}_{\mathbf{wC}}+{\mathbf{r}}_{\mathbf{C}}\right)\times {\mathbf{F}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{f}}}$

where $F_{p_{sk}}\left(t\right)$ is the instantaneous thrust force of the thruster, index $i=1,\cdots ,M$ is the layer number, index $j=1,\cdots ,N$ is the number of a thruster in a particular layer, ${\mathsf{\Phi }}_{i,j}$ is the azimuth angle of a thruster is a particular layer, ${\mathbf{r}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{j}}}$ is the vector describing the position of the layer with respect to the rocket’s base, measured from the base in the direction of the rocket’s axis. The total force and torque generated by the gasodynamic control system is given as: (41) and (42)

(41)${\mathbf{F}}_{\mathbf{sk}}=\sum _{i=1}^M\sum _{j=1}^N{\mathbf{F}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{j}}}$

(42)${\mathbf{M}}_{\mathbf{sk}}=\sum _{i=1}^M\sum _{j=1}^N{\mathbf{M}}_{{\mathbf{sk}}_{\mathbf{i},\mathbf{j}}}$

For the simulation purposes it was assumed that the projectile is equipped in a modular unit (Figure 3) composed from 32 solid propellant lateral thrusters and placed before the center of mass of the missile. These thrusters are set into a 4 arrays with 8 motors in each layer. Each of the thrusters might by used only once. The mass of the propellant in the single motor is approximately 0.005 kg so it is reasonable to assume that the ignition does not influence the mass and inertia projectile properties. The aerodynamic interference effects of the thrusters with the external flow were also omitted.

2.1.10. Inertial Parameters

The instantaneous mass of the rocket is given as:

(43)$m\left(t\right)=m_0-\frac{m_p}{I_c}{\int }_{t_0}^t{F}_p\left(t\right)dt$

where $m_0$ is the starting mass of the rocket at time $t_0$, $m_p$ is the mass of the propellant and $I_c$ is the total impulse given as:

(44)$I_c={\int }_{t_0}^{t_k}{F}_p\left(t\right)dt$

where $t_k$ is the time of propellant burnout. During the powered flight, the rocket’s mass center position vector ${\mathbf{r}}_{\mathbf{wC}}$ measured from the rocket’s base is changing according to:

(45)${\mathbf{r}}_{\mathbf{wC}}=\left[\begin{array}{ccc}{x}_{cg}\left(t\right)=x_{c{g}_0}-\frac{x_{c{g}_0}-x_{c{g}_k}}{I_c}{\int }_{t_0}^t{f}_p\left(t\right)dt& 0& 0\end{array}\right]$

where $x_{c{g}_0}$ is the center of mass position on the $O_b{x}_b$ axis during launch and $x_{c{g}_k}$ is the center of mass position on the $O_b{x}_b$ axis after the propellant burnout. The change of moments of inertia can be express as:

(46)$I_{ij}\left(t\right)=I_{i{j}_0}-\frac{I_{i{j}_0}-I_{i{j}_k}}{I_c}{\int }_{t_0}^t{F}_p\left(t\right)dt$

where $I_{i{j}_0}$ is the moment of inertia tensor component during launch and $I_{i{j}_k}$ is the moment of inertia tensor component after the propellant burnout.

2.1.11. Atmosphere Model

The air density, temperature, and the speed of sound are calculated according to the International Standard Atmosphere model [24].

(47)$\rho ={\rho }_0{\left(1-\frac{h}{44300}\right)}^{4.256}$

(48)$T=T_0-0.0065h$

(49)$a=a_0\sqrt{\frac{T}{288}}$

The reference values of these thermodynamic parameters are taken for the troposphere: ${\rho }_0=1.225$ kg/m${}^3$, $T_0=288.15$ K, $a_0=340.3$ m/s and $h=-z_n$ is the height in meters. It was assumed that the flight take place in the steady state atmosphere (wind speed was set to 0 m/s).

2.1.12. Inertial Measurement Unit Model

It was assumed that the rocket’s is equipped with the strapdown Inertial Measurement Unit with the three-axis accelerometer and three-axis gyroscope, and these are the only sources of information about rocket position, velocity and orientation. Abovementioned design requirements are quite difficult to fullfill due to errors. Pure inertial navigation has a tendency due to drift. These drift errors might be reduced using integration with GPS receivers. Also, additional sensors like might be used to improve the system acccuracy. Magnetometers measurement are imprecise because the projectile and launcher structure are made from steel alloys. The IMU is intended for a projectile that spins about the longitudinal axis of symmetry. Photodiode sensors could be used to measure the projectile roll rate.

An example of IMU that is suitable for the considered application is Micro-Electro-Mechanical-Systems based HG1930. The mass of this device is approximately 0.16 kg and power consumption less than 3 W. The operating temperature range is from $-54{}^{°}$C up to $+85{}^{°}$C. The gyroscopes ranges are up to 7200 deg/s in the X axis and 1440 deg/s in the Y and Z axes (X axis range must be significanlty higher than Y and Z due to projectile axial spin). Accelerometers operating range is up to 85 g in the X axes and 35 g in the Y and Z axes. High measurement range for X axis results from acceleration caused by main motor. This device requires supply voltage of 5 V. The IMU might be connected with the central onboard computer using military standard RS-422 serial interface. The maximum rate of data transmission for control purposes is 600 Hz. This measurement device is placed in front of the missile center of mass (between main motor unit and lateral thrusters module).

Accelerometers model

The acceleration of the rocket’s center of mass in $O_b{x}_b{y}_b{z}_b$ coordinate system is:

(50)$\mathbf{a}=\left[\begin{array}{c}{a}_x\\ a_y\\ a_z\end{array}\right]=\frac{{\mathbf{F}}_{\mathbf{b}}}{m}$

In the general case, the accelerometer position need not to coincide with the center of mass of the rocket. Therefore, the center of mass acceleration must be recalculated to the point of accelerometers’ mounting:

(51)${\mathbf{a}}_{\mathbf{IMU}}=\mathbf{a}+\mathsf{\Omega }\times \left(\mathsf{\Omega }\times {\mathbf{r}}_{\mathbf{wz}}\right)+\dot{\mathsf{\Omega }}\times {\mathbf{r}}_{\mathbf{wz}}-\mathbf{g}$

where ${\mathbf{r}}_{\mathbf{wz}}={\mathbf{r}}_{\mathbf{wC}}-{\mathbf{r}}_{\mathbf{IMU}}$ is the position of the IMU with respect to center of mass and ${\mathbf{r}}_{\mathbf{IMU}}$ is the IMU position with respect to the rocket’s base. As a next step, the model of the sensor’s errors was included, which is comprised of scale factors $s_x,s_y,s_z$, cross-coupling $c_{xy},c_{xz},c_{yz}$ and biases $b_x,b_y,b_z$:

(52)${\widehat{\mathbf{a}}}_{\mathbf{IMU}}=\left[\begin{array}{ccc}{s}_x& -c_{xy}& c_{xz}\\ c_{xy}& s_y& -c_{yz}\\ -c_{xz}& c_{yz}& s_z\end{array}\right]{\mathbf{a}}_{\mathbf{IMU}}+\left[\begin{array}{c}{b}_x\\ b_y\\ b_z\end{array}\right]$

Accelerometer is treated as a second-order dynamic system:

(53)${\mathbf{a}}_{\mathbf{meas}}=\frac{{\omega }_{n_{acc}}^2}{s^2+2{\xi }_{acc}{\omega }_{n_{acc}}s+{\omega }_{n_{acc}}^2}{\widehat{\mathbf{a}}}_{\mathbf{IMU}}$

where ${\xi }_{acc}$ is the accelerometer damping coefficient and ${\omega }_{n_{acc}}$ is the accelerometer natural frequency (${\xi }_{acc}=0.707$ and ${\omega }_{n_{acc}}=7600$). The last step was to include the sensor noise, assumed as white noise with known standard deviation and zero mean, and output saturation.

Gyroscopes Model

The gyroscope output does not depend on the gyroscopes’ position inside the rocket. Therefore, there is no need to transform its output to the center of mass. The gyroscope errors’ model includes scale factor, cross-coupling, bias, and sensitivity to accelerations, given by the gyroscopes’ sensitivity matrix $\mathbf{G}$:

(54)${\mathsf{\Omega }}_{\mathbf{IMU}}=\left[\begin{array}{ccc}{s}_x& -c_{xy}& c_{xz}\\ c_{xy}& s_y& -c_{yz}\\ -c_{xz}& c_{yz}& s_z\end{array}\right]\mathsf{\Omega }+\left[\begin{array}{c}{b}_x\\ b_y\\ b_z\end{array}\right]+\mathbf{G}\left[\begin{array}{c}{a}_x\\ a_y\\ a_z\end{array}\right]$

Gyroscope is also treated as a second-order system:

(55)${\mathsf{\Omega }}_{\mathbf{meas}}=\frac{{\omega }_{n_{gyro}}^2}{s^2+2{\xi }_{gyro}{\omega }_{n_{gyro}}s+{\omega }_{n_{gyro}}^2}{\mathsf{\Omega }}_{\mathbf{IMU}}$

where ${\xi }_{gyro}$ is the gyroscope damping coefficient and ${\omega }_{n_{gyro}}$ is the gyroscope natural frequency (it was assumed that ${\xi }_{gyro}=0.356$ and ${\omega }_{n_{gyro}}=7600$ Hz). The last step, as with accelerometers, was to include the sensor noise and output saturation. The values of ${\mathbf{a}}_{\mathbf{meas}}$ and ${\mathsf{\Omega }}_{\mathit{m}\mathit{e}\mathit{a}\mathit{s}}$ are sampled with the sensor sample frequency. From the measured angular velocity, the rocket’s orientation is calculated. To obtain the rocket’s velocity, the measured accelerations are firstly recalculated back to the center of mass position:

(56)${\mathbf{a}}_{\mathbf{meas},\mathbf{CG}}={\mathbf{a}}_{\mathbf{meas}}+{\mathsf{\Omega }}_{\mathit{m}\mathit{e}\mathit{a}\mathit{s}}\times \left({\mathsf{\Omega }}_{\mathit{m}\mathit{e}\mathit{a}\mathit{s}}\times {\mathbf{r}}_{\mathbf{wz}}\right)+{\dot{\mathsf{\Omega }}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{s}}\times {\mathbf{r}}_{\mathbf{wz}}$

and then the velocity vector can be obtained by numerical integration:

(57)${\dot{\mathbf{V}}}_{\mathbf{b},\mathbf{meas}}={\mathbf{a}}_{\mathbf{meas},\mathbf{CG}}-{\mathsf{\Omega }}_{\mathit{m}\mathit{e}\mathit{a}\mathit{s}}\times {\mathbf{V}}_{\mathbf{b},\mathbf{meas}}+\Lambda {\left[\begin{array}{ccc}0& 0& g_0\end{array}\right]}^T$

where ${\mathbf{V}}_{\mathbf{b},\mathbf{meas}}={[U_{meas},V_{meas},W_{meas}]}^T$. Velocity vector is expressed in body frame $O_b{x}_b{y}_b{z}_b$ so it must be transformed to the navigational coordinate system $O_n{x}_n{y}_n{z}_n$. The rocket’s position in the navigational coordinate system is calculated by numerical integration of velocity components:

(58)$\left[\begin{array}{c}{\dot{x}}_{n,meas}\\ {\dot{y}}_{n,meas}\\ {\dot{z}}_{n,meas}\end{array}\right]=\mathsf{\Lambda }\left[\begin{array}{c}{U}_{meas}\\ V_{meas}\\ W_{meas}\end{array}\right]$

In order to solve the mentioned earlier navigation equations the initial attitude, velocity and position of the projectile must be known. It was assumed that these parameters are perfectly estimated using initial alignment procedure before the flight.

2.1.13. Control Algorithms

Projectiles equipped in solid propellant lateral motors often have very low control authority. It means that the maneuverability of such objects is small. This fact makes the guidance process challenging. For the tests of landing dispersion analysis, two algorithms were created: Multi-Condition Control Algorithm (MCCA) and modified Proportional Navigation Guidance (mPNG). The guidance algorithm used in the MCCA is based on the reference trajectory tracking. This algorithm is discussed in detail in [33]. The main idea behind MCCA is to use the reference trajectory to minimize the hitting error just and the end of flight. Due to limited number of lateral motors it is difficult to track the trajectory along the full flight path. In MCCA approach the position error between reference trajectory and the actual projectile position is minimized after trajectory vertex, during the descending flight. It is assumed that the reference trajectory is calculated prior to launch, so the position of the target is known a priori, and implemented in the rocket’s control system prior to launch. The reference trajectory is calculated for unguided projectile in such a way the missile hits perfectly the target (miss distance at the end of nominal trajectory is 0 m). The guidance algorithm used in the mPNG algorithm is the classical Proportional Navigation Guidance [34] modified by the term accounting for the trajectory bending due to gravity. Both algorithms use the same thrusters’ ignition logic that is presented in [12,16]. Due to the rocket’s high roll angular velocity during the flight, the correction thrusters must be ignited in the right moment, which means when the rocket achieves a certain roll angle. At any moment, only one correction thruster can be ignited. The set of conditions of thrusters’ ignition, common for both algorithms are:

• Correction thruster was not used already (solid motor thrusters are single-use motors)

• The time between the last thruster ignition $t_{last}$ is greater than some limit value $\tau \in (0;\infty )$

  (59)$t-t_{last}>\tau$

• The correction thruster must be ignited so that the resultant thrust force was in the direction of the desired lateral displacement [35,36], which means that the absolute value of the difference between the error phase $\gamma$ and thruster azimuth angle ${\mathsf{\Phi }}_{i,j}$ diminished by the control prediction times ${\tau }_d$ and ${\tau }_{sk}$ multiplied by the roll angular velocity was lower than some limit value ${\gamma }_t$.

  (60)$\left|\gamma -{\mathsf{\Phi }}_{i,j}-\pi -P({\tau }_d+{\tau }_{sk})\right|\le {\gamma }_t$

• The rocket’s pitch angle must be lower or equal the threshold value ${\mathsf{\Theta }}_g$ and the time of flight must be at least equal to the threshold value $t_g$

  (61)$\mathsf{\Theta }\le {\mathsf{\Theta }}_g\wedge t\ge t_g$

Additional conditions for the MCCA algorithm:

• The distance between the rocket’s center of mass and the reference trajectory $\mathsf{\Gamma }$, measured perpendicular, is greater than some limit value ${\mathsf{\Gamma }}_t$

  (62)$\mathsf{\Gamma }>{\mathsf{\Gamma }}_t$

Additional conditions for the mPNG algorithm:

• the norm of the commanded value of the lateral acceleration ${\mathbf{a}}_{\mathbf{cmd}}$ must be greater than the threshold value ${\mathbf{a}}_{\mathbf{cmd},\mathbf{g}}$

  (63)${\mathbf{a}}_{\mathbf{cmd}}\ge {\mathbf{a}}_{\mathbf{cmd},\mathbf{g}}$

The parameters of the control laws were determined using the expert method and parametric study: $\tau =0.2$ s, ${\tau }_d=0.001$ s, ${\tau }_{sk}=0.015$ s, ${\gamma }_t=2.5$ deg, ${\mathsf{\Theta }}_g=-10$ deg (the guidance process starts after trajectory vertex), $t_g=15$ s, ${\mathsf{\Gamma }}_t=1$ m, $a_{cmd,g}=3$ m/s${}^2$.

2.2. Simulation Model

The mathematical model described in Section 2.1 was implemented in MATLAB/ Simulink 2020b environment. The main Simulink block model of the system is presented in Figure 4.

The program simulates the flight of the gasodynamically controlled rocket, calculates the loads from gravity, aerodynamics, thrust, and correction thrusters. It solves the set of ordinary differential equations for the rigid body with 6 degrees of freedom and variable mass. It includes the models of International Standard Atmosphere and Inertial Measurement Unit as well as the inertial navigation equation for determining the rocket’s position, orientation, and velocity. The equations of motion of the projectile were integrated using fixed step, third order Bogacki-Shampine method. The step size was set to 0.0001 s. Simulations were realized using Simulink build in option “Accelerator mode”. Marsenne-Twister algorithm [37] was used to generate in a pseudorandom way the disturbances for the Monte-Carlo simulation. The model was optimized to make the run time as short as possible. The simulation might be realized in a batch mode from the external MATLAB script.

3. Results and Discussion

3.1. Input Data for the Simulation Study

A generic rocket model was used for the simulations, which general data are provided in Table 1.

The missile is stabilized with four trapezoidal fins. The maximum flight velocity of this projectile is 605 m/s and maximum roll rate 4700 deg/s (these values are obtained in 3 s of flight). The missile was fired at elevation angle 25 deg. Initial velocity was set to 42 m/s and initial roll rate 1073 deg/s. Projectile was fired from the initial position (0, 0, 0) m. The reference trajectory of the projectile (unguided flight) is presented in Figure 5.

3.2. Initial Verification of Control Algorithms

As a first step, five deterministic cases were evaluated to test if the prepared control algorithms work as intended: no control, mPNG algorithm with and without IMU model, and MCCA algorithm with and without IMU model. To intentionally introduce aiming error it was assumed that the launch tube is not perfectly aligned with the demanded shoot direction. The initial heading error angle equaled 2 degrees. The position coordinates of the stationary target were set to $(9296.54,-7.29,0)$ m. Table 2 presents the results of the performed cases. The first column describes the used algorithm, with the information whether the IMU model was on. The next columns present the error between the x and y components of the rocket’s and reference trajectories at the impact point, and the distance between the rocket’s landing point and the target position given as $\Delta \mathbf{R}=\sqrt{\Delta {\mathbf{X}}^2+\Delta {\mathbf{Y}}^2}$.

The results shown in Table 2 indicate that the largest error was obtained for the uncontrolled flight. The projectile landed 324.5 m from the desired point and it is obvious that the target was not achieved (typical radius of destruction for rocket artillery projectiles is order of 20–40 m). Both algorithms work properly and that the mPNG control algorithm lowers the miss distance by about $92\%$ and MCCA algorithm by about $75\%$ on average. The IMU model errors slightly increase the miss distance for both algorithms.

3.3. Monte-Carlo Simulations

Next, to test the influence of various uncertainties on the performance of control algorithms and the resulting landing dispersion, a few sets of simulations were performed. Tested were uncertainties in aerodynamic data, thrust parameters, and initial conditions. For every uncertainty, again five cases were simulated: no control, mPNG algorithm with and without IMU model, and MCCA algorithm with and without IMU model. Every case took 1000 runs, using the Monte Carlo method, giving a total of 15 thousand runs. As a merit of accuracy, the Circular Error Probable CEP was used. It gives information about the radius of a circle inside which 50% of landing points are located.

3.3.1. Aerodynamic Parameters Uncertainties

The first set of Monte Carlo simulations consisted of uncertainties in aerodynamic data. It was assumed that the normal distribution standard deviation of all aerodynamic parameters was equal to $\sigma =0.2$. Next, the maps of impact points were obtained. The Figure 6 presents the results of the performed simulations. On the horizontal axis there is crossrange and on the vertical axis range of the projectile.

In the not controlled case, the CEP was equal to $231.14$ m. Ideal case of mPNG and MCCA algorithms (IMU model off) achieved $85.6\%$ and $83.1\%$ miss distance reduction respectively and with the IMU model on $84.5\%$ for mPNG and $82.4\%$ for MCCA. Much bigger dispersion is observed along the rocket’s flight path, because the drag coefficient uncertainties have the biggest influence on the range. Small directional dispersion is mostly caused by the IMU model errors. The achieved miss distance reduction for both algorithms were very similar.

3.3.2. Thrust Curve Uncertainties

The next set of simulations consisted of uncertainties in the thrust data. For simplicity, it was assumed that the thrust curve can be approximated by the quadrilateral comprised of four characteristic pairs of points, time-thrust, presented in Figure 7. Every point was randomly chosen, using a uniform distribution, between the maximum and minimum allowable values, with additional constraints that the 4th time had to be larger that the 3rd time and that the total impulse of the thrust should lie in between allowable values. The values of the uncertainties are presented in Table 3. In this way a set of pseudorandom thrust curves was obtained as input data for the Monte-Carlo simulations (in each simulation run a different thrust curve was used).

The Figure 8 presents the results of the performed simulations. Again, the largest dispersion of the impact points was observed for the uncontrolled projectile. The uncontrolled flight case resulted in CEP of $66.90$ m. It means, that thrust uncertainties produce smaller dispersion than uncertainties in aerodynamic parameters. This dispersion is reduced significantly in controlled shoots. Results for mPNG and MCCA algorithms without IMU achieved $92.1\%$ and $89.2\%$ miss distance reduction respectively and with the IMU model on $89.2\%$ for mPNG and $86.7\%$ for MCCA. Again, a bigger dispersion is observed in the longitudinal direction, which thrust uncertainties affect the most. Directional dispersion is a bit lower for the mPNG algorithm. Again, the results for both algorithms lied very close.

3.3.3. Initial Condition Uncertainties

The last set of simulations consisted of uncertainties in the initial conditions. The values of the initial linear velocity vector, angular velocity vector, and orientation angles were chosen randomly, using a normal distribution, using set values of standard deviation. The Figure 9 presents the results of the performed simulations. The uncontrolled case resulted in CEP of $67.09$ m and the longitudinal and directional dispersion were the same. This is a typical dispersion pattern that is obtained for rocket artillery projectiles at medium elevation angles. Results for mPNG and MCCA algorithms without IMU achieved $95.2\%$ and $80.3\%$ miss distance reduction respectively and with the IMU model on $93.7\%$ for mPNG and $79.1\%$ for MCCA. In this scenario, the miss distance reduction is in favour of the mPNG algorithm with around $10\%$ difference in results. Directional dispersion is more affected by control than the longitudinal, which may arise from the difference in longitudinal and lateral velocity of the rocket.

Several new aspects brought by the paper might be mentioned. First, the influence of measurement errors on the resulting projectile miss distance was investigated. Second, two different guidance methods intended for lateral thrusters controlled missiles were compared for idealized and realistic case. From the obtained results it might be concluded that it is possible to achieve CEP order of several meters but to realize this goal the missile must be equipped in high-accuracy IMU.

4. Conclusions

Precision guided munition become more and more important in modern military conflicts. To achieve a high direct hit probability the influence of various factors on the resulting dispersion must be understood in detail. In the article, the impact point dispersion, caused by the uncertainties in aerodynamic parameters, thrust curve, initial conditions, and on-board measuring devices, for two types of control algorithms, was presented. Monte Carlo approach was used in the simulations and as the merit of dispersion the CEP was utilized.

The results showed that the use of control algorithms greatly reduces the miss distance by more than 80% in most simulated cases. From the simulations it might be concluded that there is possible to achieve CEP smaller than 8 m. For modern guided munition this is quite a realistic result (for example, 160 mm ACCULAR projectile has declared CEP < 10 m). It means that equipping the projectile with control module composed from lateral thrusters allows effectively reduce the impact points dispersion. The uncontrolled projectile might land even 200 m from the intended location in the worst-case scenario. This issue is very important in modern military applications due to the requirement of minimizing the collateral damage.

The IMU model causes a slight increase in the dispersion of about 3% in every case. The mPNG algorithm proved to be better for all simulated cases, the greatest difference between the two algorithms was observed in the initial condition case dispersion.

The developed numerical simulation might be used in the design of new measurement systems intended for missile navigation. Parametric model allows on rapid implementation of data for other missiles and IMU-s and investigate the dispersion as a function of measurement uncertainties. In this way the overall time and cost of the projectile design might be reduced.

Further works might concentrate on flight tests of the real ground-to-ground projectile and validating the model. Also wind tunnel measurements of the missile could be evaluated to obtain the aerodynamic data for a wide range of flight conditions. The influence of wind on the projectile dispersion might be also explored in detail. Hardware-in-the-loop simulation might be also considered to investigate the influence of sensor errors on the projectile hitting accuracy.

Author Contributions

Conceptualization, R.G. and P.L.; methodology, M.J. and P.L.; software, M.J. and D.M.; validation, M.J. and P.L.; formal analysis, R.G.; investigation, D.M.; resources, D.M. and M.J. and P.L.; data curation, D.M.; writing—original draft preparation, D.M., R.G. and P.L.; writing—review and editing, P.L.; visualization, D.M.; supervision, R.G. and P.L.; project administration, P.L.; funding acquisition, P.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Centre for Reaserch and Development grant number DOB-SZAFIR/03/B/002/01/2021.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Definition of measuring and aerodynamic coordinate systems, and navigational, gravitational and body coordinate systems.

Figure 2. Aerodynamic coefficients.

Figure 3. Configuration of the lateral thrusters module.

Figure 4. Top level architecture of the Simulink simulational model.

Figure 5. Projectile trajectory.

Figure 6. Results of landing point dispersion with uncertainties of aerodynamic data: (a) no control, (b) mPNG, (c) mPNG + IMU, (d) MCCA, (e) MCCA + IMU.

Figure 6. Results of landing point dispersion with uncertainties of aerodynamic data: (a) no control, (b) mPNG, (c) mPNG + IMU, (d) MCCA, (e) MCCA + IMU.

Figure 7. Thrust curve.

Figure 8. Results of landing point dispersion with uncertainties of thrust data: (a) no control, (b) mPNG, (c) mPNG + IMU, (d) MCCA, (e) MCCA + IMU.

Figure 9. Results of landing point dispersion with uncertainties of initial conditions: (a) no control, (b) mPNG, (c) mPNG + IMU, (d) MCCA, (e) MCCA + IMU.

Figure 9. Results of landing point dispersion with uncertainties of initial conditions: (a) no control, (b) mPNG, (c) mPNG + IMU, (d) MCCA, (e) MCCA + IMU.

Table 1

Generic rocket’s parameters.

Parameter	Value	Unit
diameter	122	mm
length	$1.58$	m
initial mass	$22.14$	kg
propellant mass	$5.83$	kg
initial moment of inertia $I_{xx}$	$0.0422$	kg·m${}^2$
final moment of inertia $I_{xx}$	$0.0326$	kg·m${}^2$
initial moment of inertia $I_{yy}$	$11.223$	kg·m${}^2$
final moment of inertia $I_{yy}$	$9.513$	kg·m${}^2$
initial moment of inertia $I_{zz}$	$11.223$	kg·m${}^2$
final moment of inertia $I_{zz}$	$9.513$	kg·m${}^2$
maximum thrust	$7277.5$	N
average thrust	$3383.2$	N
burn time	$3.31$	s
total impulse	13,529	N
correction thruster’s thrust	200	N
correction thruster’s burn time	$0.03$	s
number of correction thrusters per layer	8	-
number of correction thrusters’ layers	4	-

Table 2

The errors in landing point components.

Algorithm	$\Delta \mathbf{X}$[m]	$\Delta \mathbf{Y}$[m]	$\Delta \mathbf{R}$[m]
None	$-5.4$	$324.4$	$324.5$
mPNG	$-12.3$	$19.7$	$23.2$
mPNG + IMU	$-15.3$	$20.1$	$25.3$
MCCA	$-32.1$	$69.8$	$76.8$
MCCA + IMU	$-32.0$	$75.2$	$81.8$

Table 3

Thrust curve uncertainties.

Point	Value	Time t [s]	Thrust T [N]
1	min	0	0
1	max	0	0
2	min	0	3600
2	max	$0.15$	4400
3	min	$1.9$	6900
3	max	$2.5$	8700
4	min	$2.1$	0
4	max	$2.7$	0

References

1. Liu, L.J.; Shen, Y. Three-Dimension H_∞ Guidance Law and Capture Region Analysis. Trans. Aerosp. Electron. Syst.; 2012; 48, pp. 419-429.

2. Budiyono, A.; Rachman, H. Proportional Guidance and CDM Control Sythesis for a Short-Range Homing Surface-to-Air Missile. J. Aerosp. Eng.; 2012; 25, pp. 168-177. [DOI: https://dx.doi.org/10.1061/(ASCE)AS.1943-5525.0000104]

3. Guo, Y.; Li, X.; Zhang, H.; Cai, M.; He, F. Data-Driven Method for Impact Time Control Based on Proportional Navigation Guidance. J. Guid. Control. Dyn.; 2020; 43, pp. 955-966. [DOI: https://dx.doi.org/10.2514/1.G004669]

4. Chen, X.; Wang, J. Optimal Control Based Guidance Law to Control Both Impact Time and Impact Angle. Aerosp. Sci. Technol.; 2018; 84, pp. 454-463. [DOI: https://dx.doi.org/10.1016/j.ast.2018.10.036]

5. He, S.; Lee, C.H. Gravity-Turn-Assisted Optimal Guidance Law. J. Guid. Control. Dyn.; 2018; 41, pp. 171-183. [DOI: https://dx.doi.org/10.2514/1.G002949]

6. He, S.; Shin, H.S.; Tsourdos, A.; Lee, C.H. Energy-Optimal Waypoint-Following Guidance Considering Autopilot Dynamics. Trans. Aerosp. Electron. Syst.; 2020; 56, pp. 2701-2717. [DOI: https://dx.doi.org/10.1109/TAES.2019.2954149]

7. Guelman, M.; Shinar, J. Optimal Guidance Law in the Plane. J. Guid.; 1984; 7, pp. 471-476. [DOI: https://dx.doi.org/10.2514/3.19880]

8. Burchett, B.; Costello, M. Model Predictive Lateral Pulse Jet Control of an Atmospheric Rocket. J. Guid. Control. Dyn.; 2002; 25, pp. 860-867. [DOI: https://dx.doi.org/10.2514/2.4979]

9. Gupta, S.K.; Saxena, S.; Singhal, A.; Ghosh, A.K. Trajectory Correction Flight Control System using Pulsejet on an Artillery Rocket. Def. Sci. J.; 2008; 58, pp. 15-33. [DOI: https://dx.doi.org/10.14429/dsj.58.1621]

10. Głębocki, R. Guidance Impulse Algorithms for Air Bomb Control. Bull. Pol. Acad. Sci.; 2012; 60, pp. 825-833. [DOI: https://dx.doi.org/10.2478/v10175-012-0096-4]

11. Gao, M.; Zhang, Y.; Yang, S. Firing Control Optimization of Impulse Thrusters for Trajectory Correction Projectiles. Int. J. Aerosp. Eng.; 2015; 2015, [DOI: https://dx.doi.org/10.1155/2015/781472]

12. Jitpraphai, T.; Costello, M. Dispersion Reduction of a Direct Fire Rocket Using Lateral Pulse Jets. J. Spacecr. Rocket.; 2001; 38, pp. 929-936. [DOI: https://dx.doi.org/10.2514/2.3765]

13. Trzun, Z.; Vrdoljak, M.; Cajner, H. Statistical Analysis of Missile’s Trajectory Deviation Due to Production Errors. Proc. Appl. Math. Mech.; 2020; 20, e202000059. [DOI: https://dx.doi.org/10.1002/pamm.202000059]

14. Trzun, Z.; Vrdoljak, M. Monte Carlo Simulation of Missile Trajectories Dispersion Due to Imperfectly Manufactured Warhead. Proceedings of the 31st DAAAM International Symposium; Vienna, Austria, 21–24 October 2020; Katalinic, B. DAAAM International: Vienna, Austria, 2020.

15. Mingireanu, F. Trajectory Modelling of GRAD Rocket with Low-Cost Terminal Guidance Upgrade Coupled to Range Increase Through Step-Like Thrust-Curves. Proceedings of the Aerospace Sciences & Aviation Technology; Cairo, Egypt, 28–30 May 2013.

16. Głębocki, R.; Jacewicz, M. Parametric Study of Guidance of a 160-mm Projectile Steered with Lateral Thrusters. Aerospace; 2020; 7, 61. [DOI: https://dx.doi.org/10.3390/aerospace7050061]

17. Głębocki, R.; Jacewicz, M. Sensitivity Analysis and Flight Tests Results for a Vertical Cold Launch Missile System. Aerospace; 2020; 7, 168. [DOI: https://dx.doi.org/10.3390/aerospace7120168]

18. Figueroa-Gonzales, A.; Cacciatore, F.; Haya-Ramos, R. Landing Guidance Strategy of Space Rider. J. Spacecr. Rocket.; 2021; 54, pp. 1220-1231. [DOI: https://dx.doi.org/10.2514/1.A34957]

19. Mandić, S. Dispersion Reduction of Artillery Rocket Guided by Flight Path Steering Method. Aeronaut. J.; 2016; 120, pp. 435-456. [DOI: https://dx.doi.org/10.1017/aer.2016.4]

20. De Celis, R.; Cadarso, L.; Sánchez, J. Guidance and Control for High Dynamic Rotating Artillery Rockets. Aerosp. Sci. Technol.; 2017; 64, pp. 204-212. [DOI: https://dx.doi.org/10.1016/j.ast.2017.01.026]

21. Chen, S.Y.; Xia, Q.L. A Multicontrained Ascent Guidance Method for Solid Rocket-Powered Launch Vehicles. Int. J. Aerosp. Eng.; 2016; 2016, 6346742. [DOI: https://dx.doi.org/10.1155/2016/6346742]

22. Saghafi, F.; Khalilidelshad, M. A Monte Carlo Dispersion Analysis of a Rocket Flight Simulation Software. Proceedings of the 17th European Simulation Multiconference; Nottingham, UK, 9–11 June 2003.

23. Le, V.D.T.; Nguyen, A.T.; Nguyen, L.H.; Dang, N.T.; Tran, N.D. Effectiveness Analysis of Spin Motion in Reducing Dispersion of Sounding Rocket Flight due to Thrust Misalignment. Int. J. Aeronaut. Space Sci.; 2021; 22, pp. 1194-1208. [DOI: https://dx.doi.org/10.1007/s42405-021-00383-x]

24. National Oceanic and Atmospheric Administration. U.S. Standard Atmosphere; Technical Report; Government Printing Office: Washington, DC, USA, 1976.

25. National Imagery and Mapping Agency. Department of Defense World Geodetic System 1984: Its Definition, and Relationship with Local Geodetic Systems, TR8350.2; Technical Report; Department of Defense: Washington, DC, USA, 2000.

26. Nizioł, J.; Radziszewski, B. Dynamika ukłAdów Mechanicznych; Komitet Mechaniki PAN: Warsaw, Poland, 2005.

27. Etkin, B. Dynamics of Atmospheric Flight; John Wiley & Sons: New York, NY, USA, 1972.

28. Rogers, R.M. Applied Mathematics in Integrated Navigation Systems; AIAA: Reston, VA, USA, 2000.

29. Zipfel, P.H. Modeling and Simulation of Aerospace Vehicle Dynamics; AIAA: Reston, VA, USA, 2007.

30. Allerton, D. Principles of Flight Simulation; John Wiley & Sons: Wiltshire, UK, 2009.

31. Stevens, B.L.; Lewis, F.L.; Johnson, E.N. Aircraft Control and Simulation. Dynamics, Controls Design, and Autonomous System; John Wiley & Sons: New York, NY, USA, 2016.

32. Winchenbach, G.L. Aerodynamic Testing in a Free-Flight Spark Range; Weapon Flight Mechanics Division: Valparaiso, FL, USA, 1997.

33. Ozog, R.; Głębocki, R.; Jacewicz, M. Modified Trajectory Tracking Guidance for Artillery Rocket. J. Theor. Appl. Mech.; 2020; 3, pp. 611-622. [DOI: https://dx.doi.org/10.15632/jtam-pl/121981]

34. Zarchan, P. Tactical and Strategic Missile Guidance; 6th ed. AIAA: Washington, DC, USA, 2012.

35. Xiaobing, C.; Yicen, X.; Changying, R.; Zhengya, S. Simulation of Trajectory Correction for an Impulse Control Mortar Projectile with a Strapdown Laser Seeker. Appl. Mech. Mater.; 2013; 367, pp. 377-381.

36. Lichota, P.; Jacewicz, M.; Szulczyk, J. Spinning Gasodynamic Projectile System Identification Experiment Design. Aircr. Eng. Aerosp. Technol.; 2020; 92, pp. 452-459. [DOI: https://dx.doi.org/10.1108/AEAT-06-2019-0124]

37. Matsumoto, M.; Nishimura, T. Mersenne twister: A 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Trans. Model. Comput. Simul.; 1998; 8, pp. 3-30. [DOI: https://dx.doi.org/10.1145/272991.272995]