Optimization of Allocation and Launch Conditions of Multiple Missiles for Three-Dimensional Collaborative Interception of Ballistic Targets

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.

[ProQuest: [...] denotes non US-ASCII text; see PDF]

Academic Editor:Mahmut Reyhanoglu

Aerospace Research Center, Department of Aeronautical Engineering, Istanbul Technical University, 34469 Istanbul, Turkey

Received 12 April 2016; Accepted 13 July 2016

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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; ( 1 ) probability to hit increases due to increased number of missiles targeted at the hostile; ( 2 ) 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 _{x T} = [ _{p T} , _{p T} ] ∈ ^{R 6} denote the state vector of the ballistic target and let ^{x M i} = [ ^{p M i} , ^{p M i} ] ∈ ^{R 6} , i = 1,2 , ... , _{N S A M} , denote the state vectors of the SAMs, where p denotes the position, p denotes the velocity, and _{N S A M} 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 ^{Γ 0} = { ^{γ i} ( 0 ) , i = 1 , ... , _{N S A M} } and ^{Ψ 0} = { ^{ψ i} ( 0 ) , i = 1 , ... , _{N S A M} } denote the initial launching pitch and heading angles for each SAM. Let _{t f} be the terminal time for the scenario. Define the miss distance as [figure omitted; refer to PDF] 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 [ 0 , _{t f} ] 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 ^{Γ 0} , ^{Ψ 0} : [figure omitted; refer to PDF] where _{f T} and _{f M} 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 z = { _{z i} ∈ { 0,1 } , i = 1 , ... , _{N S A M} } denote the binary decision vector, where _{z i} = 1 means that SAM i is launched for intercepting the ballistic target and _{z i} = 0 means that SAM i is not launched and will stay in its initial position throughout the mission. The modified optimization problem is formulated as [figure omitted; refer to PDF] where r > 0 is a user defined parameter that reflects the cost for launching a missile. Larger values of r favor launching lesser number of missiles. Since the optimization problem in (3) contains both continuous variables ^{Γ 0} , ^{Ψ 0} and discrete variables z and since the dynamics _{f M} and _{f T} 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] [figure omitted; refer to PDF] where x is the downrange displacement of the SAM, y is the cross-range displacement of the SAM, z is the altitude of the SAM, g is the gravity acceleration, γ is the flight path angle, _{V m} is SAM velocity, and _{a p i t c h} , _{a y a w} 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.

[figure omitted; refer to PDF]

The approximate mathematical model of the target in the reentry phase is given in [20] [figure omitted; refer to PDF] where _{F d} is the drag force, _{V T x} , _{V T y} , and _{V T z} are the velocity components, m is the target missile mass in the reentry phase, g is the gravity acceleration, _{γ T} is the reentry angle, and _{ψ T} 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.

[figure omitted; refer to PDF]

We assume that the radar estimates the initial state vector _{x ^ 0 ; R} without any error and send the initialization data to the interceptors without any delay [figure omitted; refer to PDF]

The geometric relation between the i t h interceptor and the target is given in (7). In this work, it is assumed that _{γ 0 , T} and _{a 0 , T} are estimated perfectly by the radar [figure omitted; refer to PDF]

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 omitted; refer to PDF]

Figure 3: 3D engagement geometry.

[figure omitted; refer to PDF]

The state vector of the i t h missile according to the target at time t is given in (9) where i = 1,2 , ... , n . In the information sharing mode, each missile could transmit and receive the estimated target data without any loss and delay [figure omitted; refer to PDF]

Relative kinematics between the interceptors and the ballistic missile in _{S x z} , _{S x y} , and _{S y z} planes are defined by using [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

The relative kinematics equations can be described in the discrete-time by using [figure omitted; refer to PDF] where _{x k} is state vector of the i t h missile at time _{t k} and _{f k - 1} 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 _{λ k ; i} . This measurement has a zero-mean Gaussian noise with standard deviation _{σ λ i} . LOS angle measurements are performed by each interceptor missile separately; therefore E ( _{v k ; i} , _{v k ; j} ) = 0 , ∀ i ≠ j . The i t h missile LOS measurement is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] 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): [figure omitted; refer to PDF]

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

The LOS angle of the missile j measured by the missile i is calculated by using [figure omitted; refer to PDF]

In addition, the range-to-go of the j t h missile measured by the i t h missile is obtained by using [figure omitted; refer to PDF]

In the information sharing mode, the measurement model is given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the transition matrix associated with the relative kinematics, T is sampling time, and I is the identity matrix with appropriate dimensions. _{Q k - 1} is the covariance matrix of the equivalent discrete process noise and it is calculated as shown in [figure omitted; refer to PDF] _{F k - 1} is the Jacobian matrix associated with the nonlinear relative kinematics [figure omitted; refer to PDF] The state estimation is performed by using [figure omitted; refer to PDF] where _{x ^ k | k - 1} is the predicted state vector and K is the Kalman gain as shown in [figure omitted; refer to PDF] Here, _{H k} is the measurement Jacobian matrix and _{R k} is the measurement noise covariance matrix. The covariance matrix is updated as shown in [figure omitted; refer to PDF]

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 _{V c} and LOS rate λ [figure omitted; refer to PDF] where _{n c} is the acceleration command perpendicular to the LOS, _{V c} is closing velocity, and N is navigation ratio which is generally between 3 and 5 [20]. In this 3D interception study, TPN algorithm is applied for _{S x z} , _{S x y} , and _{S y z} planes separately [21]. Geometry of relative kinematics for each different plane is displayed in Figures 4, 5, and 6.

Figure 4: Relative kinematics on _{S x z} plane.

[figure omitted; refer to PDF]

Figure 5: Relative kinematics on _{S x y} plane.

[figure omitted; refer to PDF]

Figure 6: Relative kinematics on _{S y z} plane.

[figure omitted; refer to PDF]

Acceleration commands in _{S x z} , _{S x y} , and _{S y z} are obtained as shown in [figure omitted; refer to PDF] The acceleration components of the interceptor in the x -, y -, and z -axis ( _{a m x} , _{a m y} , _{a m z} ) can be obtained from (30) by using the trigonometric relations: [figure omitted; refer to PDF] Before applying the control commands to the interceptor, vertical and horizontal components _{a p i t c h} and _{a y a w} should be calculated. Here, _{a p i t c h} is in the pitch plane and perpendicular to the velocity vector of the interceptor and _{a y a w} is perpendicular to both velocity vector and vertical acceleration vector. For TPN, these acceleration components are calculated using [figure omitted; refer to PDF]

4. Optimization Algorithm

4.1. CMA Algorithm

Consider the general form unconstrained optimization problem: [figure omitted; refer to PDF] It is well known that when f 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 f 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 f , many global methods treat f as a black box function, and hence the solution is found entirely by examining the input-output pairs ( x , f ( x ) ) .

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 ^{R n} and sample candidate solutions ( _{x i} , f ( _{x i} ) ) 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 ^{x [low *]} 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 ( 1 )). At the k t h iteration, the algorithm samples λ number of samples from a multivariate Gaussian distribution with mean _{m k} and covariance _{C k} (Line ( 4 )). Next, the samples _{x i} are sorted according to their costs _{f i} . The weighted average of top μ number of solutions is computed to find the mean vector _{m k + 1} for the next iteration (Line ( 6 )), 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 _{p σ} (Line ( 7 )) and _{p c} (Line ( 13 )), 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.

Algorithm 1: Covariance Matrix Adaptation (CMA).

Input : Objective Function f : ^{R n} [arrow right] R , Number of Samples per Iteration λ , Number of Iterations _{n i t e r} , Weights _{w 1} , ... , _{w μ}

Output : Approximate optimal solution ^{x [low *]}

// I nitialization

(1) _{m 0} , _{σ 0} , _{C 0} [arrow left] I , _{p σ} , _{p c} , k [arrow left] 0 , _{μ w} [arrow left] 1 / ^{∑ j = 1 μ w j 2} , _{c 1} [arrow left] 2 / ^{n 2} , _{c μ} [arrow left] _{μ w} / ^{n 2} , ^{c σ - 1} [arrow left] n / 3 , ^{c c - 1} [arrow left] n / 4 , α [arrow left] 1.5 μ [arrow left] λ / 2

(2) while k < _{n i t e r} do

(3) for i in ( 1 , ... , λ ) do

// Sample Candidate Solutions

(4) _{x i} ~ N ( m , ^{σ k 2} _{C k} ) , _{f i} [arrow left] f ( _{x i} )

// Sort the Candidate Solutions Based on Their Cost

(5) _{x 1 , ... , λ} [arrow left] _{x t ( 1 )} , ... , _{x t ( λ )} , such that _{f t ( 1 )} <= [...] <= _{f t ( λ )}

// Move the mean to low cost solutions

(6) _{m k + 1} [arrow left] _{m k} + ^{∑ i = 1 μ} _{w i} ( _{x i} - _{m k} )

// Update Evolution Path Variables

(7) _{p σ} [arrow left] ( 1 - _{c σ} ) _{p σ} + 1 - ( 1 - _{c σ} ^{) 2} _{μ w} ^{C k - 1 / 2} ( ( _{m k + 1} - _{m k} ) / _{σ k} )

(8) _{σ k + 1} [arrow left] _{σ k} × exp [...] ( _{c σ} ( ( _{p σ} ) / E ( N ( 0 , I ) ) - 1 ) )

// Update The Covariance Matrix

(9) if ( _{p σ} ) < α ( n ) then

(10) _{d k} [arrow left] 1

(11) else

(12) _{d k} [arrow left] 0

(13) _{c s} [arrow left] ( 1 - ^{d k 2} ) _{c 1 c c} ( 2 - _{c c} )

_{p c} [arrow left] ( 1 - _{c c} ) _{p c} + _{d k} ( 1 - _{c c} ^{) 2} _{μ w} ^{C k - 1 / 2} ( ( _{m k + 1} - _{m k} ) / _{σ k} )

(14) _{C k + 1} [arrow left] ( 1 - _{c 1} - _{c μ} + _{c s} ) _{C k} + _{c 1} ^{p c [inverted perpendicular]} _{p c} + _{c μ} ^{∑ i = 1 μ} _{w i} ( ( _{x i} - _{m k} ) / _{σ k} ) ^{( ( _{x i} - _{m k} ) / _{σ k} ) [inverted perpendicular]}

(15) k [arrow left] k + 1

// After the algorithm stops, output the best sample

(16) ^{x [low *]} [arrow left] _{x t ( 1 )}

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: [figure omitted; refer to PDF] where Z is the set of integers. The special case we are interested in is the problem where the discrete variable z is a binary vector; that is, z ∈ { 0,1 ^{} d} . This is also the case for the missile allocation problem in (3), where _{z i} = 1 refers to missile i 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 x , 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 ^{m [variant prime]} and the covariance matrix for ^{C [variant prime]} and computes a corresponding multivariate distribution with mean γ and covariance Λ by solving the equations given on Lines ( 2 ) and ( 4 ). In these equations Φ is the cumulative distribution of a univariate normal variable with zero mean and unit variance, Ψ ( x , y , z ) = _{Φ 2} ( x , y , z ) - Φ ( x ) Φ ( y ) , where _{Φ 2} ( x , y , z ) is the cumulative distribution of a bivariate normal variable with mean [ x , y ] and correlation z . After solving these equations using numerical techniques, we sample the normal variable in Line ( 8 ). Then we loop through the components of the sample and set _{z l} = 1 if the components are positive and set _{z l} = 0 otherwise. It can be shown that the multivariate sample generated via this fashion comes from a distribution with first and second moments, ^{m [variant prime]} and ^{C [variant prime]} , respectively.

Algorithm 2: Sample from a multivariate Bernoulli distribution.

Input : Mean ^{m [variant prime]} ∈ [ 0,1 ^{] d} , Covariance ^{C [variant prime]}

Output : Sample z ∈ ^{( 0,1 ) d}

// C ompute the corresponding multivariate Normal distrubtion

(1) for i in ( 1 , ... , d ) do

(2) _{γ i} [arrow left] ^{Φ - 1} ( ^{m i [variant prime]} )

(3) for j in ( 1 , ... , d ) do

(4) if i ≠ j then

(5) _{Λ i j} [arrow left] Solve ^{C i j [variant prime]} - Ψ ( _{γ i} , _{γ j} , _{Λ i j} ) = 0

(6) else

(7) _{Λ i j} = 1

// Sample from the corresponding multivariate Normal distrubtion and transform the results

(8) u ~ N ( γ , Λ )

(9) for l in ( 1 , ... , d ) do

(10) if _{u l} > 0 then

(11) _{z l} [arrow left] 1

(12) else

(13) _{z l} [arrow left] 0

The pseudocode for the complete CMA-MV algorithm is given in Algorithm 3. In the k t h iteration, algorithm fixes the value of the discrete variables and hence recovers the function f ~ ( x ) = f ( x , ^{z k} ) . Note that f ~ is a function of a continuous variable; hence we can apply the CMA algorithm (Algorithm 1) to obtain a solution (Line ( 3 )). Next, we fix the value of the continuous variable to ^{x k + 1} to recover the function f ^ ( z ) = f ( ^{x k + 1} , z ) . Then we sample ^{λ [variant prime]} solution candidates from the Bernoulli distribution with mean ^{m k [variant prime]} and ^{C k [variant prime]} using the sampling algorithm given in Algorithm 2 (Line ( 5 )). Then we sort the solution from the lowest cost to highest cost (Line ( 6 )). Next, we use a weighted average of the low cost solutions to compute the updated mean ^{m k + 1} of the Bernoulli distribution (Line ( 7 )). Similarly, we use the weighted sample covariance estimate of the low cost solution candidates to compute the updated covariance matrix ^{C k + 1 [variant prime]} . 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.

Algorithm 3: Covariance Matrix Adaptation with Mixed Variables (CMA-MV).

Input : Objective Function f : ^{R n} × ^{( 0,1 ) d} [arrow right] R ,

Number of Continuous Samples per Iteration λ ,

Number of Discrete Samples per Iteration λ ,

Number of Iterations _{n i t e r} , Weights for continuous samples _{w 1} , ... , _{w μ} , Weights for

discrete samples ^{w 1 [variant prime]} , ... , ^{w μ [variant prime] [variant prime]} , Initial guess ^{x 0} , ^{z 0}

Output : Approximate optimal solution ^{x [low *]} , ^{z [low *]}

// I nitialization

(1) _{m 0} , _{σ 0} , _{C 0} [arrow left] I , _{p σ} , _{p c} , k [arrow left] 0 , _{μ w} [arrow left] 1 / ^{∑ j = 1 μ w j 2}

(2) while k < _{n i t e r} do

// Fix the discrete variables to ^{z k} and use CMA to solve for x

(3) ^{x k + 1} [arrow left] Call CMA (Algorithm 1) with f ~ ( x ) = f ( x , ^{z k} )

// Fix the continous variables and sample from the multivariate Bernoulli distrubution

(4) for i in ( 1 , ... , ^{λ [variant prime]} ) do

// Sample

(5) _{z i} ~ B ( ^{m k [variant prime]} , ^{C k [variant prime]} ) using Algorithm 2 _{f ^ i} [arrow left] f ( ^{x k + 1} , _{z i} )

// Sort the candidate Solutions Based on Their Cost

(6) _{z 1 , ... , ^{λ [variant prime]}} [arrow left] _{z t ( 1 )} , ... , _{z t ( ^{λ [variant prime]} )} , such that _{f ^ t ( 1 )} <= [...] <= _{f ^ t ( ^{λ [variant prime]} )}

// Move the mean to low cost solutions

(7) ^{m k + 1 [variant prime]} [arrow left] ^{∑ i = 1 μ [variant prime] w i [variant prime]} _{z i}

// Update The Covariance Matrix

(8) ^{C k + 1 [variant prime]} [arrow left] 1 ^{λ [variant prime]} - 1 ^{∑ i = 1 μ [variant prime] w i [variant prime]} ( _{z i} - ^{m k + 1 [variant prime]} ) ( _{z i} - ^{m k + 1 [variant prime] ) [inverted perpendicular]}

(9) ^{z k + 1} [arrow left] ^{z t k} ( 1 )

(10) k [arrow left] k + 1

// After the algorithm stops, output the best sample

(11) ^{x [low *]} [arrow left] ^{x t ( 1 ) k} , ^{z [low *]} [arrow left] ^{z t ( 1 ) k}

5. Simulation Results

In this section we fuse our optimization algorithm (Algorithm 3) with the control and estimation methods given in Section 3 to create an integrated solution to multiple missile allocation and control for ballistic target interception. We first give detailed results for two specific missions in order to give a better understanding of how the algorithm works, and then we demonstrate the effectiveness of the algorithm by comparing its performance to Heuristic and noncollaborative methods in Monte-Carlo simulations.

In all experiments, we use the following parameters for SAM defense system:

(i) The number of SAMs _{N S A M} is 5 , and they are arranged in two parallel lines with the back line containing 3 SAMs and the front line containing 2 SAMs. The arrangement can be seen in the upper left corners of Figures 7 and 9. Radar is placed in front of the front line.

(ii) SAM velocity is set to Mach = 3.5 .

(iii): Maximum number of iterations for CMA-MV is set to 50 . The number of samples is set to 100 for both continuous and discrete variables. The rest of the parameters are tuned manually.

Figure 7: Interception of a high altitude-low velocity ballistic target. Red line depicts the trajectory of the ballistic target. Algorithm chooses only 1 missile for interception.

[figure omitted; refer to PDF]

5.1. Results for a High Altitude-Low Velocity Target

First we examine a mission where the ballistic target has relatively low kinetic energy. The ballistic target's initial conditions are set to 80000 meters of altitude and speed equivalent to Mach number 5 . This is a less challenging scenario, since the ballistic threat has relatively longer time till it hits the ground, giving the enough time for the filters of the SAM defense system to converge. Resulting trajectory of the target and the launched missile is shown in Figure 7. The optimization algorithm also recognizes that filters have enough time to converge in this case and launches only a single missile. The missile intercepts the ballistic threat with a miss distance less than 1 meter.

Figures 8(a), 8(b), and 8(c) show the estimation performance of the filter of the missile, for range-to-go, line of sight, and target acceleration estimations. It can be seen that filters converged rapidly in the terminal phase of the mission. These plots justify the decision of the algorithm to launch only a single missile; in this case algorithm recognized that a single filter would yield sufficient performance and did not choose to allocate more missiles, in order to keep the cost as close to minimum as possible. Also note that no collaborative filtering is performed in this mission, since only a single missile is launched.

Figure 8: Filter performance versus time in terminal phase for intercepting high altitude-low velocity target.

(a) Norm of the range-to-go estimation error versus time in terminal phase

[figure omitted; refer to PDF]

(b) Norm of the line of sight estimation error versus time in terminal phase

[figure omitted; refer to PDF]

(c) Norm of the target acceleration estimation error versus time in terminal phase

[figure omitted; refer to PDF]

Figure 9: Interception of a low altitude-high velocity ballistic target. Red line depicts the trajectory of the ballistic target. Algorithm launches 3 missiles to intercept the target.

[figure omitted; refer to PDF]

5.2. Results for a Low Altitude-High Velocity Target

To complement the results of the previous subsection, now we look at a mission that corresponds to high kinetic energy target. For this simulation the ballistic target's initial conditions are set to 60000 meters of altitude and speed equivalent to Mach number 7 . This scenario is much more challenging, since the target's established time of impact is much shorter. Resulting trajectory of the target and the launched missile is shown in Figure 9. In this case it is seen that the algorithm launches 3 collaborative missiles to intercept the target. The interception is achieved with a miss distance of approximately 1 meter.

Figures 10(a), 10(b), and 10(c) show the estimation performance of the filter of the missiles, for range-to-go, line of sight, and target acceleration estimations averaged over the 3 launched missiles. For comparison, performance of the individual filters is also plotted in these figures. These individual filter performances correspond to the case where missiles do not cooperate; hence no information is shared between them. Examining these plots gives us a good insight on algorithm's decision to launch 3 missiles. It is seen that individual filters did not have enough time to converge for this high kinetic energy target; hence launching a single or even two missiles would result in high miss distances. Algorithm made the necessary trade-off analysis and found out that launching 3 collaborative missiles would generate enough information flow for estimators to converge.

Figure 10: Filter performance versus time in terminal phase for intercepting low altitude-high velocity target.

(a) Average norm of the range-to-go estimation error versus time in terminal phase for individual and collaborative estimators

[figure omitted; refer to PDF]

(b) Average norm of the line of sight estimation error versus time in terminal phase for individual and collaborative estimators

[figure omitted; refer to PDF]

(c) Average norm of the target acceleration estimation error versus time in terminal phase for individual and collaborative estimators

[figure omitted; refer to PDF]

5.3. Monte-Carlo Results for Multiple Scenarios

The previous simulation results demonstrated that the algorithm yields sound decisions on selected scenarios. However, in order to truly assess the performance of the algorithm, a wide range of initial conditions that corresponds to different ballistic threat should be analyzed. Also we need to compare the performance of the algorithm to alternative methods. For this purpose we conducted a Monte-Carlo test over 100 randomly sampled initial conditions for the ballistic target. The initial altitude of the target was sampled in the interval [ 40000,80000 ] meters and the speed was sampled in the Mach number interval [ 5,8 ] . The following alternative methodologies are compared:

(i) Heuristic Collaborative Interception . In this simple algorithm no launch condition or missile allocation optimization is conducted. This method always launches the same number of missiles that are closest to the ballistic target at the beginning of the simulation. Launch conditions of the missiles are always set to 0 heading and 90 degrees of pitch angle. Missiles use collaborative filtering for interception.

(ii) Optimized Noncollaborative Interception . In this methodology optimization algorithm CMA-MV is used for optimizing the launching conditions and the allocation of the missiles. However missiles do not run collaborative filtering algorithms on-board.

(iii): Optimized Collaborative Interception . This is the approach developed in this paper. The CMA-MV algorithm (Algorithm 3) is used for optimization of launch conditions and missile allocation and the missiles run collaborative filtering algorithms.

Table 1 depicts the results of the Monte-Carlo analysis. We see that Heuristic method's performance gets better as the number of missiles used by the method increases. This is expected, since increased number of missiles translates to improved estimation performance. However, even using 3 missiles for all conditions does not reduce the average miss distance substantially. This is due to fact that Heuristic method does not optimize the launch conditions; hence the missile autopilots do not have enough time to restore the missiles into the desired trajectories in the terminal phase. On the other hand, optimized noncollaborative method yields substantially lower miss distances than Heuristic methods that use 1 or 2 missiles, while launching only 1.05 missiles on average. This is because the optimized noncollaborative method optimizes the launch conditions for the missiles, which leads to improved interception performance. However, this method is outperformed by Heuristic method that uses 3 missiles, because optimized noncollaborative method does not utilize collaborative filters; hence the algorithm can not take advantage of improved estimation performance gained by launching multiple missiles against high kinetic energy ballistic targets.

Table 1: Average number of missiles and average miss distance obtained by different methods, averaged over 100 random initial conditions of the ballistic target.

	Average number of missiles launched	Average miss distance (m)
Heuristic method 1 missile	1	5020.78 ± 600.21
Heuristic method 2 missiles collaborative	2	490.22 ± 41.08
Heuristic method 3 missiles collaborative	3	60.55 ± 11.2
Optimized noncollaborative	1.05 ± 0.23	320.55 ± 50.77
Optimized collaborative	1.56 ± 0.34	1.43 ± 0.22

Finally, we see that the approach developed in this paper, the optimized collaborative method, outperforms the compared approaches in terms of both resource management efficiency and miss distance. This is because, unlike the compared approaches, the developed method optimizes the launch conditions and missile allocation simultaneously, and hence it is able to assess the right trade-off between the number of missiles launched and attainable miss distance.

6. Conclusions and Future Work

In this work we have developed a novel probabilistic search algorithm for allocation and launch condition optimization of multiple missiles intercepting a ballistic target. Through simulation studies, we have verified that the algorithm makes sound decisions and yields more efficient resource management and lower miss distance compared to alternative approaches.

Future works involve generalization of the problem into intercepting multiple ballistic targets, developing alternative cooperative estimation/control algorithms, and investigation of theoretical properties of the CMA-MV algorithm.

References

[1] B.-Z. Naveh, A. Lorber, "Theater ballistic missile defense," Progress in Astronautics and Aeronautics , vol. 192, pp. 1-397, 2001.

[2] G. M. Siouris Missile Guidance and Control Systems , Springer, New York, NY, USA, 2004.

[3] R. H. Chen, J. L. Speyery, "Terminal and boost phase intercept of ballistic missile defense," in Proceedings of the AIAA Guidance, Navigation and Control Conference and Exhibit, Honolulu, Hawaii, USA, August 2008.

[4] T.-K. Wang, L.-C. Fu, "A guidance strategy for multi-player pursuit and evasion game in maneuvering target interception," in Proceedings of the 9th Asian Control Conference (ASCC '13), June 2013.

[5] I.-S. Jeon, J.-I. Lee, M.-J. Tahk, "Homing guidance law for cooperative attack of multiple missiles," Journal of Guidance, Control, and Dynamics , vol. 33, no. 1, pp. 275-280, 2010.

[6] E. Daughtery, Z. Qu, "Optimal design of cooperative guidance law for simultaneous strike," in Proceedings of the 53rd IEEE Annual Conference on Decision and Control (CDC '14), pp. 988-993, Los Angeles, Calif, USA, December 2014.

[7] V. Shaferman, Y. Oshman, "Cooperative interception in a multi-missile engagement," in Proceedings of the AIAA Guidance, Navigation, and Control Conference, Chicago, Ill, USA, August 2009.

[8] V. Shaferman, T. Shima, "Cooperative multiple-model adaptive guidance for an aircraft defending missile," Journal of Guidance, Control, and Dynamics , vol. 33, no. 6, pp. 1801-1813, 2010.

[9] Y. Liu, N. Qi, J. Shan, "Cooperative interception with double-line-of-sight measuring," in Proceedings of the AIAA Guidance, Navigation and Control Conference, Boston, Mass, USA, August 2013.

[10] V. Shaferman, T. Shima, "Cooperative optimal guidance laws for imposing a relative intercept angle," Journal of Guidance, Control, and Dynamics , vol. 38, no. 8, pp. 1395-1408, 2015.

[11] L. Wang, H. Fenghua, Y. Yu, "Guidance law design for two flight vehicles cooperative interception," in Proceedings of the AIAA Guidance, Navigation and Control Conference, Kissimmee, Fla, USA, January 2015.

[12] Z.-J. Lee, S.-F. Su, C.-Y. Lee, "Efficiently solving general weapon-target assignment problem by genetic algorithms with greedy eugenics," IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics , vol. 33, no. 1, pp. 113-121, 2003.

[13] P. Teng, H. Lv, J. Huang, L. Sun, "Improved particle swarm optimization algorithm and its application in coordinated air combat missile-target assignment," in Proceedings of the 7th World Congress on Intelligent Control and Automation (WCICA '08), pp. 2833-2837, IEEE, Chongqing, China, June 2008.

[14] J. Wang, X.-G. Gao, Y.-W. Zhu, H. Wang, "A solving algorithm for target assignment optimization model based on ACO," in Proceedings of the 6th International Conference on Natural Computation (ICNC '10), vol. 7, pp. 3753-3757, IEEE, Yantai, China, August 2010.

[15] J. Wang, Y.-W. Zhu, "A solving algorithm for target assignment optimization model based on SA," in Proceedings of the International Conference on Artificial Intelligence and Computational Intelligence (AICI '10), pp. 489-493, Sanya, China, October 2010.

[16] M. Wei, G. Chen, K. Pham, E. Blasch, Y. Wu, "Game theoretic target assignment approach in ballistic missile defense," in Proceedings of the SPIE Defense and Security Symposium, International Society for Optics and Photonics, 2008.

[17] B. Xin, J. Chen, Z. Peng, L. Dou, J. Zhang, "An efficient rule-based constructive heuristic to solve dynamic weapon-target assignment problem," IEEE Transactions on Systems, Man, and Cybernetics Part A:Systems and Humans , vol. 41, no. 3, pp. 598-606, 2011.

[18] N. Hansen, A. Ostermeier, "Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation," in Proceedings of the IEEE International Conference on Evolutionary Computation (ICEC '96), pp. 312-317, May 1996.

[19] Z. B. Tatas Modeling and autopilot design for a SCUD Type ballistic missile [M.S. thesis] , Department of Electrical Engineering, Hacettepe University, Ankara, Turkey, 2006.

[20] P. Zarchan Tactical and Strategic Missile Guidance , vol. 199, of Progress in Aeronautics and Astronautics, AIAA, 2002., 4th.

[21] I. Moran, D. T. Altilar, "Three plane approach for 3D true proportional navigation," in Proceedings of the AIAA Guidance, Navigation and Control Conference, San Francisco, Calif, USA, August 2005.

[22] E. K. P. Chong, H. Stanislaw An Introduction to Optimization , vol. 76, John Wiley & Sons, New York, NY, USA, 2013.

[23] N. Hansen, A. Auger, R. Ros, S. Finck, P. Posík, "Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009," in Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation (GECCO '10), ACM, July 2010.

[24] N. Hansen, S. D. Müller, P. Koumoutsakos, "Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)," Evolutionary Computation , vol. 11, no. 1, pp. 1-18, 2003.

[25] J. H. Macke, P. Berens, A. S. Ecker, A. S. Tolias, M. Bethge, "Generating spike trains with specified correlation coefficients," Neural Computation , vol. 21, no. 2, pp. 397-423, 2009.