Abstract

We consider the integrated problem of allocation and control of surface-to-air-missiles for interception of ballistic targets. Previous work shows that using multiple missile and utilizing collaborative estimation and control laws for target interception can significantly decrease the mean miss distance. However, most of these methods are highly sensitive to initial launch conditions, such as the initial pitch and heading angles. In this work we develop a methodology for optimizing selection of multiple missiles to launch among a collection of missiles with prespecified launch coordinates, along with their launch conditions. For the interception we use 3-DoF models for missiles and the ballistic target. The trajectory of the missiles is controlled using three-dimensional extensions of existing algorithms for planar collaborative control and estimation laws. Because the dynamics of the missiles and nature of the allocation problem is highly nonlinear and involves both discrete and continuous variables, the optimization problem is cast as a mixed integer nonlinear programming problem (MINP). The main contribution of this work is the development of a novel probabilistic search algorithm for efficiently solving the missile allocation problem. We verify the algorithm by performing extensive Monte-Carlo simulations on different interception scenarios and show that the developed approach yields significantly less average miss distance and more efficient use of resources compared to alternative methods.

1. Introduction

Missile systems form an integral part of a nation’s defensive capability. In particular, high precision control of surface-to-air-missiles (SAMs) is a key technology for intercepting maneuvering targets, such as hostile aircraft and ballistic missiles [1]. Successful mid-air interception of such targets might prevent significant losses and enable preparing for counter-attacks. Hence, there has been a significant amount of previous work that focuses on development of guidance, control, and estimation laws for interception of maneuvering targets by a single missile in the past years [2].

Intercepting ballistic targets is generally more challenging compared to interception of static targets or targets executing simple maneuvers (such as an aircraft in level flight), due to lack of precise dynamic models and high reentry speed of the target. Hence, even the smallest perturbations in the guidance errors might translate into large miss distances for the interception of ballistic targets. In order to reduce the miss distance, there has been a growing interest in using multiple collaborating missiles for the interception of a single target. There are two main advantages for using multiple missiles; () probability to hit increases due to increased number of missiles targeted at the hostile; () missiles can fuse their measurements to improve the estimation of the target’s state, which results in increased guidance performance and hence reduced miss distance.

That being said, the initial conditions of the launched missiles still have a significant impact on their overall performance. Ballistic target’s total speed is usually greater than the interceptors speed; hence the launching coordinates of the SAMs directly affect the miss distance. In addition, launching pitch and heading angles of the SAMs are also critical for the guidance performance. Even though errors in these angles are compensated by the missile autopilots in later phases of the flight, due to high kinetic energy of the ballistic target, having 1-2 degrees off in the initial pitch and heading angles of the SAMs might lead to significant deviations from the desired performance, even in the presence of restoring actions of the autopilot. Hence, the selection of the launching coordinates of the SAMs along with their initial launching angles should be adjusted carefully to gain optimal collaborative missile interception performance. The main contribution of this paper is the development of an optimization algorithm that drives the initial conditions of collaborative guidance and estimation laws. In particular, the missile launching parameters and the number of launched missiles are selected using a probabilistic search algorithm, which attempts to optimize an objective function that favors minimum miss distance and maximum efficiency in use of resources.

1.1. Previous Work on Collaborative Estimation and Control for Missiles

There has been a significant amount of previous work on control and estimation of multiple missiles. Chen and Speyery [3] formulated the multiple missile coordination problem as a Linear Exponential Gaussian Differential game and applied their algorithm to interception of ballistic missiles in terminal and boost phase. Wang and Fu [4] formulated the multiple missile interception problem as a multiplayer pursuit and evasion game and studied the interception of a ballistic target in three dimensions. Jeon et al. [5] developed a cooperative proportional navigation (PN) law that enables multiple missiles to close simultaneously on a stationary target. Similarly, Daughtery and Qu [6] also developed an algorithm for multiple missiles attacking simultaneously a target, and they also showed that their algorithm is robust to communication losses between the missiles.

Shaferman and Oshman [7] developed an extended Kalman filter (EKF) algorithm that fuses information gathered from multiple interceptors. They were able to show that using cooperative estimation algorithms yield improved guidance and control performance. Shaferman and Shima [8] combined adaptive control laws and multiple model filtering algorithms for collaborative interception. Liu et al. [9] also considered using multiple EKFs to improve guidance performance; in particular they discovered that the estimation performance improves as the relative line of sight (LOS) between two intercepting missiles gets larger and designed a control law that enforces separation between two missiles. Recently, Shaferman and Shima [10] also developed guidance laws for enforcing relative intercept angles. Another recent development was provided by Wang et al. [11], where the authors used a probabilistic framework to maximize hit-to-kill probability for two-missile cooperative interception.

These previous works showed that fusing the estimation process between multiple missiles almost certainly leads to improved control performance and hence reduced miss distance. Most of these works only study the scenario where missiles are about to close on the target; the launch conditions are largely ignored. However, as explained in the beginning of this section, estimation and control performance might vary significantly under different launch conditions, such as different initial heading and elevation angles of the missiles. Moreover, launch conditions could also require adjustments based on relative heading and velocity of the ballistic target. In particular, the sensitivity of miss distance based on different conditions would be higher for targets with high speeds, such as ballistic missiles.

1.2. Previous Work on Missile Allocation

The existing algorithms for multiple missile interception assume that a fixed number of missiles have been assigned for the target interception. In real scenarios, usually a larger set of missiles are available on the ground, and a certain subset of these missiles should be allocated and launched for interception. Deciding on which missiles to allocate depends on several factors, such as the altitude and the velocity of the target. Since increasing the number of allocated missiles improves the estimation performance and hence the probability of successful interception, it might be desirable to launch many missiles as possible. However, launching more missiles than necessary would result in inefficient use of resources; hence a trade-off exists between the kill probability and the number of launched missiles.

To the best of our knowledge, there are no previous works that directly address the allocation problem for multiple missiles that utilize collaborative estimation and control algorithms. However, previous researches on weapon and task allocation methods are certainly relevant to this problem, and they are reviewed in this subsection.

Weapon and task allocation algorithms generally address the problem of assigning a finite number of weapons to multiple targets. The search space for such problems are usually huge and the objective functions lack useful properties such as linearity or convexity; hence evolutionary algorithms and probabilistic search approaches are very popular for these problems. Lee et al. [12] used a genetic algorithm based on greedy eugenics to solve the assignment problem, whereas Teng et al. [13] utilized particle swarm optimization approach for the coordinated air combat missile-target assignment. Ant colony optimization [14] and simulated annealing [15] have also been applied to the target assignment problem. Besides genetic algorithms, game theoretic [16] and rule-based [17] approaches have also been applied to missile allocation problems.

These existing works either completely ignore the vehicle dynamics or use simplified missile and target models for solving the allocation problem. However, utilization of collaborative estimation and control algorithms is integral to the interception of ballistic target problem studied in this paper; hence the dynamics of estimation and control algorithms cannot be ignored by the allocation algorithm. Existing works usually remedy this problem by assigning a probability of success to each possible assignment; however overestimating or underestimating these probabilities might lead to serious performance problems.

1.3. Overview and Contribution of This Work

The main idea behind the developed algorithm is to cast the missile allocation problem as a mixed integer nonlinear program (MINP), where the discrete decision variables correspond to selection of which missiles to launch, and the continuous decision variables correspond to launch conditions such as elevation and heading angles for each missile. In order to solve this optimization problem, a novel probabilistic search method is developed, which extends the well-known continuous optimization algorithm Covariance Matrix Adaptation (CMA) [18] to solve problems that contain both continuous and discrete variables. The optimization problem also embeds the full dynamics of the missiles, where we have extended an existing collaborative EKF algorithm [7, 10] to work in three-dimensional settings.

With respect to the state of the art, our algorithm offers the following contributions:(i)The previous work on missile control assumes fixed initial launch conditions for collaborative targets. Our algorithm automates the process of selecting missile heading and elevation angles, which impacts the estimation and control performance in later stages and hence reduces the average miss distance.(ii)Besides optimizing the launch conditions, the developed algorithm also handles the selection of the subset of missiles for interception among the available larger set of missiles by examining the trade-off between the missile launch cost and achievable miss distance. This feature is captured in MINP formulation, and the results show that the developed algorithm yields lesser miss distances using lesser number of missiles, compared to alternative approaches.(iii)The developed optimization algorithm works by propagating missile and target dynamics equipped with estimation and control laws. Hence unlike existing missile allocation algorithms, there are no simplifications or approximations in dynamical models, which translates into better prediction of overall system performance.(iv)As a minor contribution, we would also like to note that most existing collaborative missile control papers consider two-dimensional (planar) missile and target dynamics. In this work we have extended the collaborative EKF to work in three dimensions. Although this is a straightforward extension, having the results in three dimensions leads to more realistic miss distances and hence enables a better assessment of the fitness of the developed algorithms on real battlefield operations.

The paper is structured as follows. In Section 2, the ballistic target interception problem is formulated and the missile modes for both the interceptors and the ballistic target are provided. Section 3 explains the collaborative estimation and control algorithms used in this work, and Section 4 gives the details of the novel optimization algorithm used for solving the allocation problem. Finally, Section 5 studies the performance of the algorithm by examining the results of Monte-Carlo simulations that involve various different target initial conditions.

2. Problem Formulation and Missile Models

In this section we provide the complete formulation of the optimization problem for the multiple missile ballistic target interception. Since the problem formulation involves the dynamics of the missiles, the model of the interceptor missiles and the ballistic target are also given in this section.

2.1. Multiple Missile Ballistic Target Interception Problem

Let denote the state vector of the ballistic target and let , , denote the state vectors of the SAMs, where denotes the position, denotes the velocity, and denotes the number of missiles. It is assumed that SAMs initial positions are fixed at the launch site and the ballistic target is detected by a radar in its terminal phase. We assume that closed loop guidance and control algorithms are implemented on SAMs (see Section 3) and the decision maker only needs to decide on the initial launching conditions.

Let and denote the initial launching pitch and heading angles for each SAM. Let be the terminal time for the scenario. Define the miss distance as In other words, for a given set of launch conditions, miss distance is computed by propagating the dynamics of the SAMs and ballistic target over the horizon and finding the minimum distance achieved between the ballistic target and all launched SAMs.

The objective of the multiple missile ballistic target interception problem is to minimize the miss distance by optimizing over the launch conditions : where and are the dynamical models of the ballistic target and interceptor missiles (see Section 2.2). For the sake of simplicity, we assume that all SAMs share the same model.

Although the optimization problem above reflects the main objective of the mission, it assumes that all of the SAMs will be used for interception. However, in many scenarios only a strict subset of the available SAMs might be sufficient for successful interception. Hence we should modify the objective function to reflect the cost of launching a missile in order to yield solutions that minimize the miss distance while avoiding inefficient use of resources.

Let denote the binary decision vector, where means that SAM is launched for intercepting the ballistic target and means that SAM is not launched and will stay in its initial position throughout the mission. The modified optimization problem is formulated as where is a user defined parameter that reflects the cost for launching a missile. Larger values of favor launching lesser number of missiles. Since the optimization problem in (3) contains both continuous variables and discrete variables and since the dynamics and are nonlinear, it can be classified as a nonlinear mixed integer programming (MINP) problem. MINP are known to be very challenging to solve, and the major contribution of the paper is the development of a probabilistic search method (explained in Section 4) for solving the optimization problem in (3).

Note that we are only interested in deciding the allocation of the missiles and their launch conditions, and no further trajectory optimization is needed, since it is assumed that SAMs have closed loop estimation and guidance laws embedded in their dynamics. These estimation and control laws are further detailed in Section 3. However, before examining the control laws, we describe the open loop dynamics of the SAMs and ballistic target in the next subsection.

2.2. Missile Models

2.2.1. Interceptor SAM Model

The SAM interceptor used in this study is a 3-degree-of-freedom (DoF) point mass model of a missile, which is controlled via aerodynamic forces generated by the fins. Three-dimensional equations of motion of the interceptor can be stated as shown in [2]where is the downrange displacement of the SAM, is the cross-range displacement of the SAM, is the altitude of the SAM, is the gravity acceleration, is the flight path angle, is SAM velocity, and are the vertical and horizontal load factors.

2.2.2. Ballistic Target Model

In the 3D interception scenario, SCUD-B short range ballistic missile (SRBM) is used as the incoming ballistic threat. Simulation starts from a given apogee altitude, velocity, reentry angle, and heading. Aerodynamic data of the SCUD-B is obtained by using Missile DATCOM software [19].

The drag force and the gravity acceleration components that act on the ballistic target during reentry phase are shown in Figure 1. Drag force acts in the opposite direction of the velocity vector. Hence, if the effect of the drag force is greater than that of gravity, the ballistic target decelerates during its flight. The perpendicular component of the deceleration vector to the line of sight is shown as a target maneuver to the pursuing interceptor. Thus, the interceptor guidance law must take this deceleration maneuver of the ballistic threat into account [20].

Figure 1 

Forces acting on the ballistic target model.

The approximate mathematical model of the target in the reentry phase is given in [20]where is the drag force, , , and are the velocity components, is the target missile mass in the reentry phase, is the gravity acceleration, is the reentry angle, and is the heading angle.

3. Control and Estimation Algorithms

3.1. Radar Measurement for Initialization

For initialization of EoMs of the missiles and the estimation process, a ground based radar is used. The radar is fixed in a prespecified position. Relative positions of the interceptors and the target according to the radar in the 3D space are shown in Figure 2.

Figure 2 

Radar geometry.

We assume that the radar estimates the initial state vector without any error and send the initialization data to the interceptors without any delay

The geometric relation between the interceptor and the target is given in (7). In this work, it is assumed that and are estimated perfectly by the radar

3.2. Relative Kinematics

3D engagement geometry between the target and interceptors is shown in Figure 3. It is assumed that each interceptor missile can measure its own inertial state vector as shown in (8) and measurements are exact

Figure 3 

3D engagement geometry.

The state vector of the missile according to the target at time is given in (9) where . In the information sharing mode, each missile could transmit and receive the estimated target data without any loss and delay

Relative kinematics between the interceptors and the ballistic missile in , , and planes are defined by usingwhere

The relative kinematics equations can be described in the discrete-time by usingwhere is state vector of the missile at time and is obtained by integrating the relative kinematics EoMs in (10) [7].

3.3. Measurement Model

Each interceptor is equipped with the infrared seeker that measures the line of sight (LOS) angle . This measurement has a zero-mean Gaussian noise with standard deviation . LOS angle measurements are performed by each interceptor missile separately; therefore , . The missile LOS measurement is where The interceptor missiles can operate in two modes. The first mode is information nonsharing mode. In this mode, each interceptor measures only its own LOS angle and uses it in the estimation process of the relative states. The measurement vector in the information nonsharing mode is given as shown in (13).

The second mode is information sharing mode in which each interceptor not only measures its own LOS angle but also calculates the LOS angle of the other interceptors as shown in (15) [7]. Also, to improve the estimation quality, each interceptor calculates the range-to-go distance by using position and LOS angle data as shown in (16):

Here, is the missile which performs measurement and is the missile observed by the missile .

The LOS angle of the missile measured by the missile is calculated by using

In addition, the range-to-go of the missile measured by the missile is obtained by using

In the information sharing mode, the measurement model is given bywhere

3.4. Extended Kalman Filter

Because of the nonlinear relative kinematics between the interceptors and the ballistic target, an extended Kalman filter (EKF) is used for estimating the unmeasured data and to filter the noisy LOS angle measurements [7]. The prediction error covariance matrix is given in where is the transition matrix associated with the relative kinematics, is sampling time, and is the identity matrix with appropriate dimensions. is the covariance matrix of the equivalent discrete process noise and it is calculated as shown in is the Jacobian matrix associated with the nonlinear relative kinematicsThe state estimation is performed by usingwhere is the predicted state vector and is the Kalman gain as shown inHere, is the measurement Jacobian matrix and is the measurement noise covariance matrix. The covariance matrix is updated as shown in

3.5. 3D True Proportional Navigation Algorithm

The proportional algorithm is one of the most common and effective guidance techniques because of its simple structure and implementation. The true proportional navigation (TPN) system generates the acceleration command perpendicular to the LOS. As shown in (29), the acceleration command is a function of closing velocity and LOS rate where is the acceleration command perpendicular to the LOS, is closing velocity, and is navigation ratio which is generally between 3 and 5 [20]. In this 3D interception study, TPN algorithm is applied for , , and planes separately [21]. Geometry of relative kinematics for each different plane is displayed in Figures 4, 5, and 6.

Figure 4 

Relative kinematics on plane.

Figure 5 

Relative kinematics on plane.

Figure 6 

Relative kinematics on plane.

Acceleration commands in , , and are obtained as shown in The acceleration components of the interceptor in the -, -, and -axis () can be obtained from (30) by using the trigonometric relations: Before applying the control commands to the interceptor, vertical and horizontal components and should be calculated. Here, is in the pitch plane and perpendicular to the velocity vector of the interceptor and is perpendicular to both velocity vector and vertical acceleration vector. For TPN, these acceleration components are calculated using

4. Optimization Algorithm

4.1. CMA Algorithm

Consider the general form unconstrained optimization problem: It is well known that when possess a certain structure (such as being continuous, linear, or convex), there are variety of local search algorithms that can be applied to solve this optimization problem efficiently. However, when does not possess these desirable properties, local search methods either fail to find an answer or get stuck in local minima. Global search methods [22] remedy this problem by generalizing the search over the entire state space. Although global methods can also exploit the structure of , many global methods treat as a black box function, and hence the solution is found entirely by examining the input-output pairs .

Covariance Matrix Adaptation (CMA) [18] is a popular global search method that usually ranks among the best solvers in global search benchmarks [23]. The basic idea behind CMA is to place a multivariate normal distribution over the search space and sample candidate solutions from this distribution. The mean vector and covariance matrix of the distribution are incrementally updated at each step based on the values of the sampled solutions. The objective is to eventually steer the mean vector to the optimal solution and shrink the covariance matrix to identity matrix; hence in the limit the distribution will yield the optimal solution when it is sampled.

For completeness, we provide the pseudocode for the CMA algorithm in Algorithm 1, taken from [24]. The algorithm starts by initializing its internal parameters (Line ()). At the iteration, the algorithm samples number of samples from a multivariate Gaussian distribution with mean and covariance (Line ()). Next, the samples are sorted according to their costs . The weighted average of top number of solutions is computed to find the mean vector for the next iteration (Line ()), which moves the mean of the distribution towards samples with lower costs. Next, algorithm updates the covariance matrix with the help of the evolution path variables which are (Line ()) and (Line ()), which ensures that the adaptation steps are conjugate directions. The interested reader is referred to [24] for the full derivation of the algorithm and the intuition for updating the path parameters.

4.2. CMA-MV Algorithm

Unfortunately, CMA algorithm is only applicable to continuous optimization problems; hence we cannot use it to solve the missile launch condition setting problem given in (3), since the allocation of missiles is determined by the integer variables.

To overcome this issue, we develop a novel algorithm named Covariance Matrix Adaptation with Mixed Variables (CMA-MV), which extends the classical CMA algorithm to work on nonlinear mixed integer optimization problems. The generic nonlinear mixed integer programming problem is of the form: where is the set of integers. The special case we are interested in is the problem where the discrete variable is a binary vector; that is, . This is also the case for the missile allocation problem in (3), where refers to missile being launched.

The main idea behind the CMA-MV algorithm is to define and update two probability distributions for sampling continuous and discrete variables. For the continuous variables , we use a multivariate normal distribution and we use the exact same procedure followed in the CMA algorithm (Algorithm 1) to update the mean and covariance of the distribution. For the discrete variables, we use a multivariate Bernoulli distribution and update the mean and covariance of this distribution based on the costs of sampled variables.

However, sampling from the multivariate Bernoulli distribution is not as straightforward as sampling from a multivariate normal distribution. We use the method described in [25] for this purpose. The pseudocode for the sampling process is given in Algorithm 2. The algorithm takes the given mean vector and the covariance matrix for and computes a corresponding multivariate distribution with mean and covariance by solving the equations given on Lines () and (). In these equations is the cumulative distribution of a univariate normal variable with zero mean and unit variance, , where is the cumulative distribution of a bivariate normal variable with mean and correlation . After solving these equations using numerical techniques, we sample the normal variable in Line (). Then we loop through the components of the sample and set if the components are positive and set otherwise. It can be shown that the multivariate sample generated via this fashion comes from a distribution with first and second moments, and , respectively.

The pseudocode for the complete CMA-MV algorithm is given in Algorithm 3. In the iteration, algorithm fixes the value of the discrete variables and hence recovers the function . Note that is a function of a continuous variable; hence we can apply the CMA algorithm (Algorithm 1) to obtain a solution (Line ()). Next, we fix the value of the continuous variable to to recover the function . Then we sample solution candidates from the Bernoulli distribution with mean and using the sampling algorithm given in Algorithm 2 (Line ()). Then we sort the solution from the lowest cost to highest cost (Line ()). Next, we use a weighted average of the low cost solutions to compute the updated mean of the Bernoulli distribution (Line ()). Similarly, we use the weighted sample covariance estimate of the low cost solution candidates to compute the updated covariance matrix . After each update, the distribution puts more mass on low cost solution candidates and hence with each iteration probability of sampling the optimal solution increases.