# A Modified Pareto Ant Colony Optimization Approach to Solve Biobjective Weapon-Target Assignment Problem

Authors: You Li; Yingxin Kou; Zhanwu Li; An Xu; Yizhe Chang
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

The weapon-target assignment (WTA) problem, known as an NP-complete problem, aims at seeking a proper assignment of weapons to targets. The biobjective WTA (BOWTA) optimization model which maximizes the expected damage of the enemy and minimizes the cost of missiles is designed in this paper. A modified Pareto ant colony optimization (MPACO) algorithm is used to solve the BOWTA problem. In order to avoid defects in traditional optimization algorithms and obtain a set of Pareto solutions efficiently, MPACO algorithm based on new designed operators is proposed, including a dynamic heuristic information calculation approach, an improved movement probability rule, a dynamic evaporation rate strategy, a global updating rule of pheromone, and a boundary symmetric mutation strategy. In order to simulate real air combat, the pilot operation factor is introduced into the BOWTA model. Finally, we apply the MPACO algorithm and other algorithms to the model and compare the data. Simulation results show that the proposed algorithm is successfully applied in the field of WTA which improves the performance of the traditional P-ACO algorithm effectively and produces better solutions than the two well-known multiobjective optimization algorithms NSGA-II and SPEA-II.


---

## Body

## 1. Introduction

Among air power, land power, and sea power, air power plays the most decisive role in the outcome of any war [1]. With the development of military air combat, the problem of WTA has received considerable attention [2–6]. Because the goal of air power is to maximize the military effectiveness, it becomes imperative to find a proper assignment of weapons to targets. From the perspective of the quantity of objective functions, the WTA problem can be classified into single-objective WTA problem and multiple-objective WTA problem. When taking the time factor into account, the WTA problem can be divided into dynamic WTA problem and static WTA problem. The current research status of various WTA is summarized in Table 1. In this paper, we mainly focus on the static biobjective WTA problem which aims at finding a proper static assignment of maximizing the expected damage of the enemy and minimizing the cost of missiles.Table 1
Summary of variant metaheuristic algorithms and implementation of various WTA.
Researchers
Year
Metaheuristic algorithm
Implementation (WTA)
Lee et al. [40]
2002
IS + ACO
Static single-objective
Lee et al. [41]
2002
GA
Static single-objective
Z.-J. Lee and W.-L. Lee [15]
2003
GA + ACO
Static single-objective
Galati and Simaan [14]
2007
Tabu
Dynamic single-objective
Lee [13]
2010
VLSN
Static single-objective
Xin et al. [42]
2010
VP + Tabu
Dynamic single-objective
Li and Dong [16]
2010
DPSO + SA
Dynamic single-objective
Xin et al. [43]
2011
Rule-based heuristic
Dynamic multiobjective
Chen et al. [44]
2010
SA
Static single-objective
Bogdanowicz et al. [10]
2013
GA
Static single-objective
Fei et al. [12]
2012
Auction algorithm
Static single-objective
Liu et al. [18]
2013
MOPSO
Static multiobjective
Zhang et al. [19]
2014
MOEA/D
Static multiobjective
Ahner and Parson [45]
2015
Dynamic programming
Dynamic multiobjective
Li et al. [20]
2015
NSGA-II, MOEA/D
Static multiobjective
Dirik et al. [46]
2015
MILP
Dynamic multiobjective
Liang and Kang [47]
2016
CSA
Static single-objective
Li et al. [48]
2016
MDE
Dynamic multiobjectiveThe WTA problem belongs to NP-complete problems (NP: nondeterministic polynomial), so the increment of weapons and targets will result in an explosive growth of the searching space. Although the enumeration method ensures a global solution, its computational requirement grows exponentially as the scale of WTA problem grows. Therefore, seeking better weapons for targets distribution scheme and improving the operational efficiency especially for the large-scale WTA problem are still an important issue in the command and control problem (C2).The study of WTA problem can be traced back to the 1950s and 1960s when Manne [7] and Day [8] built the model of WTA problem. Hosein and Athans [9] classified the WTA problem into two classes: the single-objective WTA problem and the multiple-objective WTA problem. Genetic algorithm [10], ACO algorithm [11], auction algorithm [12], VLSN algorithm [13], Tabu search [14], and other hybrid algorithms [15–17] have been used to optimize single-objective WTA model by many scholars. In contrast to single-objective WTA, multiple-objective optimization can take different criterions into consideration and is more in line with the real combat decision-making. So, it has aroused wide attention from scholars. Liu et al. [18] proposed an improved multiobjective particle swarm optimization (MOPSO) algorithm to solve multiple-objective WTA problem and apply it to a simple example including 7 platforms and 10 targets. However, they did not apply the proposed algorithm to different scales WTA problems for contrastive study. Zhang et al. [19] designed a WTA mathematic model and proposed a decomposition-based evolutionary multiobjective optimization algorithm. But the algorithm has not been tested on large-scale WTA problem and it has a low convergence speed. Li et al. [20] adopted NSGA-II (domination-based) and MOEA/D (decomposition-based) to solve the multiple-objective WTA problem and carried out tests on three different BOWTA problems. They only applied the proposed adaptive mechanism to the WTA problems, but they did not verify the behavior of the proposed adaptive mechanism on standard problems.Different from other algorithms, ACO algorithm is a class of reactive search optimization (RSO) methods adopting the principle of “learning while optimizing” [21, 22]. Since it was introduced in 1992, many variant ACO algorithms have been presented, including ant colony system (ACS) [23] and MAX-MIN ant system (MMAS) [24]. Meanwhile, ACO algorithms have been intensively investigated and successfully applied to deal with multiobjective problems such as travelling salesmen problem (TSP) [25], scheduling problem [26], vehicle routing problem [27], portfolio selection problem [28], network optimization problem [29], and some others problems [30–32]. The main applications of ACO algorithms are shown in Figure 1. In 2004, Doerner et al. proposed the P-ACO algorithm [28] which combined traditional ant system with ant colony optimization. Because P-ACO has the ability to balance local search and global search, it will avoid premature convergence during the solution construction phase. The characteristic is especially suitable for solving portfolio selection problem. Therefore, we choose P-ACO as the approach for solving the BOWTA problem. However, to the best of our knowledge, there has not been any research conducted about P-ACO for BOWTA problem.Figure 1
The main applications of ACO algorithm in recent years.Basic multiobjective optimization algorithms have such shortcomings as slow convergence to the Pareto front, no efficient termination criterion, and a lack of theoretical convergence proof. In this paper, we have proposed a modified P-ACO for tackling the above problems. The proposed algorithm is used to seek better compromise solutions between maximizing the expected damage and minimizing the cost. Based on the framework of the traditional P-ACO, the proposed algorithm is devised with several attractive features to enhance the optimization performance, including a dynamic heuristic information calculation approach, an improved movement probability rule, a dynamic evaporation rate strategy, a global updating rule of pheromone, and a boundary symmetric mutation strategy.Dynamic Heuristic Information Calculation Approach. Because heuristic information seriously affects the efficiency of ACO algorithm, the problem-specific heuristic information function is one of the most important problems that should be solve. Hence, we design a dynamic heuristic information calculation approach for providing necessary guidance for an efficient search.Improved Movement Probability Rule. In order to avoid P-ACO plunge into the local optimum which is caused by the stagnation phenomenon, we proposed a dynamic movement probability rule to enhance local exploitation and alleviate the premature convergence.Dynamic Evaporation Rate Strategy. When the algorithm is used to solve large-scale WTA problem, a static pheromone evaporation rate will reduce global searching ability and convergence speed. Therefore, a dynamic evaporation rate strategy is proposed in this paper.Global Updating Rule of Pheromone. The traditional static pheromone updating rule makes the algorithm easy to fall into local optimum. Hence, a global pheromone updating rule is introduced to solve the conflict between expanding search and finding optimal solution.Boundary Symmetric Mutation Strategy. When solving large-scale problem, the distributions of targets are complex. In order to avoid the search being stuck in local optimum and diversify the solutions, a boundary symmetric mutation strategy is proposed to improve the mutation efficiency and the quality of solutions.The rest of this paper is organized as follows. In Section2, a mathematical model of the BOWTA problem and a problem-specific assumption description are presented. Section 3 provides algorithmic descriptions of P-ACO algorithm and presents the application of ACO algorithm for solving WTA problem. Detailed improvements of the proposed algorithm are introduced in Section 4. In Section 5, experimental evaluation is carried out and the analysis of the results is presented. Firstly, in order to validate effectiveness and show performance, the proposed algorithm is applied in a specific instance. Secondly, the proposed algorithm and others, like NSGA-II [33], SPEA-II [34], and traditional P-ACO [28], are tested on several different scales WTA problems. Section 6 provides some conclusion of the paper.


## 2. Problem Formulation

The WTA problem focuses on how to allocate attack units to targets and can be illustrated in Figure2. WTA assumptions and mathematical model are introduced as follows in this section.Figure 2
Illustration of the WTA problem.


### 2.1. Assumption Description

In order to establish a reasonable WTA mathematical model, the following assumptions can be defined.Assumption 1.
We assume that the mathematical model is composed ofF fighters, M missiles, and T targets. The fighters and targets are not necessarily equal in quantity (each fighter is equivalent to one platform which takes different kinds and quantities of missiles).Assumption 2.
Each fighter can use different missiles to attack one target (each missile can only attack one target).Assumption 3.
All allocated missiles cannot exceed the number of missiles resources in military air operation.Assumption 4.
We assume the probability of kill (qij) between the missile (ith unit of M) and the target (jth unit of T).Assumption 5.
If the target is in work area, then a missile can be assigned effectively. On the other hand, the missile is not.


### 2.2. Mathematical Model

The multiobjective optimization is to make the number of weapons to targets reach the maximum and reduce missiles consumption as much as possible. The definitions about the optimization model are shown as follows.Definition 6.
The multiobjective functions of WTA are to maximize the total effectiveness of attack and minimize the cost of missiles. The optimization mathematical model is described as(1)f1=max⁡∑j=1n1-∏i=1m1-ρkqijxij,f2=min⁡∑i=1m∑j=1ncixij,f=min⁡1f1,f2.Definition 7.
q
i
j
∈
0,1 represents the effectiveness of the missile unit. It is known as the probability of kill for each T target (j=1,2,…,n) and M missile (i=1,2,…,m).Definition 8.
The decision tableX=xijM×T can be shown as Table 2.
Here,xij is a Boolean value (xij=0,1) and represents whether i missile is assigned to j target. If the i missile of k fighter is assigned to the j target, xij=1; otherwise, xij is equal to 0.Table 2
The decision table of WTA.
T
1
T
2
T
3
T
4
⋯
T
n
M
1
x
11
x
12
x
13
x
14
⋯
x
1
n
M
2
x
21
x
22
x
23
x
24
⋯
x
2
n
M
3
x
31
x
32
x
33
x
34
⋯
x
3
n
M
4
x
41
x
42
x
43
x
44
⋯
x
4
n
M
5
x
51
x
52
x
53
x
54
⋯
x
5
n
M
6
x
61
x
62
x
63
x
64
⋯
x
6
n
M
7
x
71
x
72
x
73
x
74
⋯
x
7
n
M
8
x
81
x
82
x
83
x
84
⋯
x
8
n
⋮
⋮
⋮
⋮
⋮
x
i
j
⋮
M
m
x
l
1
x
l
2
x
l
3
x
l
4
⋯
x
m
nDefinition 9.
c
i denotes the cost of i missile and is given a constant value in this paper.Definition 10.
In combat, different fighters are equipped with different numbers of missiles, so one assumes that the number of missiles carried by each fighter is not necessarily equal (the number of missiles carried by each plane is no more than 4).Definition 11.
According to the talent, training time, and operation stability of the pilot, all the above factors may affect the attack performance of the missile. One assumes thatρk is pilot operation factor ρk∈0,1.Definition 12.
There are some constraints that the above function variables must satisfy.
Missile Assignment Constraint. For simplicity, we assume that at least one missile unit is assigned to one target. This constraint can be expressed as(2)∑i=1mxij≥1∀j=1,2,…,n.Quantity Assignment Constraint. We assume that one missile can only attack one target. The corresponding constraint is(3)∑j=1nxij≤1∀i=1,2,…,m.Parameters Constraint. In this paper, we will use some parameters like xij. The corresponding constraint is defined as follows:(4)∑i=1m∑j=1nxij≤m.


## 2.1. Assumption Description

In order to establish a reasonable WTA mathematical model, the following assumptions can be defined.Assumption 1.
We assume that the mathematical model is composed ofF fighters, M missiles, and T targets. The fighters and targets are not necessarily equal in quantity (each fighter is equivalent to one platform which takes different kinds and quantities of missiles).Assumption 2.
Each fighter can use different missiles to attack one target (each missile can only attack one target).Assumption 3.
All allocated missiles cannot exceed the number of missiles resources in military air operation.Assumption 4.
We assume the probability of kill (qij) between the missile (ith unit of M) and the target (jth unit of T).Assumption 5.
If the target is in work area, then a missile can be assigned effectively. On the other hand, the missile is not.


## 2.2. Mathematical Model

The multiobjective optimization is to make the number of weapons to targets reach the maximum and reduce missiles consumption as much as possible. The definitions about the optimization model are shown as follows.Definition 6.
The multiobjective functions of WTA are to maximize the total effectiveness of attack and minimize the cost of missiles. The optimization mathematical model is described as(1)f1=max⁡∑j=1n1-∏i=1m1-ρkqijxij,f2=min⁡∑i=1m∑j=1ncixij,f=min⁡1f1,f2.Definition 7.
q
i
j
∈
0,1 represents the effectiveness of the missile unit. It is known as the probability of kill for each T target (j=1,2,…,n) and M missile (i=1,2,…,m).Definition 8.
The decision tableX=xijM×T can be shown as Table 2.
Here,xij is a Boolean value (xij=0,1) and represents whether i missile is assigned to j target. If the i missile of k fighter is assigned to the j target, xij=1; otherwise, xij is equal to 0.Table 2
The decision table of WTA.
T
1
T
2
T
3
T
4
⋯
T
n
M
1
x
11
x
12
x
13
x
14
⋯
x
1
n
M
2
x
21
x
22
x
23
x
24
⋯
x
2
n
M
3
x
31
x
32
x
33
x
34
⋯
x
3
n
M
4
x
41
x
42
x
43
x
44
⋯
x
4
n
M
5
x
51
x
52
x
53
x
54
⋯
x
5
n
M
6
x
61
x
62
x
63
x
64
⋯
x
6
n
M
7
x
71
x
72
x
73
x
74
⋯
x
7
n
M
8
x
81
x
82
x
83
x
84
⋯
x
8
n
⋮
⋮
⋮
⋮
⋮
x
i
j
⋮
M
m
x
l
1
x
l
2
x
l
3
x
l
4
⋯
x
m
nDefinition 9.
c
i denotes the cost of i missile and is given a constant value in this paper.Definition 10.
In combat, different fighters are equipped with different numbers of missiles, so one assumes that the number of missiles carried by each fighter is not necessarily equal (the number of missiles carried by each plane is no more than 4).Definition 11.
According to the talent, training time, and operation stability of the pilot, all the above factors may affect the attack performance of the missile. One assumes thatρk is pilot operation factor ρk∈0,1.Definition 12.
There are some constraints that the above function variables must satisfy.
Missile Assignment Constraint. For simplicity, we assume that at least one missile unit is assigned to one target. This constraint can be expressed as(2)∑i=1mxij≥1∀j=1,2,…,n.Quantity Assignment Constraint. We assume that one missile can only attack one target. The corresponding constraint is(3)∑j=1nxij≤1∀i=1,2,…,m.Parameters Constraint. In this paper, we will use some parameters like xij. The corresponding constraint is defined as follows:(4)∑i=1m∑j=1nxij≤m.


## 3. Ant Colony Optimization Algorithm

### 3.1. Pareto Ant Colony Optimization Algorithm (P-ACO)

The ACO is a heuristic algorithm which has emerged in recent years and has been applied to solve hard combinatorial optimization problems [35]. Scientists have developed ACO by imitating the behaviors of real ant colonies in the process of searching for food. When ants in the colony find a short path for gathering food by frequently travelling between their nest and the food source, they deposit a chemical calledpheromone trail on the path which can be followed by other ants to search for the food source [36].Pareto ant colony optimization (P-ACO) algorithm is one of the multiple-objective ant colony optimization (MOACO) algorithms that was proposed by Doerner et al. [28] to be originally applied to optimize the multiobjective portfolio selection. The main objective of P-ACO algorithm is to find a set of solutions ordered by fronts under the concept of Pareto front. Although it is based on the classical ACS, the update process ofglobal pheromone trail is implemented by using two different ants (the best and the second-best solutions) which are generated by each objective g in current iteration. In P-ACO, τg denotes the pheromone matrices of objective g. During each iteration, wg represents the weight of objective g and is used to combine the pheromone information and heuristic information. In the route, an ant selects the next node s to be visited from node r, and it uses the ACS transition rule like formula (5). ζ is a random number from 0 to 1. ζ0 is a fixed parameter (ζ0∈0,1). When ζ≤ζ0, the available knowledge is exploited for selecting the best option through pheromone trail heuristic information. The formula is shown below: (5)j=argmaxj∈Ω∑g=1Gwg·τijgα·ηijβ,ifζ≤ζ0j^,otherwise,where g=1,2,…,G is a set of objectives, ηij represents the heuristic information of edge eij, and j^ is a probabilistic node. α is the weight coefficient of pheromone information and β is the weight coefficient of heuristic information. α and β belong to R. Ω is the current feasible neighborhood of the ant.Whenζ>ζ0, pij is implemented for finding j^ and the formula is expressed as(6)pij=∑g=1Gwg·τijgα·ηijβ∑u∈Ω∑g=1Gwg·τiugα·ηiuβ,ifj∈Ω0,otherwise.When an ant has travelled an edgeeij in the algorithm,local pheromone update (online pheromone trail update) is implemented for each objective g by the rule(7)τijg=1-ρ·τijg+ρ·τ0,where ρ is the pheromone evaporation rate (ρ∈0,1) and τ0 is the value of initial pheromone trail. In this paper, we use the same method as the literature [37] to set the value of τ0.When each ant has constructed its solution,global pheromone update (offline pheromone trail update) is implemented for each objective g. The update formula is shown below:(8)τijg=1-ρ·τijg+ρ·Δτijg,where Δτijg follows the problem-specific values:(9)Δτijg=10ifedgeeij∈bestsolution,5ifedgeeij∈second-bestsolution,0otherwise.In the process, the nondominated solutions are stored in an external setS.The pseudocode of the P-ACO algorithm can be demonstrated as shown in Algorithm1.Algorithm 1
(0) Pre-processing;
(1) Initialization: generate the number of ants A, initialize pheromone vectors with τ0,
set the values ofα and β;
(2) repeat until Termination conditions has not been met
(3)  for Ant = 1 to A
determine the lifespanΘ of the ant randomly of the interval 1,2,…,N;
createx=0,…,0/∗built empty assignment solution∗/;
determine the objective weightwg for each objective g;
indicates the number of targets to be selectedξ=Θ;
(4)  while  ξ>0  and  ∃i with ηix>0
select an assignmentij using formula (5) and add it to x;
update local pheromone information;
decrementξ;
check feasibilityx;
(5)    if  x is feasible
check efficiencyx;
(6)    if  x is efficient
storex and remove dominated ones;
(7) for each objective g
determine best and second-best solutionΔτijg for each objective g;
update global pheromone information using best and second-best solution;
using formula (8);
.


### 3.2. Application of ACO Algorithm for Solving WTA Problem

At the beginning, each missile has a list of possible target assignments (Figure3). The number of assignments is independent of the conditions which are set at first. Each ant builds a tour (WTA solution) while moving from one to another list until all of them have been visited. Figure 3 gives an example with two ants, four missiles, and six targets during an iteration. M=1,2,…,m is a set of missiles. The size of each list is fixed to the number of targets (T={1,2,…,n}). To construct a WTA solution (a tour), the 1st ant visits Target 2 (missile 1), Target 5 (missile 2), Target 4 (missile 3), and Target 3 (missile 4) by using a transition rule which takes the pheromone trails of visited nodes and the heuristic information into account. On the contrary, the 2nd ant visits Target 4 (missile 1), Target 3 (missile 2), Target 6 (missile 3), and Target 1 (missile 4). If the 1st solution is better than the 2nd solution, then the 1st solution is the best solution at current iteration and the algorithm enters the next iteration.Figure 3
Example of WTA solutions (2 ants, 4 missiles, and 6 targets).


## 3.1. Pareto Ant Colony Optimization Algorithm (P-ACO)

The ACO is a heuristic algorithm which has emerged in recent years and has been applied to solve hard combinatorial optimization problems [35]. Scientists have developed ACO by imitating the behaviors of real ant colonies in the process of searching for food. When ants in the colony find a short path for gathering food by frequently travelling between their nest and the food source, they deposit a chemical calledpheromone trail on the path which can be followed by other ants to search for the food source [36].Pareto ant colony optimization (P-ACO) algorithm is one of the multiple-objective ant colony optimization (MOACO) algorithms that was proposed by Doerner et al. [28] to be originally applied to optimize the multiobjective portfolio selection. The main objective of P-ACO algorithm is to find a set of solutions ordered by fronts under the concept of Pareto front. Although it is based on the classical ACS, the update process ofglobal pheromone trail is implemented by using two different ants (the best and the second-best solutions) which are generated by each objective g in current iteration. In P-ACO, τg denotes the pheromone matrices of objective g. During each iteration, wg represents the weight of objective g and is used to combine the pheromone information and heuristic information. In the route, an ant selects the next node s to be visited from node r, and it uses the ACS transition rule like formula (5). ζ is a random number from 0 to 1. ζ0 is a fixed parameter (ζ0∈0,1). When ζ≤ζ0, the available knowledge is exploited for selecting the best option through pheromone trail heuristic information. The formula is shown below: (5)j=argmaxj∈Ω∑g=1Gwg·τijgα·ηijβ,ifζ≤ζ0j^,otherwise,where g=1,2,…,G is a set of objectives, ηij represents the heuristic information of edge eij, and j^ is a probabilistic node. α is the weight coefficient of pheromone information and β is the weight coefficient of heuristic information. α and β belong to R. Ω is the current feasible neighborhood of the ant.Whenζ>ζ0, pij is implemented for finding j^ and the formula is expressed as(6)pij=∑g=1Gwg·τijgα·ηijβ∑u∈Ω∑g=1Gwg·τiugα·ηiuβ,ifj∈Ω0,otherwise.When an ant has travelled an edgeeij in the algorithm,local pheromone update (online pheromone trail update) is implemented for each objective g by the rule(7)τijg=1-ρ·τijg+ρ·τ0,where ρ is the pheromone evaporation rate (ρ∈0,1) and τ0 is the value of initial pheromone trail. In this paper, we use the same method as the literature [37] to set the value of τ0.When each ant has constructed its solution,global pheromone update (offline pheromone trail update) is implemented for each objective g. The update formula is shown below:(8)τijg=1-ρ·τijg+ρ·Δτijg,where Δτijg follows the problem-specific values:(9)Δτijg=10ifedgeeij∈bestsolution,5ifedgeeij∈second-bestsolution,0otherwise.In the process, the nondominated solutions are stored in an external setS.The pseudocode of the P-ACO algorithm can be demonstrated as shown in Algorithm1.Algorithm 1
(0) Pre-processing;
(1) Initialization: generate the number of ants A, initialize pheromone vectors with τ0,
set the values ofα and β;
(2) repeat until Termination conditions has not been met
(3)  for Ant = 1 to A
determine the lifespanΘ of the ant randomly of the interval 1,2,…,N;
createx=0,…,0/∗built empty assignment solution∗/;
determine the objective weightwg for each objective g;
indicates the number of targets to be selectedξ=Θ;
(4)  while  ξ>0  and  ∃i with ηix>0
select an assignmentij using formula (5) and add it to x;
update local pheromone information;
decrementξ;
check feasibilityx;
(5)    if  x is feasible
check efficiencyx;
(6)    if  x is efficient
storex and remove dominated ones;
(7) for each objective g
determine best and second-best solutionΔτijg for each objective g;
update global pheromone information using best and second-best solution;
using formula (8);
.


## 3.2. Application of ACO Algorithm for Solving WTA Problem

At the beginning, each missile has a list of possible target assignments (Figure3). The number of assignments is independent of the conditions which are set at first. Each ant builds a tour (WTA solution) while moving from one to another list until all of them have been visited. Figure 3 gives an example with two ants, four missiles, and six targets during an iteration. M=1,2,…,m is a set of missiles. The size of each list is fixed to the number of targets (T={1,2,…,n}). To construct a WTA solution (a tour), the 1st ant visits Target 2 (missile 1), Target 5 (missile 2), Target 4 (missile 3), and Target 3 (missile 4) by using a transition rule which takes the pheromone trails of visited nodes and the heuristic information into account. On the contrary, the 2nd ant visits Target 4 (missile 1), Target 3 (missile 2), Target 6 (missile 3), and Target 1 (missile 4). If the 1st solution is better than the 2nd solution, then the 1st solution is the best solution at current iteration and the algorithm enters the next iteration.Figure 3
Example of WTA solutions (2 ants, 4 missiles, and 6 targets).


## 4. Improvements of the MPACO Algorithm

### 4.1. Dynamic Heuristic Information Calculation Approach

The heuristic informationηij denotes the probability of the ant moving from i to j and is an optional weighting function calculated by heuristic approach. Because ηij must be calculated for all ants in each iteration, it may seriously affect the efficiency of the ACO algorithm. How to design the problem-specific heuristic information function is one of the important problems that should be solved. There are two main calculation approaches of heuristic information at present: (1) ηij is calculated at first and remains unchanged in each iteration until the algorithm finishes [23] and (2) compared with the first approach ηij is given a dynamic value according to the current iteration [38]. In the first approach, ηij is more efficient but cannot completely indicate the probability. The second approach can satisfy precise estimation of the probability but will relatively reduce the efficiency of computation. In order to solve the WTA problem, an improved calculating approach for heuristic information which considers the efficiency of calculation and the quality of solution is proposed. The values of biobjective functions are considered in the approach. The formulas can be shown as follows:(10)ηij1=∑j=1n1-∏i=1m1-ρkqijxij,ηij2=∑i=1m∑j=1ncixij,ηij=ϑ·ηij1ηij2,where ηij1 and ηij2 dividedly denote the first objective function and the second one. ηij is the total probability heuristic information which is represented as ηij1 multiplied by ηij2. ϑ is a positive constant.


### 4.2. Improved Movement Probability Rule

In practice, ACO algorithm is easy to plunge into the local optimum which is caused by stagnation phenomenon. So, the dynamic movement probability rule based on essential selection (q≤q0) and random selection (q>q0) is proposed to improve the strategy of global search in the process. The phenomenon τij of the paths which have been visited is enhanced. The selected probability of the paths will increase in the next iteration. Because few ants visit the path, some better solutions cannot be stored at the beginning iterations. If more ants can visit the path less traveled in the process, it will be beneficial to the solutions space and effectively perfect the P-ACO algorithm. The probability pij can be computed as(11)pij=∑g=1Gpg·τijgα·ηijβ·zij∑u∈Ω∑g=1Gpg·τiukα·ηiuβ·ziu,ifj∈Ω0,otherwise,zij=A·IcA·Ic+δ·Asf·ηij/ηmax,where A denotes the number of ants according to the specific problem, Ic represents the current iteration, δ is a constant, Asf is the total number of ants that have visited the same node so far, and ηmax is the maximum value of ηij. The value of zij takes into account Asf and η simultaneously. If the current iteration tends to the sub-Pareto solution, the pheromone τij of sub-Pareto solution will be enhanced. But the number of Asf increased and the value of zij decreased. This strategy can restrain the excessive growth of pheromone and escape from the local optimum.


### 4.3. Dynamic Evaporation Rate Strategy

In classical P-ACO algorithm, the pheromone evaporation rateρ is a fixed parameter (ρ∈0,1). The value of ρ is one of the main factors which have an influence on global search ability and convergence speed. When P-ACO algorithm is applied to solve large-scale problem, the value of pheromone τ is close to 0 and global search ability declines. If the value of pheromone becomes very large, the selection probability of the passed nodes increases and global search ability is also directly affected. How to design the value of pheromone is crucial to the release and evaporation of pheromone. In order to improve the global search ability, the common strategy is to set a large value of ρ at the beginning of the algorithm. But the evaporation rate will also decay in the process, so that the P-ACO algorithm can quickly fall into local optimum. Two advantages are brought about by the dynamic evaporation rate strategy: (1) better global search ability and (2) faster convergence speed. The mathematical expression of decay model is shown below:(12)ρIc=Imax·Ic·∑g=1Gwgτmaxg-τmingImax-1+∑g=1GwgImax·τming-τmaxgImax-1,where Imax is the maximum iteration and is set at the beginning and Ic is the current iteration. τmaxg denotes the maximum of pheromone and τming denotes the minimum of pheromone.


### 4.4. Global Updating Rule of Pheromone

A fixed pheromone update rule is mainly used in the classical P-ACO algorithm. The traditional strategy may flung the algorithm into local optimum and has three disadvantages: (1) the distribution features of solution are ignored; (2) there is the trend of stagnation; (3) the convergence rate is slow. In addition to the pheromone local updating rule, the proposed MPACO adopts an adaptive adjustment global updating rule which can effectively deal with the conflict between expanding search and finding Pareto solution. In this rule, the pheromone information of each objective is updated by the best and second-best solution. The rule also provides some instructive guidance to improve the quality of the solutions. The pheromoneΔτijg can be computed as(13)Δτijg=1fbgifxbest,ijg=11fsbgifxsecond-best,ijg=10otherwise,where xij is an assignment of a population best (or second-best) ant according to objective g. The strategy maintains the Pareto exploitation and balances the random search and the function about pheromone evaporation or augmentation. In order to reduce the otherness of pheromone on the local Pareto solution and expand the search scope, the positive feedback of P-ACO algorithm is restrained by adding a little amount of negative feedback pheromone in the process. If the algorithm falls into a local Pareto point during a period of iterations, then the adaptive adjustment strategy will be adopted to decrease the amount of pheromone for escaping from local Pareto front. The two expressions are similar in format, but the objective function is different.


### 4.5. Boundary Symmetric Mutation Strategy

Based on the mathematical statistics theory, the distribution of the samples will be close to the normal distribution when the number of samples tends to infinity. In solving WTA, targets coordinates according to the values are sorted. When a larger scale problem is solved, the distributions of targets are complex and a majority of targets are on the central region. The path of ants is boundary-center-boundary in the initial stage, so that all ants will follow the original track mode.In order to solve the above problem, the mutation strategy is used for boundary mutation which can improve the mutation efficiency and the quality of solutions. Based on the quartile idea in probability theory, the mutation strategy is applied to mutate 1/4 solutions in both sides of the path. Due to the symmetry of the boundary, crossover and mutation only take place in the selected boundary which did not happen among others.According to the above improvements, a modified Pareto ant colony optimization (MPACO) algorithm based on dynamic calculation approach of heuristic information, improved movement probability rules, improvement strategy of evaporation factor, improved global updating rules of pheromone, and boundary symmetric mutation strategy is proposed in this paper.


### 4.6. The Structure of MPACO Algorithm

According to the idea of the proposed MPACO algorithm, the structure of MPACO algorithm in solving BOWTA problem can be shown as follows.BeginStep
0 (preprocessing). These nodes of decision table are divided into center nodes and boundary nodes according to qij which represents the effectiveness of the missile. The classification result is saved and sorted in order to prepare for the subsequent run of the MPACO algorithm.Step
1 (initialization of parameters for the proposed MPACO). These parameters include the colony of ants A, the maximum iteration Imax, the pheromone information factor α, heuristic information factor β, and initial uniform probability ζ0.While (Termination conditions has not been met)doForI=1,2,…,Imax.ForA=1,2,…,Amax.Step
2. The A ants randomly visit targets of each missile by constraints (2), (3), and (4), and these assignments are added into the decision list of the ant Xa=xij1,xij2,…,xijm.Step
3. Due to the dynamic calculation approach of heuristic information, formulas (10) are used to calculate coupling ηij before Step  4.Step
4. According to (11) of improved movement probability rules, the probability that the ant selects the next node is calculated (Section 4.2). Then, this node is added into the decision list, and the pheromone is locally updated.Step
5. After the ants have completed a selection, the fitness function f will be calculated. Then, the decision list is modified. Repeat Step  4 until the ant completed the tour (namely, all missiles are allocated). The current best and second-best solutions are saved and the global best solution is updated for the current iteration.End forStep
6. Based on (8) and (13), the pheromone is globally updated.Step
7. Update Pareto optimal front set.Step
8. Set the iterative counter I=I+1 for the iteration control. If I<Imax, return to Step  2.End whileStep
9. The proposed MPACO algorithm is terminated; output Pareto solutions.


## 4.1. Dynamic Heuristic Information Calculation Approach

The heuristic informationηij denotes the probability of the ant moving from i to j and is an optional weighting function calculated by heuristic approach. Because ηij must be calculated for all ants in each iteration, it may seriously affect the efficiency of the ACO algorithm. How to design the problem-specific heuristic information function is one of the important problems that should be solved. There are two main calculation approaches of heuristic information at present: (1) ηij is calculated at first and remains unchanged in each iteration until the algorithm finishes [23] and (2) compared with the first approach ηij is given a dynamic value according to the current iteration [38]. In the first approach, ηij is more efficient but cannot completely indicate the probability. The second approach can satisfy precise estimation of the probability but will relatively reduce the efficiency of computation. In order to solve the WTA problem, an improved calculating approach for heuristic information which considers the efficiency of calculation and the quality of solution is proposed. The values of biobjective functions are considered in the approach. The formulas can be shown as follows:(10)ηij1=∑j=1n1-∏i=1m1-ρkqijxij,ηij2=∑i=1m∑j=1ncixij,ηij=ϑ·ηij1ηij2,where ηij1 and ηij2 dividedly denote the first objective function and the second one. ηij is the total probability heuristic information which is represented as ηij1 multiplied by ηij2. ϑ is a positive constant.


## 4.2. Improved Movement Probability Rule

In practice, ACO algorithm is easy to plunge into the local optimum which is caused by stagnation phenomenon. So, the dynamic movement probability rule based on essential selection (q≤q0) and random selection (q>q0) is proposed to improve the strategy of global search in the process. The phenomenon τij of the paths which have been visited is enhanced. The selected probability of the paths will increase in the next iteration. Because few ants visit the path, some better solutions cannot be stored at the beginning iterations. If more ants can visit the path less traveled in the process, it will be beneficial to the solutions space and effectively perfect the P-ACO algorithm. The probability pij can be computed as(11)pij=∑g=1Gpg·τijgα·ηijβ·zij∑u∈Ω∑g=1Gpg·τiukα·ηiuβ·ziu,ifj∈Ω0,otherwise,zij=A·IcA·Ic+δ·Asf·ηij/ηmax,where A denotes the number of ants according to the specific problem, Ic represents the current iteration, δ is a constant, Asf is the total number of ants that have visited the same node so far, and ηmax is the maximum value of ηij. The value of zij takes into account Asf and η simultaneously. If the current iteration tends to the sub-Pareto solution, the pheromone τij of sub-Pareto solution will be enhanced. But the number of Asf increased and the value of zij decreased. This strategy can restrain the excessive growth of pheromone and escape from the local optimum.


## 4.3. Dynamic Evaporation Rate Strategy

In classical P-ACO algorithm, the pheromone evaporation rateρ is a fixed parameter (ρ∈0,1). The value of ρ is one of the main factors which have an influence on global search ability and convergence speed. When P-ACO algorithm is applied to solve large-scale problem, the value of pheromone τ is close to 0 and global search ability declines. If the value of pheromone becomes very large, the selection probability of the passed nodes increases and global search ability is also directly affected. How to design the value of pheromone is crucial to the release and evaporation of pheromone. In order to improve the global search ability, the common strategy is to set a large value of ρ at the beginning of the algorithm. But the evaporation rate will also decay in the process, so that the P-ACO algorithm can quickly fall into local optimum. Two advantages are brought about by the dynamic evaporation rate strategy: (1) better global search ability and (2) faster convergence speed. The mathematical expression of decay model is shown below:(12)ρIc=Imax·Ic·∑g=1Gwgτmaxg-τmingImax-1+∑g=1GwgImax·τming-τmaxgImax-1,where Imax is the maximum iteration and is set at the beginning and Ic is the current iteration. τmaxg denotes the maximum of pheromone and τming denotes the minimum of pheromone.


## 4.4. Global Updating Rule of Pheromone

A fixed pheromone update rule is mainly used in the classical P-ACO algorithm. The traditional strategy may flung the algorithm into local optimum and has three disadvantages: (1) the distribution features of solution are ignored; (2) there is the trend of stagnation; (3) the convergence rate is slow. In addition to the pheromone local updating rule, the proposed MPACO adopts an adaptive adjustment global updating rule which can effectively deal with the conflict between expanding search and finding Pareto solution. In this rule, the pheromone information of each objective is updated by the best and second-best solution. The rule also provides some instructive guidance to improve the quality of the solutions. The pheromoneΔτijg can be computed as(13)Δτijg=1fbgifxbest,ijg=11fsbgifxsecond-best,ijg=10otherwise,where xij is an assignment of a population best (or second-best) ant according to objective g. The strategy maintains the Pareto exploitation and balances the random search and the function about pheromone evaporation or augmentation. In order to reduce the otherness of pheromone on the local Pareto solution and expand the search scope, the positive feedback of P-ACO algorithm is restrained by adding a little amount of negative feedback pheromone in the process. If the algorithm falls into a local Pareto point during a period of iterations, then the adaptive adjustment strategy will be adopted to decrease the amount of pheromone for escaping from local Pareto front. The two expressions are similar in format, but the objective function is different.


## 4.5. Boundary Symmetric Mutation Strategy

Based on the mathematical statistics theory, the distribution of the samples will be close to the normal distribution when the number of samples tends to infinity. In solving WTA, targets coordinates according to the values are sorted. When a larger scale problem is solved, the distributions of targets are complex and a majority of targets are on the central region. The path of ants is boundary-center-boundary in the initial stage, so that all ants will follow the original track mode.In order to solve the above problem, the mutation strategy is used for boundary mutation which can improve the mutation efficiency and the quality of solutions. Based on the quartile idea in probability theory, the mutation strategy is applied to mutate 1/4 solutions in both sides of the path. Due to the symmetry of the boundary, crossover and mutation only take place in the selected boundary which did not happen among others.According to the above improvements, a modified Pareto ant colony optimization (MPACO) algorithm based on dynamic calculation approach of heuristic information, improved movement probability rules, improvement strategy of evaporation factor, improved global updating rules of pheromone, and boundary symmetric mutation strategy is proposed in this paper.


## 4.6. The Structure of MPACO Algorithm

According to the idea of the proposed MPACO algorithm, the structure of MPACO algorithm in solving BOWTA problem can be shown as follows.BeginStep
0 (preprocessing). These nodes of decision table are divided into center nodes and boundary nodes according to qij which represents the effectiveness of the missile. The classification result is saved and sorted in order to prepare for the subsequent run of the MPACO algorithm.Step
1 (initialization of parameters for the proposed MPACO). These parameters include the colony of ants A, the maximum iteration Imax, the pheromone information factor α, heuristic information factor β, and initial uniform probability ζ0.While (Termination conditions has not been met)doForI=1,2,…,Imax.ForA=1,2,…,Amax.Step
2. The A ants randomly visit targets of each missile by constraints (2), (3), and (4), and these assignments are added into the decision list of the ant Xa=xij1,xij2,…,xijm.Step
3. Due to the dynamic calculation approach of heuristic information, formulas (10) are used to calculate coupling ηij before Step  4.Step
4. According to (11) of improved movement probability rules, the probability that the ant selects the next node is calculated (Section 4.2). Then, this node is added into the decision list, and the pheromone is locally updated.Step
5. After the ants have completed a selection, the fitness function f will be calculated. Then, the decision list is modified. Repeat Step  4 until the ant completed the tour (namely, all missiles are allocated). The current best and second-best solutions are saved and the global best solution is updated for the current iteration.End forStep
6. Based on (8) and (13), the pheromone is globally updated.Step
7. Update Pareto optimal front set.Step
8. Set the iterative counter I=I+1 for the iteration control. If I<Imax, return to Step  2.End whileStep
9. The proposed MPACO algorithm is terminated; output Pareto solutions.


## 5. Experiment Results and Analysis

### 5.1. Parameter Settings

In order to select the initial value of pheromoneτ0, we use the formula in the literature [39], which can be shown as follows:(14)τ0=G∑g=1GfgXg,where Xg is a solution of the gth objective computed by a greedy algorithm, G is the number of objectives in problem-specific function, and fgXg is the fitness function of the solution Xg according to the gth objective.Based on the numbers of targets and missiles in a different BOWTA scale, the maximum number of iterationsImax is set. Some parameters of MPACO algorithm are shown in Table 3.Table 3
Parameters of MPACO algorithm (part).
Parameter
Value
A
20
α
1
β
2
w
1,2
0.5
τ
0
Formula (14)
ζ
0
0.98We can see from Table4 that the numerical experiments have been performed on a PC with Intel(R) Core(TM) i5-4460T CPU @ 1.90 GHz, with 8.00 GB of RAM, running Windows 7 64 bits as the operating system. The codes of all algorithms and optimization models have been written in C++ under Microsoft Visual C++ 6.0.Table 4
Parameters of the simulation environment.
Computer specifications
Intel(R) Core(TM) i5-4460T CPU @ 1.90 GHz with 8.00 GB RAM
Operating system
Windows 7 (x64) Operation System
Language
C++
Software
Microsoft Visual C++ 6.0


### 5.2. Numerical Experiments and Analysis

#### 5.2.1. Computation Results on Specific Example

A specific example is designed to verify performances of the algorithm. Assume that there are 4 fighters that carry different numbers of missiles (12 missile units in total) and 10 targets in the simulation. Missiles damage probabilitypij  i=1,…,12;j=1,…,10, pilot operation factor ρk  k=1,…,4, and cost of missiles ci are shown in Tables 5–7.Table 5
The table of missiles damage probability.
Missile unit
Target
1
2
3
4
5
6
7
8
9
10
1
0.82
0.85
0.78
0.75
0.52
0.88
0.44
0.76
0.72
0.56
2
0.56
0.72
0.88
0.46
0.64
0.47
0.68
0.45
0.48
0.75
3
0.45
0.62
0.54
0.73
0.84
0.76
0.78
0.42
0.53
0.65
4
0.56
0.42
0.76
0.84
0.73
0.83
0.86
0.62
0.78
0.82
5
0.45
0.58
0.81
0.44
0.63
0.59
0.78
0.77
0.65
0.70
6
0.46
0.61
0.55
0.68
0.75
0.83
0.73
0.66
0.82
0.48
7
0.66
0.71
0.65
0.44
0.86
0.79
0.59
0.85
0.53
0.56
8
0.56
0.42
0.76
0.84
0.73
0.72
0.44
0.75
0.48
0.47
9
0.88
0.78
0.44
0.67
0.56
0.86
0.58
0.65
0.73
0.42
10
0.56
0.88
0.68
0.45
0.75
0.73
0.61
0.76
0.84
0.78
11
0.84
0.54
0.44
0.42
0.65
0.56
0.71
0.55
0.45
0.88
12
0.83
0.76
0.84
0.62
0.82
0.75
0.42
0.68
0.57
0.54Table 6
Pilot operation factor table.
Fighter
Missile unit
ρ
k
F
1
M
1
0.96
M
2
M
3
F
2
M
4
0.95
M
5
F
3
M
6
0.98
M
7
M
8
M
9
F
4
M
10
0.93
M
11
M
12Table 7
The cost of each missile table.
Missile unit
M
1
M
2
M
3
M
4
M
5
M
6
M
7
M
8
M
9
M
10
M
11
M
12
c
i
0.62
0.63
0.69
0.80
0.72
0.90
0.96
0.68
0.72
0.65
0.66
0.65According to the literature [1], we present the missiles cost, as shown in Table 7.Firstly, an enumeration method is employed to the specific example instance and a set of evenly distributed true Pareto solutions (PSs) are obtained. Secondly, the proposed algorithm is used to find PSs.Figure4 presents the distribution of the true PSs solved by enumeration method and final PSs obtained by the MPACO in the objective space. We can see from the figure that the MPACO can find close and even PSs in the objective space. It is evident that MPACO can guarantee the quality of solutions. So, the MPACO for optimizing the BOWTA is well verified to be feasible.Figure 4
Plots of the true PS and final PSs found by MPACO on the specific example instance.We can get some results from Figure5 and Table 8 as follows: when funds of national defense are enough, we can choose Scheme 1 which cost the most money and obtain the greatest expected damage of the enemy. When funds are insufficient, we should choose one scheme among Schemes 34, 35, and 36 on the basis of threat coefficient of each target. Considering that all targets must be allocated and the cost of missiles should be minimized, we can only choose Scheme 4.Table 8
The value of PSs and results of WTA distribution.
1/f1
f2
M
1
M
2
M
3
M
4
M
5
M
6
M
7
M
8
M
9
M
10
M
11
M
12
1
0.056798
8.68
6
3
7
7
8
9
5
4
1
9
10
2
2
0.058623
7.99
6
3
0
7
8
9
5
4
1
9
10
2
3
0.059
7.72
6
3
5
7
8
9
0
4
1
9
10
2
4
0.060966
7.34
6
3
0
7
8
9
5
4
1
2
10
0
5
0.061662
7.00
6
3
5
7
0
9
0
4
1
9
10
2
6
0.063811
6.62
6
3
0
7
0
9
5
4
1
2
10
0
7
0.064259
6.35
6
3
5
7
0
9
0
4
1
2
10
0
8
0.067327
5.96
6
3
0
7
0
9
5
4
1
2
0
0
9
0.067848
5.67
6
3
5
7
0
9
0
0
1
2
10
0
10
0.068893
5.42
6
3
5
7
0
0
0
0
1
9
10
2
11
0.071974
5.05
0
3
5
7
0
9
0
0
1
2
10
0
12
0.072395
4.79
6
3
5
7
0
0
0
4
1
9
0
0
13
0.073822
4.73
6
3
0
7
8
0
0
0
0
9
10
5
14
0.07645
4.42
0
3
5
7
0
9
0
4
1
0
0
0
15
0.077215
4.07
6
3
5
7
0
0
0
4
0
9
0
0
16
0.081585
3.74
0
3
5
7
0
9
0
0
1
0
0
0
17
0.081702
3.64
6
3
5
7
0
9
0
0
0
0
0
0
18
0.082457
3.39
6
3
5
7
0
0
0
0
0
9
0
0
19
0.085299
3.24
6
3
7
0
0
0
0
0
0
9
0
5
20
0.087759
3.02
0
3
5
7
0
9
0
0
0
0
0
0
21
0.08863
2.77
0
3
5
7
0
0
0
0
0
9
0
0
22
0.089675
2.73
0
3
0
7
0
0
0
0
0
9
0
5
23
0.09589
2.59
9
3
7
0
0
0
0
0
0
0
0
5
24
0.103038
2.39
0
0
5
7
0
9
0
0
0
0
0
0
25
0.104241
2.14
0
0
5
7
0
0
0
0
0
9
0
0
26
0.108826
1.99
0
0
7
0
0
0
0
0
0
9
0
5
27
0.116404
1.97
0
3
7
0
0
0
0
0
0
9
0
0
28
0.119689
1.90
9
7
0
0
0
0
0
0
0
0
0
5
29
0.13725
1.70
0
0
0
7
0
9
0
0
0
0
0
0
30
0.139392
1.45
0
0
0
7
0
0
0
0
0
9
0
0
31
0.144902
1.34
0
0
7
0
0
0
0
0
0
9
0
0
32
0.15344
1.28
0
7
0
0
0
0
0
0
0
9
0
0
33
0.164821
1.25
9
7
0
0
0
0
0
0
0
0
0
0
34
0.24888
0.90
0
0
0
0
0
9
0
0
0
0
0
0
35
0.256016
0.65
0
0
0
0
0
0
0
0
0
9
0
0
36
0.289352
0.62
9
0
0
0
0
0
0
0
0
0
0
0Figure 5
Plots of PSs best performance in 30 runs of four algorithms on Scale1.


#### 5.2.2. Comparison Study

In order to verify the applicability and feasibility of the proposed algorithm, we test MPACO, P-ACO, NSGA-II, and SPEA-II on Scale1, Scale 2, and Scale 3 instances. Three examples represent, respectively, small-, medium-, and large-scale BOWTA instances under the same condition. The quantitative indexes, such as Pareto figure and time performance, are chosen to evaluate the performance on C++. The value of three instances can be seen from Table 9.Table 9
The values of Scale1, Scale 2, and Scale 3 instances.
Scale1
Scale2
Scale3
Target
10
48
98
Missile
10
53
100In what follows, Scale1 instance is firstly taken as an example to illustrate the performance of MPACO, P-ACO, NSGA-II, and SPEA-II. Then, Scale 2 and Scale 3 instances are introduced to further demonstrate the performances of both algorithms. The statistical results are shown in Figures 5–7.Figure 6
Plots of PSs best performance in 30 runs of four algorithms on Scale2.Figure 7
Plots of PSs best performance in 30 runs of four algorithms on Scale3.Figures5–7 present the distribution of PSs obtained in 30 runs with the best performance of each algorithm for each test instance in objective space. When especially tested on Scale 2 and Scale 3, the MPACO can obtain the best solutions among the four algorithms. Due to increments of dynamic calculation approach of heuristic information and improved movement probability rules, the calculation accuracy of MPACO algorithm is improved. Boundary symmetric mutation strategy is also good for the algorithm to solve large-scale BOWTA problem.The primary aim of solving BOWTA problems is to lay the foundation for the DWTA problems whose goal is to provide a set of near-optimal or acceptable real-time decisions in air combat. So, in practice, the time performance of algorithms is also an important index. In the end, we test four algorithms on three instances in 30 runs and record the literation time of each algorithm. The statistical results of time performances are shown in Figures8–10.Figure 8
Boxplot of time performance on Scale1 (in secs).Figure 9
Boxplot of time performance on Scale2 (in secs).Figure 10
Boxplot of time performance on Scale3 (in secs).Although MPACO does not have time advantage in solving Scale1 and Scale 2 compared with other algorithms, it runs as fast as SPEA-II with the same number of function evaluations on Scale 3. Boundary symmetric mutation strategy improves the mutation efficiency and the quality of solutions and new evaporation factor is good for solving large-scale problems, but dynamic heuristic information and probability rules affect the iteration rate. In real air combat situations, pilots often make deadly decisions within seconds or even within milliseconds. Although the algorithm structure and implementation process are relatively complex and increase the computational complexity, high-performance hardware platform can be adopted to reduce the processing time and satisfy the real-time requirement.In this section, we do some work as follows. Firstly, we test the MPACO on a specific example. We verify the feasibility after comparing two kinds of PSs and show the corresponding distribution results. Secondly, to give a comparison study among the four algorithms, three different scale instances are introduced to evaluate each algorithm. Plots of the best PSs in 30 runs of each algorithm are presented to have an intuitive comparison. Lastly, we show the time performance of the four algorithms on each instance. To summarize, MPACO has been proven to be an effective technique for BOWTA optimization problem and obviously is the best among the four algorithms.


## 5.1. Parameter Settings

In order to select the initial value of pheromoneτ0, we use the formula in the literature [39], which can be shown as follows:(14)τ0=G∑g=1GfgXg,where Xg is a solution of the gth objective computed by a greedy algorithm, G is the number of objectives in problem-specific function, and fgXg is the fitness function of the solution Xg according to the gth objective.Based on the numbers of targets and missiles in a different BOWTA scale, the maximum number of iterationsImax is set. Some parameters of MPACO algorithm are shown in Table 3.Table 3
Parameters of MPACO algorithm (part).
Parameter
Value
A
20
α
1
β
2
w
1,2
0.5
τ
0
Formula (14)
ζ
0
0.98We can see from Table4 that the numerical experiments have been performed on a PC with Intel(R) Core(TM) i5-4460T CPU @ 1.90 GHz, with 8.00 GB of RAM, running Windows 7 64 bits as the operating system. The codes of all algorithms and optimization models have been written in C++ under Microsoft Visual C++ 6.0.Table 4
Parameters of the simulation environment.
Computer specifications
Intel(R) Core(TM) i5-4460T CPU @ 1.90 GHz with 8.00 GB RAM
Operating system
Windows 7 (x64) Operation System
Language
C++
Software
Microsoft Visual C++ 6.0


## 5.2. Numerical Experiments and Analysis

### 5.2.1. Computation Results on Specific Example

A specific example is designed to verify performances of the algorithm. Assume that there are 4 fighters that carry different numbers of missiles (12 missile units in total) and 10 targets in the simulation. Missiles damage probabilitypij  i=1,…,12;j=1,…,10, pilot operation factor ρk  k=1,…,4, and cost of missiles ci are shown in Tables 5–7.Table 5
The table of missiles damage probability.
Missile unit
Target
1
2
3
4
5
6
7
8
9
10
1
0.82
0.85
0.78
0.75
0.52
0.88
0.44
0.76
0.72
0.56
2
0.56
0.72
0.88
0.46
0.64
0.47
0.68
0.45
0.48
0.75
3
0.45
0.62
0.54
0.73
0.84
0.76
0.78
0.42
0.53
0.65
4
0.56
0.42
0.76
0.84
0.73
0.83
0.86
0.62
0.78
0.82
5
0.45
0.58
0.81
0.44
0.63
0.59
0.78
0.77
0.65
0.70
6
0.46
0.61
0.55
0.68
0.75
0.83
0.73
0.66
0.82
0.48
7
0.66
0.71
0.65
0.44
0.86
0.79
0.59
0.85
0.53
0.56
8
0.56
0.42
0.76
0.84
0.73
0.72
0.44
0.75
0.48
0.47
9
0.88
0.78
0.44
0.67
0.56
0.86
0.58
0.65
0.73
0.42
10
0.56
0.88
0.68
0.45
0.75
0.73
0.61
0.76
0.84
0.78
11
0.84
0.54
0.44
0.42
0.65
0.56
0.71
0.55
0.45
0.88
12
0.83
0.76
0.84
0.62
0.82
0.75
0.42
0.68
0.57
0.54Table 6
Pilot operation factor table.
Fighter
Missile unit
ρ
k
F
1
M
1
0.96
M
2
M
3
F
2
M
4
0.95
M
5
F
3
M
6
0.98
M
7
M
8
M
9
F
4
M
10
0.93
M
11
M
12Table 7
The cost of each missile table.
Missile unit
M
1
M
2
M
3
M
4
M
5
M
6
M
7
M
8
M
9
M
10
M
11
M
12
c
i
0.62
0.63
0.69
0.80
0.72
0.90
0.96
0.68
0.72
0.65
0.66
0.65According to the literature [1], we present the missiles cost, as shown in Table 7.Firstly, an enumeration method is employed to the specific example instance and a set of evenly distributed true Pareto solutions (PSs) are obtained. Secondly, the proposed algorithm is used to find PSs.Figure4 presents the distribution of the true PSs solved by enumeration method and final PSs obtained by the MPACO in the objective space. We can see from the figure that the MPACO can find close and even PSs in the objective space. It is evident that MPACO can guarantee the quality of solutions. So, the MPACO for optimizing the BOWTA is well verified to be feasible.Figure 4
Plots of the true PS and final PSs found by MPACO on the specific example instance.We can get some results from Figure5 and Table 8 as follows: when funds of national defense are enough, we can choose Scheme 1 which cost the most money and obtain the greatest expected damage of the enemy. When funds are insufficient, we should choose one scheme among Schemes 34, 35, and 36 on the basis of threat coefficient of each target. Considering that all targets must be allocated and the cost of missiles should be minimized, we can only choose Scheme 4.Table 8
The value of PSs and results of WTA distribution.
1/f1
f2
M
1
M
2
M
3
M
4
M
5
M
6
M
7
M
8
M
9
M
10
M
11
M
12
1
0.056798
8.68
6
3
7
7
8
9
5
4
1
9
10
2
2
0.058623
7.99
6
3
0
7
8
9
5
4
1
9
10
2
3
0.059
7.72
6
3
5
7
8
9
0
4
1
9
10
2
4
0.060966
7.34
6
3
0
7
8
9
5
4
1
2
10
0
5
0.061662
7.00
6
3
5
7
0
9
0
4
1
9
10
2
6
0.063811
6.62
6
3
0
7
0
9
5
4
1
2
10
0
7
0.064259
6.35
6
3
5
7
0
9
0
4
1
2
10
0
8
0.067327
5.96
6
3
0
7
0
9
5
4
1
2
0
0
9
0.067848
5.67
6
3
5
7
0
9
0
0
1
2
10
0
10
0.068893
5.42
6
3
5
7
0
0
0
0
1
9
10
2
11
0.071974
5.05
0
3
5
7
0
9
0
0
1
2
10
0
12
0.072395
4.79
6
3
5
7
0
0
0
4
1
9
0
0
13
0.073822
4.73
6
3
0
7
8
0
0
0
0
9
10
5
14
0.07645
4.42
0
3
5
7
0
9
0
4
1
0
0
0
15
0.077215
4.07
6
3
5
7
0
0
0
4
0
9
0
0
16
0.081585
3.74
0
3
5
7
0
9
0
0
1
0
0
0
17
0.081702
3.64
6
3
5
7
0
9
0
0
0
0
0
0
18
0.082457
3.39
6
3
5
7
0
0
0
0
0
9
0
0
19
0.085299
3.24
6
3
7
0
0
0
0
0
0
9
0
5
20
0.087759
3.02
0
3
5
7
0
9
0
0
0
0
0
0
21
0.08863
2.77
0
3
5
7
0
0
0
0
0
9
0
0
22
0.089675
2.73
0
3
0
7
0
0
0
0
0
9
0
5
23
0.09589
2.59
9
3
7
0
0
0
0
0
0
0
0
5
24
0.103038
2.39
0
0
5
7
0
9
0
0
0
0
0
0
25
0.104241
2.14
0
0
5
7
0
0
0
0
0
9
0
0
26
0.108826
1.99
0
0
7
0
0
0
0
0
0
9
0
5
27
0.116404
1.97
0
3
7
0
0
0
0
0
0
9
0
0
28
0.119689
1.90
9
7
0
0
0
0
0
0
0
0
0
5
29
0.13725
1.70
0
0
0
7
0
9
0
0
0
0
0
0
30
0.139392
1.45
0
0
0
7
0
0
0
0
0
9
0
0
31
0.144902
1.34
0
0
7
0
0
0
0
0
0
9
0
0
32
0.15344
1.28
0
7
0
0
0
0
0
0
0
9
0
0
33
0.164821
1.25
9
7
0
0
0
0
0
0
0
0
0
0
34
0.24888
0.90
0
0
0
0
0
9
0
0
0
0
0
0
35
0.256016
0.65
0
0
0
0
0
0
0
0
0
9
0
0
36
0.289352
0.62
9
0
0
0
0
0
0
0
0
0
0
0Figure 5
Plots of PSs best performance in 30 runs of four algorithms on Scale1.


### 5.2.2. Comparison Study

In order to verify the applicability and feasibility of the proposed algorithm, we test MPACO, P-ACO, NSGA-II, and SPEA-II on Scale1, Scale 2, and Scale 3 instances. Three examples represent, respectively, small-, medium-, and large-scale BOWTA instances under the same condition. The quantitative indexes, such as Pareto figure and time performance, are chosen to evaluate the performance on C++. The value of three instances can be seen from Table 9.Table 9
The values of Scale1, Scale 2, and Scale 3 instances.
Scale1
Scale2
Scale3
Target
10
48
98
Missile
10
53
100In what follows, Scale1 instance is firstly taken as an example to illustrate the performance of MPACO, P-ACO, NSGA-II, and SPEA-II. Then, Scale 2 and Scale 3 instances are introduced to further demonstrate the performances of both algorithms. The statistical results are shown in Figures 5–7.Figure 6
Plots of PSs best performance in 30 runs of four algorithms on Scale2.Figure 7
Plots of PSs best performance in 30 runs of four algorithms on Scale3.Figures5–7 present the distribution of PSs obtained in 30 runs with the best performance of each algorithm for each test instance in objective space. When especially tested on Scale 2 and Scale 3, the MPACO can obtain the best solutions among the four algorithms. Due to increments of dynamic calculation approach of heuristic information and improved movement probability rules, the calculation accuracy of MPACO algorithm is improved. Boundary symmetric mutation strategy is also good for the algorithm to solve large-scale BOWTA problem.The primary aim of solving BOWTA problems is to lay the foundation for the DWTA problems whose goal is to provide a set of near-optimal or acceptable real-time decisions in air combat. So, in practice, the time performance of algorithms is also an important index. In the end, we test four algorithms on three instances in 30 runs and record the literation time of each algorithm. The statistical results of time performances are shown in Figures8–10.Figure 8
Boxplot of time performance on Scale1 (in secs).Figure 9
Boxplot of time performance on Scale2 (in secs).Figure 10
Boxplot of time performance on Scale3 (in secs).Although MPACO does not have time advantage in solving Scale1 and Scale 2 compared with other algorithms, it runs as fast as SPEA-II with the same number of function evaluations on Scale 3. Boundary symmetric mutation strategy improves the mutation efficiency and the quality of solutions and new evaporation factor is good for solving large-scale problems, but dynamic heuristic information and probability rules affect the iteration rate. In real air combat situations, pilots often make deadly decisions within seconds or even within milliseconds. Although the algorithm structure and implementation process are relatively complex and increase the computational complexity, high-performance hardware platform can be adopted to reduce the processing time and satisfy the real-time requirement.In this section, we do some work as follows. Firstly, we test the MPACO on a specific example. We verify the feasibility after comparing two kinds of PSs and show the corresponding distribution results. Secondly, to give a comparison study among the four algorithms, three different scale instances are introduced to evaluate each algorithm. Plots of the best PSs in 30 runs of each algorithm are presented to have an intuitive comparison. Lastly, we show the time performance of the four algorithms on each instance. To summarize, MPACO has been proven to be an effective technique for BOWTA optimization problem and obviously is the best among the four algorithms.


## 5.2.1. Computation Results on Specific Example

A specific example is designed to verify performances of the algorithm. Assume that there are 4 fighters that carry different numbers of missiles (12 missile units in total) and 10 targets in the simulation. Missiles damage probabilitypij  i=1,…,12;j=1,…,10, pilot operation factor ρk  k=1,…,4, and cost of missiles ci are shown in Tables 5–7.Table 5
The table of missiles damage probability.
Missile unit
Target
1
2
3
4
5
6
7
8
9
10
1
0.82
0.85
0.78
0.75
0.52
0.88
0.44
0.76
0.72
0.56
2
0.56
0.72
0.88
0.46
0.64
0.47
0.68
0.45
0.48
0.75
3
0.45
0.62
0.54
0.73
0.84
0.76
0.78
0.42
0.53
0.65
4
0.56
0.42
0.76
0.84
0.73
0.83
0.86
0.62
0.78
0.82
5
0.45
0.58
0.81
0.44
0.63
0.59
0.78
0.77
0.65
0.70
6
0.46
0.61
0.55
0.68
0.75
0.83
0.73
0.66
0.82
0.48
7
0.66
0.71
0.65
0.44
0.86
0.79
0.59
0.85
0.53
0.56
8
0.56
0.42
0.76
0.84
0.73
0.72
0.44
0.75
0.48
0.47
9
0.88
0.78
0.44
0.67
0.56
0.86
0.58
0.65
0.73
0.42
10
0.56
0.88
0.68
0.45
0.75
0.73
0.61
0.76
0.84
0.78
11
0.84
0.54
0.44
0.42
0.65
0.56
0.71
0.55
0.45
0.88
12
0.83
0.76
0.84
0.62
0.82
0.75
0.42
0.68
0.57
0.54Table 6
Pilot operation factor table.
Fighter
Missile unit
ρ
k
F
1
M
1
0.96
M
2
M
3
F
2
M
4
0.95
M
5
F
3
M
6
0.98
M
7
M
8
M
9
F
4
M
10
0.93
M
11
M
12Table 7
The cost of each missile table.
Missile unit
M
1
M
2
M
3
M
4
M
5
M
6
M
7
M
8
M
9
M
10
M
11
M
12
c
i
0.62
0.63
0.69
0.80
0.72
0.90
0.96
0.68
0.72
0.65
0.66
0.65According to the literature [1], we present the missiles cost, as shown in Table 7.Firstly, an enumeration method is employed to the specific example instance and a set of evenly distributed true Pareto solutions (PSs) are obtained. Secondly, the proposed algorithm is used to find PSs.Figure4 presents the distribution of the true PSs solved by enumeration method and final PSs obtained by the MPACO in the objective space. We can see from the figure that the MPACO can find close and even PSs in the objective space. It is evident that MPACO can guarantee the quality of solutions. So, the MPACO for optimizing the BOWTA is well verified to be feasible.Figure 4
Plots of the true PS and final PSs found by MPACO on the specific example instance.We can get some results from Figure5 and Table 8 as follows: when funds of national defense are enough, we can choose Scheme 1 which cost the most money and obtain the greatest expected damage of the enemy. When funds are insufficient, we should choose one scheme among Schemes 34, 35, and 36 on the basis of threat coefficient of each target. Considering that all targets must be allocated and the cost of missiles should be minimized, we can only choose Scheme 4.Table 8
The value of PSs and results of WTA distribution.
1/f1
f2
M
1
M
2
M
3
M
4
M
5
M
6
M
7
M
8
M
9
M
10
M
11
M
12
1
0.056798
8.68
6
3
7
7
8
9
5
4
1
9
10
2
2
0.058623
7.99
6
3
0
7
8
9
5
4
1
9
10
2
3
0.059
7.72
6
3
5
7
8
9
0
4
1
9
10
2
4
0.060966
7.34
6
3
0
7
8
9
5
4
1
2
10
0
5
0.061662
7.00
6
3
5
7
0
9
0
4
1
9
10
2
6
0.063811
6.62
6
3
0
7
0
9
5
4
1
2
10
0
7
0.064259
6.35
6
3
5
7
0
9
0
4
1
2
10
0
8
0.067327
5.96
6
3
0
7
0
9
5
4
1
2
0
0
9
0.067848
5.67
6
3
5
7
0
9
0
0
1
2
10
0
10
0.068893
5.42
6
3
5
7
0
0
0
0
1
9
10
2
11
0.071974
5.05
0
3
5
7
0
9
0
0
1
2
10
0
12
0.072395
4.79
6
3
5
7
0
0
0
4
1
9
0
0
13
0.073822
4.73
6
3
0
7
8
0
0
0
0
9
10
5
14
0.07645
4.42
0
3
5
7
0
9
0
4
1
0
0
0
15
0.077215
4.07
6
3
5
7
0
0
0
4
0
9
0
0
16
0.081585
3.74
0
3
5
7
0
9
0
0
1
0
0
0
17
0.081702
3.64
6
3
5
7
0
9
0
0
0
0
0
0
18
0.082457
3.39
6
3
5
7
0
0
0
0
0
9
0
0
19
0.085299
3.24
6
3
7
0
0
0
0
0
0
9
0
5
20
0.087759
3.02
0
3
5
7
0
9
0
0
0
0
0
0
21
0.08863
2.77
0
3
5
7
0
0
0
0
0
9
0
0
22
0.089675
2.73
0
3
0
7
0
0
0
0
0
9
0
5
23
0.09589
2.59
9
3
7
0
0
0
0
0
0
0
0
5
24
0.103038
2.39
0
0
5
7
0
9
0
0
0
0
0
0
25
0.104241
2.14
0
0
5
7
0
0
0
0
0
9
0
0
26
0.108826
1.99
0
0
7
0
0
0
0
0
0
9
0
5
27
0.116404
1.97
0
3
7
0
0
0
0
0
0
9
0
0
28
0.119689
1.90
9
7
0
0
0
0
0
0
0
0
0
5
29
0.13725
1.70
0
0
0
7
0
9
0
0
0
0
0
0
30
0.139392
1.45
0
0
0
7
0
0
0
0
0
9
0
0
31
0.144902
1.34
0
0
7
0
0
0
0
0
0
9
0
0
32
0.15344
1.28
0
7
0
0
0
0
0
0
0
9
0
0
33
0.164821
1.25
9
7
0
0
0
0
0
0
0
0
0
0
34
0.24888
0.90
0
0
0
0
0
9
0
0
0
0
0
0
35
0.256016
0.65
0
0
0
0
0
0
0
0
0
9
0
0
36
0.289352
0.62
9
0
0
0
0
0
0
0
0
0
0
0Figure 5
Plots of PSs best performance in 30 runs of four algorithms on Scale1.


## 5.2.2. Comparison Study

In order to verify the applicability and feasibility of the proposed algorithm, we test MPACO, P-ACO, NSGA-II, and SPEA-II on Scale1, Scale 2, and Scale 3 instances. Three examples represent, respectively, small-, medium-, and large-scale BOWTA instances under the same condition. The quantitative indexes, such as Pareto figure and time performance, are chosen to evaluate the performance on C++. The value of three instances can be seen from Table 9.Table 9
The values of Scale1, Scale 2, and Scale 3 instances.
Scale1
Scale2
Scale3
Target
10
48
98
Missile
10
53
100In what follows, Scale1 instance is firstly taken as an example to illustrate the performance of MPACO, P-ACO, NSGA-II, and SPEA-II. Then, Scale 2 and Scale 3 instances are introduced to further demonstrate the performances of both algorithms. The statistical results are shown in Figures 5–7.Figure 6
Plots of PSs best performance in 30 runs of four algorithms on Scale2.Figure 7
Plots of PSs best performance in 30 runs of four algorithms on Scale3.Figures5–7 present the distribution of PSs obtained in 30 runs with the best performance of each algorithm for each test instance in objective space. When especially tested on Scale 2 and Scale 3, the MPACO can obtain the best solutions among the four algorithms. Due to increments of dynamic calculation approach of heuristic information and improved movement probability rules, the calculation accuracy of MPACO algorithm is improved. Boundary symmetric mutation strategy is also good for the algorithm to solve large-scale BOWTA problem.The primary aim of solving BOWTA problems is to lay the foundation for the DWTA problems whose goal is to provide a set of near-optimal or acceptable real-time decisions in air combat. So, in practice, the time performance of algorithms is also an important index. In the end, we test four algorithms on three instances in 30 runs and record the literation time of each algorithm. The statistical results of time performances are shown in Figures8–10.Figure 8
Boxplot of time performance on Scale1 (in secs).Figure 9
Boxplot of time performance on Scale2 (in secs).Figure 10
Boxplot of time performance on Scale3 (in secs).Although MPACO does not have time advantage in solving Scale1 and Scale 2 compared with other algorithms, it runs as fast as SPEA-II with the same number of function evaluations on Scale 3. Boundary symmetric mutation strategy improves the mutation efficiency and the quality of solutions and new evaporation factor is good for solving large-scale problems, but dynamic heuristic information and probability rules affect the iteration rate. In real air combat situations, pilots often make deadly decisions within seconds or even within milliseconds. Although the algorithm structure and implementation process are relatively complex and increase the computational complexity, high-performance hardware platform can be adopted to reduce the processing time and satisfy the real-time requirement.In this section, we do some work as follows. Firstly, we test the MPACO on a specific example. We verify the feasibility after comparing two kinds of PSs and show the corresponding distribution results. Secondly, to give a comparison study among the four algorithms, three different scale instances are introduced to evaluate each algorithm. Plots of the best PSs in 30 runs of each algorithm are presented to have an intuitive comparison. Lastly, we show the time performance of the four algorithms on each instance. To summarize, MPACO has been proven to be an effective technique for BOWTA optimization problem and obviously is the best among the four algorithms.


## 6. Conclusion

We apply ACO to WTA problem and propose the MPACO to solve BOWTA in this paper. The main contributions of the thesis are summarized as follows. On the one hand, the expected damage of enemy and the cost of missiles are taken into account from a practical viewpoint; on the other hand, pilot operation factor and missiles constraints are introduced to make the model in line with the real air combat. In this aspect, five strategies are proposed and embedded into traditional P-ACO algorithm to improve performances of MPACO algorithm. The experiments have shown that MPACO can find better Pareto solutions than the other three algorithms for all instances. More importantly, MPACO is more suitable for solving large-scale problems from the time performance viewpoint.However, we mainly have studied the static WTA problem; few researches focused on dynamic problems which are more instructive to real combat. In recent years, more and more researches began to pay attention to DWTA problems. A further study on this topic is one of our future tasks.



---

*Source: 1746124-2017-03-16.xml*