Variational Calculus Approach to Optimal Interception Task of a Ballistic Missile in 1D and 2D Cases

: The paper presents the application of variational calculus to achieve the optimal design of the open-loop control law in the process of anti-ballistic missile interception task. It presents the analytical results in the form of appropriate Euler–Lagrange equations for three different performance indices, with a simple model of the rocket and the missile, based on the conservation of momentum principle. It also presents the software program enabling rapid simulation of the interception process with selected parameters, parametric analysis, as well as easy potential modiﬁcation by other researchers, as it is written in open code as m-function of Matlab.


Introduction
Interception of ballistic missiles has been a research topic for years.Calculation of interception trajectories is usually performed off-line by using two-point boundary-value problem formulation.The paper [1] employs the proposed guidance law in real-world conditions, where interception trajectories are sought in R 6 , comprising of: three NED (north-east-down) coordinates, speed, path angle and heading of the interceptor.This usually calls for the use of optimization techniques to obtain the reference trajectories to follow.In [2], there are four different guidance laws considered, showing that the interception of the missile is possible if one has perfect measurements and accurate estimation of the target maneuvers.To obtain the optimal interception point, one can use multiple approaches, ranging from classical variational calculus as in this paper, or others, such as metaheuristics, e.g., a genetic algorithm, see [3], which is used there to establish the optimal interception point, minimizing certain functional.The functional may refer to flight time of the intercepting missile, or distance between the interceptor and the ballistic missile.Guidance of an intercepting rocket is a problem which has multiple solutions.One can distinguish passive and active guidance approaches [4].Several types of active guidance methods can be distinguished: • line of sight guidance (also called a two-point guidance), where the intercepting rocket always points to the target, • deviated pursuit, where the intercepting missile points ahead of the target, • proportional navigation, where the line of sight angle of the target is estimated, • optimal guidance by minimizing a certain functional.
Among the last, one can consider many approaches, such as LQR [5] embedded into variational calculus framework, LQI (LQ integrated) control law [6], applied to linearized systems, predictive control and nonlinear predictive control, which fall into the optimization strand [7], or performance-based approaches to PID tuning, [8].This paper adopts the optimization technique to generate the optimal open-loop control signal to reach the desired point minimizing a certain performance index.A more comprehensive view on a general framework of a ballistic missile defense systems, comprising of a launcher and a radar, please refer to [9].Another study of tracking a ballistic missile can also be found in [10].
Optimization techniques have been developed for ages-it was in 17th century when such mathematicians as Leibnitz, Bernoulli, Newton of d'Hospital pushed their horizons further.They are required in production processes, business/management, or in military applications, since everybody is interested in gathering as much as possible, at the least expense, in the shortest time, etc.The aim of this paper is to use this approach in the optimal interception problem of a ballistic rocket with a known trajectory, using the intercepting missile.
As has been already remarked, an optimization approach is attractive in many areas, and can include, e.g., quadratic cost function in an unbounded problem, leading to an easy solution.In real world applications, however, constraints are ubiquitous, and need to be taken into consideration during the optimization process.In this work, optimization is based on variational calculus that could be, e.g., presented as the problem of finding, among all the continuous functions x(t), the one that satisfies terminal conditions [11]: achieving local extremum of the functional The latter task has fixed terminal conditions, and the necessary condition of existence of such an extremum is satisfying the Euler-Lagrange (E-L) equation From now on, whenever applicable without misunderstandings, the time indexes will be omitted.In the case of a variable initial point lying on the curve φ(x), the function x(t) is sought that crosses (t 0 , x 0 ) and intersects φ(x) at time t f .This function satisfies, in addition to E-L equation, the following transversality conditions: If the functional depends on x, then the necessary conditions are formed as a set of E-L equations.In the considered interception task, φ might depict the route traced by the nuclear train from what the missile is launched at t 0 , or the path of the rocket to be hit at t f .
Finally, when constraints are present and can be put in the form g(t, x, u) = b, where u denotes the control signal, and g is some function satisfying continuity conditions, the functional can be put in the form The variational approach presented above will be adopted for finding optimal control signals of a missile to intercept a ballistic rocket, for which model and parameters will be presented in the next Section.The anti-ballistic missile is abbreviated in the paper as ABM, and is usually a surface-to-air missile designed to counter ballistic missiles [12,13].
As optimal interception tasks have already been solved and considered in the literature, presenting various scenarios and solutions, though as a set of pure results only, and form a continuously-developed research field, this paper presents the software in open text which can be used not only to simulate interception task, taking multiple system parameters into account during analysis, but also to perform parametric analysis on selected performance indices.All the issues will be addressed in the paper, which is the main motivation and contribution of this work.

Simplified Intercepting Missile Model (1D)
Prior to giving a model of a rocket the question of model accuracy should be addressed.It would seem advantageous to find the most accurate model, but in such a case, computer calculations would take too much time, and mathematical description would disable any actions.The model should be simplified enough to mirror the dynamics of the missile, but not be oversimplified.
The missile is assumed to move on the basis of propulsion force, where gases leaving the rocket move it in the opposite direction [14][15][16].In order to enable mathematical description, the momentum conservation rule has been used, which states that the total momentum of the rocket (p rocket (t) + p gas (t)) always adds to zero, where with m being the mass of the rocket, and v(t) being its velocity.The effect of the drag force is harder to quantify, nevertheless, since in many important applications drag effects are very small, the drag force is omitted here.Potential drag effects, where the density of air changes with altitude, usually result in equations that are impossible to be integrated analytically (numerical integrations are needed), and leads to, in typical flight conditions of a 12,000 kg rocket flying at 700 m/s, a drag force at the level of 2% of gravity force.The considered control system should work in the following way: the control signal w(t) is the momentum of the gas, which when divided by m gives the velocity of the rocket.By integrating this velocity one can obtain information about the position of the rocket that should be compared with the information about the trajectory of the rocket that must be intercepted.The location is measured by a sensor, represented by pure feedback gain, and this signal should be subtracted from the reference signal w to obtain tracking error.The controller itself is not needed, since the open-loop control signal will be obtained as the result of optimization.
Assuming that m r is the mass of he missile, k 2 is the gain in the inner feedback loop of the rocket, w(t) is the reference signal (optimized input signal to the control system), and y(t) is its output, one can propose the following model with k 1 = 1 m r (see Figure 1): The time constant of the rocket is assumed to be approximately several minutes, and its mass approximately 1000 kg [13,17,18].

Simplified Rocket Model (1D)
Typical ballistic rocket is flying approximately horizontally the middle phase of its flight, plotting parabolic ∩-shaped function.The parabola should have a flat graph with some maximum attitude assumed, as in this paper, i.e., 100 km.Typical parameters of the intercepting missiles are given in the Tables 1 and 2 [12].The trajectory of the rocket might be defined by the following canonical quadratic function where p [h] is the time in which the missile reaches its maximum altitude, a is the slope, and q [km] is the maximum altitude [14,16].
The missile used to intercept the ballistic rocket does not have any information about the time when the latter has been launched, thus t = 0 denotes the launch time of the intercepting missile.Assuming, as an example, that p = 0.25 h with arbitrary chosen a = −10, the default canonical ballistic rocket equation in 1D becomes Such a 1D case can be interpreted as a situation when the interception takes place directly above the place of installation of intercepting rockets, below expected ballistic missile trajectory.

Intercepting Missile Model (2D)
The 2D task is the actual point of interest, where translation of the missile takes place in a two-dimensional space (distance x vs. altitude y).It is assumed that the observer is in parallel to the hyperplane on which the translation of the missile takes place, as well as the translation of the rocket, and their observation angle is the right angle, as depicted in Figure 2. By assuming a specified localization of the observer, there is no need to present the translation in 3D space, thus z axis information can be omitted.
Generalizing the previous 1D momentum equation to the 2D case [15], where velocity vector v is decomposed into xand y-components.
With the base vectors e 1 and e 2 , the momentum p = p x e 1 + p y e 2 .

Rocket Model (2D)
The ballistic rocket that should be intercepted is assumed to move as in the 1D case on a paraboloid, and its range is assumed to be 2000 km by default and maximum height in y coordinate-100 km, i.e., it is a medium-range missile according to the American division or Russian classification [12], see Tables 3 and 4.Among American ABMs currently used one can find: MIM-104 Patriot (PAC-1 and PAC-2), U.S. Navy Aegis combat system using RIM-161 Standard Missile 3, or Russian Federation's missiles, as ABM-1 Galosh, ABM-3 Gazelle, ABM-4 Gorgon, and rockets of SAM system, like S-300P (SA-10) to S-400 (SA-21).
A rocket's trajectory along the y axis is described as in the 1D case, i.e., by the canonical form of a quadratic equation with the same description, whereas along the x axis as where g is the velocity of the missile along the x axis.As in a 2D case, the to-be-presented plots depict the interception phase taking place in the middle phase of flight, where missile trajectory is relatively flat.Just as in the case of a 1D task, in order to generate control action, no information concerning its launch time is available, thus t = 0 denotes launch of our missile time again, all parameters of flight along the y axis remain the same as in the 1D case.The parameter g is set as default to 400 kmph.
The interception process should take place when the ballistic rocket is in its descending phase of flight, whereas the missile continuously increases its altitude, as presented in Figure 3.

Performance Indices
In order to determine transversality conditions, and to mimic different behaviour of the system, three performance indices have been considered.
At all times, it is assumed that the optimal reference (input) w to the system is generated, to exert appropriate changes in x, y, to intercept a rocket.The performance indices presented below, refer to energy consumption or to the time-of-interception, which both have well-defined meanings.

Index J 1
The first performance index refers to minimization of energy of the input signals to the guidance system, since our control is open-loop based, and the problem is of generation of optimal input.In such a case, where s is the weight between the components.The input signal should satisfy dynamic equations of the rocket, and in both dimensions it should be the same, thus: where w x denotes the position along the x axis, and w y along the y axis, T is a time constant, and k is the gain.The function under the integral takes the form: Euler equations derivatives: giving the E-L equation which solution becomes . By performing similar calculations in the y axis, one obtains y(t) = C 3 exp − t T +C 4 exp t T .Since we have five unknowns (C 1 , C 2 , C 3 , C 4 , t f ) the following boundary conditions must be used: Finally, a single equation is obtained 3.2.Index J 2 In the second performance index, only energy along the y axis is considered, since along the x axis the trajectory moves according to the linear equation of motion x = gt, thus where By performing similar calculations as in the case of J 1 one obtains the E-L equation and thus y(t) = C 1 exp − t T +C 2 exp t T .The border conditions become: Finally, a single equation is obtained This index is connected to minimum-time control that requires the velocity constraint ẏ = v max to be introduced, otherwise the input command becomes infinite.The performance index is of the form There are no constraints imposed on the translation in x axis, since the movement in x axis is proportional to time (ϕ x (t) = gt), thus minimum time implies minimal value of ϕ x .Taking this constraint into consideration, the performance index takes the form Euler equations derivatives for y: and for λ: yield in the following set of E-L equations: with the solutions: Border conditions become: and, as a result, From the latter equation, one can obtain the interception time.In the x axis, it is assumed that the equation of motion becomes with border conditions: x(0) = 0 , (51) where t f is already calculated for y.

Simulation Parameters
In all the considered interception problems, a single equation binding interception time has been obtained that needed to be solved, using numerical methods with sampling time T S = 0.3 s, resulting from observation of time constants of the obtained plots and on the basis of multiple simulations, and observations of performance indices behaviour.
The following parameters have been set as default values: It is to be stressed that parameter s has been chosen in such a way so as to obtain interception time in minutes, to avoid false results.

Simulation Results
Below, the optimal curves for performance indices considered are included: In Figures 4-6, the visualisation of the obtained results is presented in y axis, to show the interception takes place at the calculated time instants, whereas in Figures 7-9, the reference inputs are shown, respectfully.

Software Program to Calculate Optimal Interception Time
The software program with the GUI in Matlab has been created that solves using Symbolic Toolbox the equations defining interception time instants t f , calculates integration constants, and visualises the results.
Its main window is presented in Figure 10, where it is possible to define the task giving p, q, T, V max , choosing the appropriate performance index to minimize (J 1 , J 2 , J 3 ) and dimensionality of the problem (1D, 2D).After pressing Simulation button, the calculations are executed, with the final plot drawn at the end, see Figure 11.Optionally, it is also possible to perform parametric analysis.where in different windows one can choose simulations for three different values of parameters from the list: p, q, T, V x , for which appropriate plots of x, y, w x , w y are automatically drawn, with performance indices and interception time instants printed.
The software program will be made available in the open code, as m-files, in the MMD Database (see [19]), to enable further development of the program by interested researchers, and on the web page of Applied Control Techniques research group (http://act.put.poznan.pl).It can also be found at the link in the Supplementary section.
In order to make use of the derived equations and optimal trajectories, one can perform a sample parametric analysis.For example, Figure 12 presents results from 2D case with T x = 1 min, T y = 2 min, and J 2 index.In all the cases, a = −10, k = 0.25, g = 400 kmph, s = 0.45.

Figure 1 .
Figure 1.Model of the rocket control system.

Figure 10 .
Figure 10.Main menu of the software program.

Figure 11 .
Figure 11.Visualisation of the interception process for default parameters.

Figure 12 .
Figure 12.Interception time t f as a function of p and q for minimized J 2 index (2D task).

Table 1 .
Standard parameters of interception missiles.

Table 2 .
Selected parameters of sample interception missiles.