# Simulation-Based Early Prediction of Rocket, Artillery, and Mortar Trajectories and Real-Time Optimization for Counter-RAM Systems

Authors: Arash Ramezani; Hendrik Rothe
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The threat imposed by terrorist attacks is a major hazard for military installations, for example, in Iraq and Afghanistan. The large amounts of rockets, artillery projectiles, and mortar grenades (RAM) that are available pose serious threats to military forces. An important task for international research and development is to protect military installations and implement an accurate early warning system against RAM threats on conventional computer systems in out-of-area field camps. This work presents a method for determining the trajectory, caliber, and type of a projectile based on the estimation of the ballistic coefficient. A simulation-based optimization process is presented that enables iterative adjustment of predicted trajectories in real time. Analytical and numerical methods are used to reduce computing time for out-of-area missions and low-end computer systems. A GUI is programmed to present the results. It allows for comparison between predicted and actual trajectories. Finally, different aspects and restrictions for measuring the quality of the results are discussed.


---

## Body

## 1. Introduction

Field camps are military facilities which provide living and working conditions in out-of-area missions. During an extended period of deployment abroad, they have to ensure safety and welfare for soldiers.Current missions in Iraq or Afghanistan have shown that the safety of military camps and air bases is not sufficient. A growing threat to these military facilities is the use of unguided rockets, artillery projectiles, and mortar grenades. Damage with serious consequences has occurred increasingly often in the past few years.This paper focuses on mortars and rockets because they are more and more used by irregular forces, where they have easy access to a large amount of these weapons. Further reasons are the small radar cross-section, the short firing distance, and the thick cases made of steel or cast-iron, which makes mortar projectiles and rockets hard to detect and destroy.The challenge is to establish an early warning system for different projectiles using analytical and numerical methods to reduce computing time and improve simulation results compared to similar systems. An appropriate estimation of the ballistic coefficient and the associated calculation of unknown parameters is the central issue in this field of research.Up to now, only a few approaches have been published. Khalil et al. [1] presented a trajectory prediction for the special field of fin stabilized artillery rockets. Chusilp et al. [2] compared 6-DOF trajectory simulations of a short range rocket using aerodynamic coefficients. A very good overview of modeling and simulation of aerospace vehicle dynamics is given by Zipfel [3].An et al. [4] used a fitting coefficient setting method to modify their point mass trajectory model. Chusilp and Charubhun [5] estimated the impact points of an artillery rocket fitted with a nonstandard fuze. Scheuermann et al. [6] characterized a microspoiler system for supersonic finned projectiles. Wang et al. [7] established a guidance and control design for a class of spin-stabilized projectiles with a two-dimensional trajectory correction fuze. Lee and Jun [8] developed guidance algorithm for projectile with rotating canards via predictor-corrector approach. Fresconi et al. [9] developed a practical assessment of real-time impact point estimators for smart weapons.This paper is based on Ramezani et al. [10]. Real-time prediction of trajectories and continuous optimization is one of the main aims of this work. With the aid of graphical solutions, it is possible to differentiate between several objects and determine firing locations as well as points of impact. The goal is to provide active protection of stationary assets in today’s crisis regions. Therefore, a modern counter-RAM system with a clear GUI must be developed and will then be employed for most threats.


## 2. Ballistic Model

The projectile is to be expected as a point mass: that is, the entire projectile mass is located in the center of gravity. Rotation is irrelevant in this case, so we regard a ballistic model with 3-DOF.The Earth can be regarded as a static sphere with infinite radius and represents an inertial system. Based on an Earth-fixed Cartesian coordinate system, the force of inertia is applied in a single direction.Different projectiles have to be considered in order to set up a mathematical model. While rockets can be regarded as spin-stabilized projectiles, which have a short phase of thrust and are particularly suitable for long distances up to 20 km, mortar grenades are arrow-stabilized and fired on short distances up to approximately 8 km.Other mathematical models for typical fin stabilized artillery rockets are presented in [11–16].


### 2.1. Exterior Ballistics

The ballistic model is principally based on Newton’s law and the equations of motion are considered to be under the effect of air drag and the force of gravity only. Additionally, rockets have a thrust vector impelling the projectile for a few seconds (generally, combustion gases have a velocity range of 1800–4500 m/s [18]). Anyhow, rockets as well as mortars have ballistic trajectories and the object is to identify the threat on the basis of different flight characteristics.Letg→ denote a reference acceleration (acceleration of gravity at sea level on Earth), with(1)g=9.80665m/s2, taking effect on the point mass in vertical direction.The air dragD→ can have different values, depending on the design of the projectile, that is,(i)
muzzle velocityv0,(ii)
weight,(iii)
aerodynamics, and the properties of air, for example, (i)
density,(ii)
temperature,(iii)
wind,(iv)
speed of sound.Considering the general formula(2)D→=12·Cd·A·ρ·v·v→,Cx=Cd·CA·Bcontaining all parameters named above with (i)
A: cross-section area of the projectile,(ii)
ρ: air density,(iii)
v: velocity of the projectile,(iv)
Cd: air drag coefficient,(v)
CA: environmental properties,(vi)
B: ballistic coefficient, it is operative to find an appropriate approximation, so that the projectile can be specified. The parameters A, ρ, Cd, CA, and B are unknown, whereas v can be defined precisely from the measured radar data.The air drag coefficientCd for instance depends on the critical velocity ratio, pictured in Figure 1. Since the drag coefficient does not vary in a simple manner with Mach number, this makes the analytic solutions inaccurate and difficult to accomplish.Figure 1
Characteristics of the air drag coefficientCd.One can see from this figure that there is no simple analytic solution to this variation. With computer power nowadays, we usually solve or approximate the exact solutions numerically, doing the quadratures by breaking the area under the curve into quadrilaterals and summing the areas. In general, there are three forms of the drag coefficient:(1)
ConstantCd that is useful for the subsonic flight regime: Ma<1(2)
Cd inversely proportional to the Mach number that is characteristic of the high supersonic flight regime: in this case, Ma≫1(3)
Cd inversely proportional to the square root of the Mach number that is useful in the low-supersonic flight regime: Ma≥1 Carlucci and Jacobson [19] give a detailed description of the air drag coefficient.Another coefficient in common use in ballistics is the ballistic coefficientB, which is defined as(3)B=md2, where m and d are the mass and diameter of the projectile [19]. Section 3.2 deals with the problem of estimating the unknown parameters.


### 2.2. Equations of Motion

The aerodynamics and ballistics literature are quite diverse and terminology is far from consistent. This has particular significance in the coordinate systems used to define the equations of motion. Nevertheless, this field of research has a long history and a lot of approaches. More details are discussed in [20–24].In this paper, an Earth-bounded coordinate system is used. The Earth-bounded coordinate systemi,j,k is centered in the muzzle, with the axes i, j, k pointing to fixed directions in space. Axes i is tangent to the Earth, j is orthogonal to i and runs against the gravity, and k is orthogonal to both i and j, setting up a right-handed trihedron. The model is illustrated in Figure 2.Figure 2
Mass point model with 3-DOF.With the aforementioned parameters, the equilibrium of forces in this case can be described with the formula(4)mdv→dt=g→+D→, where m is the total mass of the projectile.For setting up the system of equations, letx,y denote the projectile position and u,w the velocity, with u determining the horizontal and w the vertical projection of the velocity vector.Lett denote the time, 0≤t≤tf, with t=0 the initial time and t=tf the final time.The system of equations can be written as(5)dxdt=u,dydt=v,dudt=-Cxv2cos⁡φ,dwdt=-g-Cxv2sin⁡φ, where(6)v=u2+w2 is the radial velocity and φ is the angle between the thrust vector and the x-axis: particularly(7)φ=dydx.


## 2.1. Exterior Ballistics

The ballistic model is principally based on Newton’s law and the equations of motion are considered to be under the effect of air drag and the force of gravity only. Additionally, rockets have a thrust vector impelling the projectile for a few seconds (generally, combustion gases have a velocity range of 1800–4500 m/s [18]). Anyhow, rockets as well as mortars have ballistic trajectories and the object is to identify the threat on the basis of different flight characteristics.Letg→ denote a reference acceleration (acceleration of gravity at sea level on Earth), with(1)g=9.80665m/s2, taking effect on the point mass in vertical direction.The air dragD→ can have different values, depending on the design of the projectile, that is,(i)
muzzle velocityv0,(ii)
weight,(iii)
aerodynamics, and the properties of air, for example, (i)
density,(ii)
temperature,(iii)
wind,(iv)
speed of sound.Considering the general formula(2)D→=12·Cd·A·ρ·v·v→,Cx=Cd·CA·Bcontaining all parameters named above with (i)
A: cross-section area of the projectile,(ii)
ρ: air density,(iii)
v: velocity of the projectile,(iv)
Cd: air drag coefficient,(v)
CA: environmental properties,(vi)
B: ballistic coefficient, it is operative to find an appropriate approximation, so that the projectile can be specified. The parameters A, ρ, Cd, CA, and B are unknown, whereas v can be defined precisely from the measured radar data.The air drag coefficientCd for instance depends on the critical velocity ratio, pictured in Figure 1. Since the drag coefficient does not vary in a simple manner with Mach number, this makes the analytic solutions inaccurate and difficult to accomplish.Figure 1
Characteristics of the air drag coefficientCd.One can see from this figure that there is no simple analytic solution to this variation. With computer power nowadays, we usually solve or approximate the exact solutions numerically, doing the quadratures by breaking the area under the curve into quadrilaterals and summing the areas. In general, there are three forms of the drag coefficient:(1)
ConstantCd that is useful for the subsonic flight regime: Ma<1(2)
Cd inversely proportional to the Mach number that is characteristic of the high supersonic flight regime: in this case, Ma≫1(3)
Cd inversely proportional to the square root of the Mach number that is useful in the low-supersonic flight regime: Ma≥1 Carlucci and Jacobson [19] give a detailed description of the air drag coefficient.Another coefficient in common use in ballistics is the ballistic coefficientB, which is defined as(3)B=md2, where m and d are the mass and diameter of the projectile [19]. Section 3.2 deals with the problem of estimating the unknown parameters.


## 2.2. Equations of Motion

The aerodynamics and ballistics literature are quite diverse and terminology is far from consistent. This has particular significance in the coordinate systems used to define the equations of motion. Nevertheless, this field of research has a long history and a lot of approaches. More details are discussed in [20–24].In this paper, an Earth-bounded coordinate system is used. The Earth-bounded coordinate systemi,j,k is centered in the muzzle, with the axes i, j, k pointing to fixed directions in space. Axes i is tangent to the Earth, j is orthogonal to i and runs against the gravity, and k is orthogonal to both i and j, setting up a right-handed trihedron. The model is illustrated in Figure 2.Figure 2
Mass point model with 3-DOF.With the aforementioned parameters, the equilibrium of forces in this case can be described with the formula(4)mdv→dt=g→+D→, where m is the total mass of the projectile.For setting up the system of equations, letx,y denote the projectile position and u,w the velocity, with u determining the horizontal and w the vertical projection of the velocity vector.Lett denote the time, 0≤t≤tf, with t=0 the initial time and t=tf the final time.The system of equations can be written as(5)dxdt=u,dydt=v,dudt=-Cxv2cos⁡φ,dwdt=-g-Cxv2sin⁡φ, where(6)v=u2+w2 is the radial velocity and φ is the angle between the thrust vector and the x-axis: particularly(7)φ=dydx.


## 3. Concept

The purpose of the software is the calculation of trajectories. It receives the measured position of the projectile from the tracking radar and returns the predicted trajectory.A C code was written for the simulation and a GUI eases the handling of the results. Radar data can be read in and will be plotted for comparison.This chapter gives an overview of the methods used in this paper. An integration method for differential equations is introduced, which is used to solve the equations of motion in the previous section.


### 3.1. Integration Method

There are several integration methods implemented, all providing better results compared to the analytical methods used in [25].In this paper, the equations of motion are basically calculated with explicit, fixed step-size Runge-Kutta integration techniques. The advantage of this scheme over other schemes is that the approximating problems that result can be solved very efficiently and accurately. More details are discussed by Ramezani [26].Knowinghi=ti+1-ti the algorithm can be programmed on the analogy of [27](8)xi+1=xi+16ti+1-tik1+4·k3+k4 with(9)k1≔fxi,k2≔fxi+hi2k1,k3≔fxi+hi4k1+k2,k4≔fxi-hik2+2k3. With a global discretization error Oh4, the algorithm offers a tradeoff between high computing speed and best possible results [28].One may also select the Euler method in the program. Comparing to Runge-Kutta, the results are less precise due to a lower order of consistency. Anyhow, the Euler method achieves a significant improvement of computing time in the majority of cases with a global discretization errorO(h) [29].


### 3.2. Iterative Optimization

In the mathematical model which has been described in Section2.2, there are a number of parameters missing. The other variables are given and can be easily obtained through the measured trajectory elements. In order to determine the air drag D→ with the most accurate precision, the following algorithm was developed.The air drag is chosen in a way so that the exterior ballistic model complies to the measured trajectory of the projectile in the best possible way. This implies that the sum of the deviations between the calculated and the measured mortar positions should be minimal:(10)ε=min⁡fCx=min⁡∑i=1Nxim-xia2+yim-yia2+zim-zia2. The index m refers to the coordinates which are measured by the radar, while the index a belongs to the coordinates which are calculated by using numerical methods. The total amount of measurements is called N.Consequently, this is a nonlinear optimization. The objective function contains parameterCx. In order to find the optimum, one of the fastest methods of one-dimensional optimization, the so-called “Golden Section Search,” is applied. It only needs one value of the objective function for each step of the calculation. The second value is taken from the preceding iteration step. This method possesses a robust and linear convergence speed to find the minimum of a unimodal continuous function over an interval without using derivatives.The method chooses two pointsu1<u2 on the section a,b considering Golden Section:(11)u1=a+b-a·3-52,(12)u2=a+b-a·5-12. If the inequality f(u1)<f(u2) is complied, the minimum is in the interval a,u1. In any other case, it will be found on the stretch u1,b. When this procedure is repeated, the interval can be shortened again. In case of a new partition a,u2, there are new boundaries u1∗, u2∗ with u2∗=u1. Therefore, only two values of the goal functional are needed to be measured during the first step of the calculation [30].The goal is an optimal reduction factor for the search interval. Additionally, a minimal number of function calls are necessary [31].Golden Section Search enables an iterative adjustment of the trajectory in each step by using the calculated parameterCx for every previous iteration. Therefore, prediction gets more precise in the course of time.The programming flowchart is illustrated in Figure3.Figure 3
Programming flowchart.


## 3.1. Integration Method

There are several integration methods implemented, all providing better results compared to the analytical methods used in [25].In this paper, the equations of motion are basically calculated with explicit, fixed step-size Runge-Kutta integration techniques. The advantage of this scheme over other schemes is that the approximating problems that result can be solved very efficiently and accurately. More details are discussed by Ramezani [26].Knowinghi=ti+1-ti the algorithm can be programmed on the analogy of [27](8)xi+1=xi+16ti+1-tik1+4·k3+k4 with(9)k1≔fxi,k2≔fxi+hi2k1,k3≔fxi+hi4k1+k2,k4≔fxi-hik2+2k3. With a global discretization error Oh4, the algorithm offers a tradeoff between high computing speed and best possible results [28].One may also select the Euler method in the program. Comparing to Runge-Kutta, the results are less precise due to a lower order of consistency. Anyhow, the Euler method achieves a significant improvement of computing time in the majority of cases with a global discretization errorO(h) [29].


## 3.2. Iterative Optimization

In the mathematical model which has been described in Section2.2, there are a number of parameters missing. The other variables are given and can be easily obtained through the measured trajectory elements. In order to determine the air drag D→ with the most accurate precision, the following algorithm was developed.The air drag is chosen in a way so that the exterior ballistic model complies to the measured trajectory of the projectile in the best possible way. This implies that the sum of the deviations between the calculated and the measured mortar positions should be minimal:(10)ε=min⁡fCx=min⁡∑i=1Nxim-xia2+yim-yia2+zim-zia2. The index m refers to the coordinates which are measured by the radar, while the index a belongs to the coordinates which are calculated by using numerical methods. The total amount of measurements is called N.Consequently, this is a nonlinear optimization. The objective function contains parameterCx. In order to find the optimum, one of the fastest methods of one-dimensional optimization, the so-called “Golden Section Search,” is applied. It only needs one value of the objective function for each step of the calculation. The second value is taken from the preceding iteration step. This method possesses a robust and linear convergence speed to find the minimum of a unimodal continuous function over an interval without using derivatives.The method chooses two pointsu1<u2 on the section a,b considering Golden Section:(11)u1=a+b-a·3-52,(12)u2=a+b-a·5-12. If the inequality f(u1)<f(u2) is complied, the minimum is in the interval a,u1. In any other case, it will be found on the stretch u1,b. When this procedure is repeated, the interval can be shortened again. In case of a new partition a,u2, there are new boundaries u1∗, u2∗ with u2∗=u1. Therefore, only two values of the goal functional are needed to be measured during the first step of the calculation [30].The goal is an optimal reduction factor for the search interval. Additionally, a minimal number of function calls are necessary [31].Golden Section Search enables an iterative adjustment of the trajectory in each step by using the calculated parameterCx for every previous iteration. Therefore, prediction gets more precise in the course of time.The programming flowchart is illustrated in Figure3.Figure 3
Programming flowchart.


## 4. Simulation Results

An example is calculated for a rocket Type 63 HE on a common Intel x86. The data specification of the rocket is listed in Table1.Table 1
Rocket Type 63 HE specification [17].
Maximum range
8.5 km
Overall length
839.0 mm
Caliber
107 mm
Cross-sectional area
7.72 cm2
Weight
18.84 kg
Lateral moment of inertia
0.98122 kgm2
Longitudinal moment of inertia
0.03135 kgm2
Position of center of gravity
395.8 mm
Standard empty weight
8.496 kg
Combustion duration
0.6 sec
Impulse
6.7 kNs
Flight time
21.5 secThe trajectory was recorded by a Weibel radar of the type MFDR-2100/35. It detects the RAM target at high accuracy. It is designated to the RAM target with information received from a search radar. The accuracy is listed in Table2. The Kalman filter illustrates the track error over time.Table 2
Weibel radar accuracy [17].
Accuracy 2.0 km
Accuracy 4.0 km
Time (ms)
Max. range (km)
Rg (m)
Az (mills)
El (mills)
Rg (m)
Az (mills)
El (mills)
10
8.9
0.20
0.491
0.491
0.81
1.965
1.965
20
13.2
0.09
0.220
0.220
0.36
0.879
0.879
Kalman filter
500
23.6
0.03
0.069
0.069
0.12
0.278
0.278Lett0=3 s be the time at which the radar detects the bullet. The final flight time is reached after tf=21.5 s.The duration of calculation is adjustable at will. More tracking points will certainly help to get better results, but sometimes a fast interception of the RAM threat is indispensable. Starting the forecast with a 3-second period of calculation, there will be a mean square deviation of 566.9 m between the calculated and the real trajectory. By now it is possible to identify different RAM targets by regarding the predicted trajectory characteristics. This example is illustrated in Figure4. After 6 seconds of calculation, the mean square deviation is reduced to 181.8 m. There are inaccuracies in all axes of coordinates. The estimation of the calculated parameter Cx needs more iterative calculation steps at this point.Figure 4
Forecast with the calculation period 3–6 s.Finally, after 12 seconds of calculation, the mean square deviation is reduced to 32.2 m and there is still enough time for the command and control system to initiate all necessary steps, for example, warning and defending. The results are shown in Figure5. It is quite obvious that the simulated altitude is overestimated. Mortar grenades have a strong change in altitude on their trajectory, a challenge for simulation-based early prediction systems.Figure 5
Forecast with the calculation period 3–15 s.The prediction of the trajectory allows the calculation of the point of impact. The future position of the projectile is calculated through extrapolation of the measured values.It is clear that the prediction gets significantly better with each iteration. Thus, certain areas in the field camp can be warned partially and a counterattack can be initiated. The more the radar data available for the analysis, the closer the prediction to the measured trajectory. More tracking points will certainly help to get better results, but sometimes a fast interception of the RAM threat is indispensable.With a prediction of 3 seconds into the future, for example, which corresponds to an intercept range of almost 3000 m, the computational error at the point of impact(13)εp=xpm-xpa2+ypm-ypa2+zpm-zpa2 is smaller than 3 m. Here, the index p refers to the coordinates that are measured and calculated at the point of impact. Details and more examples are discussed by Ramezani et al. [10].


## 5. Summary and Outlook

This paper introduces an algorithm for early warning systems used for command and control applications in out-of-area missions and is based on the MONARC (modular naval artillery concept) project. The basic methods have been tested successfully and they are used in fire guidance solutions for German frigates of type 124 and 125.The most important aspect is that one can distinguish between different projectiles in order to predict the trajectories and hit points more accurately. To calculate their trajectories, different flight phases are analyzed and the designs of the projectiles are estimated by the use of iterative optimization methods for approximating environmental and ballistic properties.Future work concentrates on giving the user specific information of the projectile data. Further work has also to be done on a 3-dimensional simulation.At the end, sophisticated simulation software will be established through which it will be possible to show and evaluate a real-time battlefield scenario.



---

*Source: 8157319-2017-08-07.xml*