# A Hybrid Multiobjective Discrete Particle Swarm Optimization Algorithm for Cooperative Air Combat DWTA

Authors: Guang Peng; Yangwang Fang; Shaohua Chen; Weishi Peng; Dandan Yang
Journal: Journal of Optimization
 (2017)

---

## Abstract

A hybrid multiobjective discrete particle swarm optimization (HMODPSO) algorithm is proposed to solve cooperative air combat dynamic weapon target assignment (DWTA). First, based on the threshold of damage probability and time window constraints, a new cooperative air combat DWTA multiobjective optimization model is presented, which employs the maximum of the target damage efficiency and minimum of ammunition consumption as two competitive objective functions. Second, in order to tackle the DWTA problem, a mixed MODPSO and neighborhood search algorithm is proposed. Furthermore, the repairing operator is introduced into the mixed algorithm, which not only can repair infeasible solutions but also can improve the quality of feasible solutions. Besides, the Cauchy mutation is adopted to keep the diversity of the Pareto optimal solutions. Finally, a typical two-stage DWTA scenario is performed by HMODPSO and compared with three other state-of-the-art algorithms. Simulation results verify the effectiveness of the new model and the superiority of the proposed algorithm.


---

## Body

## 1. Introduction

The weapon target assignment (WTA) is a typical NP-complete constrained combinatorial optimization problem [1], which can be classified into two categories: static WTA (SWTA) and dynamic WTA (DWTA) [2, 3]. In SWTA, all the weapons attack targets in a single stage. In contrast, DWTA is much more complicated than SWTA, which takes the time window and resource constraints into account [4]. Besides, DWTA needs to deal with the new incoming targets and assesses the outcome of each engagement.Most of the previous researches on WTA are focused on SWTA [4]. However, DWTA has begun to gain more attention of researchers since it was put forward by Hosein and Athans in 1990 [3]. Cai et al. [4] provided a survey of the research on DWTA problem and introduced some basic concepts on DWTA. Khosla [5] proposed a hybrid approach, which combines genetic algorithm (GA) with simulated annealing (SA) to solve a target-based DWTA problem. Liu et al. [6] analyzed the time and space restriction of the mathematical model of DWTA and proposed an adaptive memetic algorithm to obtain the suboptimum solution step by step. Chen et al. [7] established a generic asset-based DWTA model which incorporates four categories of constraints, namely, capability constraints, strategy constraints, resource constraints, and engagement feasibility constraints. Based on the asset-based DWTA model, Xin et al. [8] proposed a new technique for constraint handling. Moreover, Xin et al. [9] proposed an efficient rule-based heuristic, which uses the domain knowledge of DWTA in the form of three crucial rules, to solve asset-based DWTA problems. The proposed method has obvious advantage over the Monte Carlo method (MCM) with regards to solution quality and computation time. Wang et al. [10] applied intuitionistic fuzzy entropy of discrete particle swarm optimization (IFDPSO) algorithm to solve DWTA problem. Wang et al. [11] firstly established a multicombat step DWTA game model of UAV aerial combat and then presented a clonal selection optimization algorithm to solve the model.However, the above WTA problems are focused on one objective (i.e., operational effects), ignoring the operational cost, while in actual combat situations, apart from considering the maximum of the damage to targets, the ammunition consumption should be also taken into account. Clearly, the two competitive objectives are conflicting, which implies that the WTA problem is a multiobjective optimization problem. Up to now, there have been several studies on multiobjective optimization for DWTA problems. Liu et al. [12] and Zhou et al. [13] proposed an improved MOPSO algorithm to solve the multiobjective programming model of SWTA, respectively. Lötter and Van Vuuren [14] used NSGA-II to solve a triobjective DWTA model for surface-based air defense. Li et al. [15, 16] established deterministic and uncertain multiobjective optimization models of multistage WTA (MWTA) problem and modified two multiobjective optimizers, NSGA-II and MOEA/D, by adding an adaptive mechanism for solving the NWTA models. However, MOPSO easily falls into the local optimum; NSGA-II and MOEA/D have complex computation for DWTA. DWTA problem has higher requirements for the real-time performance and the convergence. In order to meet the real-time performance and the convergence accuracy simultaneously, this paper proposed an efficient HMODPSO algorithm to solve the DWTA multiobjective optimization problem. The proposed HMODPSO algorithm can generate obviously better DWTA decisions without the cost of overmuch extra computation time, which can improve the cooperative air combat effectiveness.The rest of this paper is organized as follows. In Section2, the cooperative air combat DWTA multiobjective optimization model based on the threshold of damage probability is formulated. Section 3 presents the structure of HMODPSO algorithm. The simulation results based on the proposed algorithm for a typical two-stage DWTA scenario are discussed in Section 4, and also the comparisons with three other state-of-the-art algorithms are conducted in this section. Conclusions and future work will be drawn in Section 5.


## 2. The Cooperative Air Combat Multiobjective Optimization Model for DWTA

Definition 1 (time window of target).
The time window of target (tT) is the exposure time of a target which the weapon can attack efficiently.Definition 2 (time window of algorithm).
The time window of algorithm (tA) is the running time for solving the DWTA.A complete cooperative air combat is a multistage offensive and defensive process. So the DWTA model can be regarded as the repetition of the SWTA model with the damage assessment. There are at mostS stages of DWTA, which means no targets or weapons left after final stage assignment. In the cooperative air combat DWTA model, the “shoot-look-shoot” engagement policy is adopted. When entering the t+1 stage, it is necessary to determine a set of alive targets and remaining weapons. Therefore, it needs to observe the outcome of the t stage engagement and reformulate air combat situation assessment. The schematic diagram of cooperative air combat DWTA model is shown in Figure 1.Figure 1
The schematic diagram of cooperative air combat DWTA model.


### 2.1. The DWTA Multiobjective Optimization Model

Assuming in the cooperative air combat that there areF flights in blue formation and each flight carries Mf(f=1,2,…,F) missiles, the total number of weapons is M=∑f=1F‍Mf. At a certain time, the formation detects N targets, which can be attacked by the weapons. pij is the damage probability that the i-th missile attacks j-th target. wj is the threat value of target j. Obviously, pij and wj can be obtained by C4I system according to the weapons’ performance and air combat situation.In order to describe the DWTA problem, a Boolean type decision matrix is introduced:(1)X=x11x12⋯x1Nx21x22⋯x2N⋮xM1xM2⋯xMN,where xij=1 if the weapon i is assigned to the target j and xij=0 otherwise. After all the weapons attack the target j cooperatively, the joint damage probability of the target j can be expressed as(2)Pj=1-∏i=1M1-pijxij.In order to avoid wasting weapon resources in the single WTA optimization model, we define the maximum of the target damage efficiency and minimum of using weapon units as two objective functions. Additionally, aiming at targets with different threat values, the new model should guarantee the threshold of damage probability of each target. Hence, the formulation of the objective functions for the staget is constructed:(3)maxft=∑j=1Ntwjt1-∏i=1Mt1-pijtxijtmingt=∑j=1Nt∑i=1Mtxijt,where t is the stage index; M(t) and N(t) are the numbers of existing weapons and targets at the stage t, respectively; pij(t) is the damage probability that the i-th missile attacks the j-th target at the stage t; xij(t) is the Boolean type decision variable at the stage t.The following three categories of constraints are incorporated in the cooperative air combat DWTA multiobjective optimization model:(4)Pjt≥Pdjt,∀t∈1,2,…,S,∀j∈1,2,…,N(5)∑j=1Ntxijt≤1,∀t∈1,2,…,S,∀i∈1,2,…,M(6)tjTt≥tAt,∀t∈1,2,…,S,∀j∈1,2,…,N.Constraint set (4) represents the threshold of damage probability of each target, the value of which depends on the decision-makers or the command system based on the current air combat situation. If one of the targets does not satisfy the threshold, the weapon target assignment is regarded as invalid assignment. Constraint set (5) reflects the capability of weapons attacking targets at the same time. In fact, a missile can only shoot one target each time. Constraint set (6) is very important for DWTA model, which takes the time windows of target and algorithm into account, influencing the engagement feasibility of weapons. In this case, it makes particularly high requirements for operational efficiency of the algorithm.


### 2.2. Constraints Handling

The DWTA multiobjective optimization model is a typical constrained nonlinear combinational optimization problem. To ensure the feasibility of solutions generated by the proposed algorithm, we make some preliminary treatments to constraints (4)–(6). First, we utilize the penalty function method to handle the constraint set (4). The function is defined as(7)hjt=min⁡0,Pjt-Pdjt.For constraint set (5), it can be satisfied by solution encoding. And the constraint set (6) is satisfied by adopting the proposed algorithm.Above all, the multiobjective optimization model for the cooperative air combat DWTA problem aforementioned can be formulated as(8)maxft=∑j=1Ntwjt1-∏i=1Mt1-pijtxijt-σ∑j=1Nhjt2mingt=∑j=1Nt∑i=1Mtxijt,where σ is the penalty parameter and σ=1000 is selected in this paper.


## 2.1. The DWTA Multiobjective Optimization Model

Assuming in the cooperative air combat that there areF flights in blue formation and each flight carries Mf(f=1,2,…,F) missiles, the total number of weapons is M=∑f=1F‍Mf. At a certain time, the formation detects N targets, which can be attacked by the weapons. pij is the damage probability that the i-th missile attacks j-th target. wj is the threat value of target j. Obviously, pij and wj can be obtained by C4I system according to the weapons’ performance and air combat situation.In order to describe the DWTA problem, a Boolean type decision matrix is introduced:(1)X=x11x12⋯x1Nx21x22⋯x2N⋮xM1xM2⋯xMN,where xij=1 if the weapon i is assigned to the target j and xij=0 otherwise. After all the weapons attack the target j cooperatively, the joint damage probability of the target j can be expressed as(2)Pj=1-∏i=1M1-pijxij.In order to avoid wasting weapon resources in the single WTA optimization model, we define the maximum of the target damage efficiency and minimum of using weapon units as two objective functions. Additionally, aiming at targets with different threat values, the new model should guarantee the threshold of damage probability of each target. Hence, the formulation of the objective functions for the staget is constructed:(3)maxft=∑j=1Ntwjt1-∏i=1Mt1-pijtxijtmingt=∑j=1Nt∑i=1Mtxijt,where t is the stage index; M(t) and N(t) are the numbers of existing weapons and targets at the stage t, respectively; pij(t) is the damage probability that the i-th missile attacks the j-th target at the stage t; xij(t) is the Boolean type decision variable at the stage t.The following three categories of constraints are incorporated in the cooperative air combat DWTA multiobjective optimization model:(4)Pjt≥Pdjt,∀t∈1,2,…,S,∀j∈1,2,…,N(5)∑j=1Ntxijt≤1,∀t∈1,2,…,S,∀i∈1,2,…,M(6)tjTt≥tAt,∀t∈1,2,…,S,∀j∈1,2,…,N.Constraint set (4) represents the threshold of damage probability of each target, the value of which depends on the decision-makers or the command system based on the current air combat situation. If one of the targets does not satisfy the threshold, the weapon target assignment is regarded as invalid assignment. Constraint set (5) reflects the capability of weapons attacking targets at the same time. In fact, a missile can only shoot one target each time. Constraint set (6) is very important for DWTA model, which takes the time windows of target and algorithm into account, influencing the engagement feasibility of weapons. In this case, it makes particularly high requirements for operational efficiency of the algorithm.


## 2.2. Constraints Handling

The DWTA multiobjective optimization model is a typical constrained nonlinear combinational optimization problem. To ensure the feasibility of solutions generated by the proposed algorithm, we make some preliminary treatments to constraints (4)–(6). First, we utilize the penalty function method to handle the constraint set (4). The function is defined as(7)hjt=min⁡0,Pjt-Pdjt.For constraint set (5), it can be satisfied by solution encoding. And the constraint set (6) is satisfied by adopting the proposed algorithm.Above all, the multiobjective optimization model for the cooperative air combat DWTA problem aforementioned can be formulated as(8)maxft=∑j=1Ntwjt1-∏i=1Mt1-pijtxijt-σ∑j=1Nhjt2mingt=∑j=1Nt∑i=1Mtxijt,where σ is the penalty parameter and σ=1000 is selected in this paper.


## 3. HMODPSO Algorithm for DWTA

### 3.1. Particle Position Encoding

Decimal encoding for particles is adopted in this paper. The length of the particle position encoding (denoted byD) is the total number of weapons. Each weapon is treated as a dimension of particle position encoding, and the value of each dimension indicates the number of the targets to which the weapon is assigned. Figure 2 provides an example to explain the encoding method (M=6,N=4,t=1), and X is the corresponding 0-1 decision matrix. Obviously, such a particle position encoding method can guarantee that all the generated solutions satisfy the constraint set (5).Figure 2
Particle position encoding.


### 3.2. The Leader Particle Selecting

In PSO [17], each particle in the swarm corresponds to a potential solution of the optimization problem. In a D-dimensional search space, each particle has a position Xi(k)=[Xi,1(k),Xi,2(k),…,Xi,D(k)] and a velocity Vi(k)=[Vi,1(k),Vi,2(k),…,Vi,D(k)]. During the movement, the personal best position (Pbest) for each particle is recorded as Pi(k)=[Pi,1(k),Pi,2(k),…,Pi,D(k)], and the global best particle (Gbest) of the swarm is referred to as Gk=[G1(k),G2(k),…,GD(k)]. For specific DWTA problem, the velocity and position of the particle are updated by the following equations, respectively:(9)Vijk+1=wVijk+c1r1Pijk-Xijk+c2r2Gjk-Xijk(10)Xijk+1=roundXijk+Vijk+1(11)Vijk+1⟶max⁡min⁡Vijk+1,Vmax,Vmin(12)Xijk+1=max⁡min⁡Xijk+1,Xmax,Xmin,where i represents the i-th particle of the swarm, j represents the j-th dimension in the search space, k is the number of current iterations, w is the inertia weight, c1,c2 are the acceleration coefficients, and r1,r2∈[0,1] are uniformly distributed random variables. Xmax and Xmin are the lower and upper boundaries of the position, respectively; Vmax and Vmin are the lower and upper boundaries of the velocity, respectively. round is said to be an integer operator. In MOPSO, selecting the global best position or the leader (Gbest) randomly from the external archive is a popular way [12]. However, the roulette method may damage evolution direction of the particles. Considering the weakness, a new method of selecting the leader particle is presented. First, the square root distance (SRD) [18] of particles X1 and X2 is defined as follows:(13)SRDX1,X2=∑i=1mfiX1-fiX2,where m represents the number of the objective functions. During the process of iteration, first, calculate all the particles’ SRD with the nondominated solutions in the external archive; then, select the nondominated solution which has the minimum SRD as the leader particle; that is,(14)min∑i=1mfiY1-fiX,…,∑i=1LfiYL-fiX,where Y1,Y2,…,YL represents the L nondominated solutions in the external archive and X represents the current position of particle. Particles can choose closer nondominated solutions as their leader particles through the new method, which can help the algorithm overcome the uncertainty of random selection and improve the convergence accuracy.


### 3.3. Repairing Operator

Clearly, unfeasible solutions which do not satisfy the constraint set (4) may be generated for solving the DWTA. To cope with this issue, a repairing operator is introduced [19], which includes two parts.(
1) Deleting the Redundant Allocation. First, select the unfeasible solutions and mark the targets satisfying the threshold of the damage probability. Then, try to delete the missile attacking the target with minimum damage probability. After deleting, if the target still satisfies the threshold, confirm the deletion; otherwise, retain the missile.(
2) Supplementing the Insufficient Allocation. First, select the targets which do not satisfy threshold of damage probability. Then mark the target with the minimum joint damage probability, and select one missile attacking the target with the maximum damage probability. Repeat the operation until no weapons are left or all the targets satisfy the threshold of damage probability.The corresponding pseudocode of the repairing operator is shown in Algorithm1.Algorithm 1:The pseudocode of the repairing operator.
Input. Unfeasible solutions
Output. New solutions
(
1
)Deleting the redundant allocation
(1)for j=1:1:N
(2) Flagj≔i∣xij=1,i=1,2,…,M
(3) while Pj>Pdj
(4)  k≔argmini∈FlagPij
(5)  Pj′=(Pj-Pkj)/(1-Pkj)
(6)  if Pj′>Pdj
(7)    Flagj:xkj=0
(8)    Pj=Pj′
(9)  end if
(10) end while
(11)end for
(
2
)Supplementing the insufficient allocation
(12)  M′≔i∣∑‍xij=0,i=1,2,…,M
(13)  T′≔j∣Pj<Pdj,j=1,2,…,N
(14)  P′′=argminj∈T′Pj
(15)while M′>0&&P′′>Pdj
(16)   l≔argminj∈T′Pj
(17)   k≔argmaxi∈M′Pil
(18)   xkl=1
(19)   M′=M′-1
(20)  P′′=argminj∈T′Pj
(21)end while


### 3.4. Cauchy Mutation Operator

In order to further maintain the diversity of particles, the Cauchy mutation [20] is introduced into HMODPSO algorithm. When the current particle is dominated by Pbest, the Cauchy mutation is applied to disturb the current particle. The Cauchy mutation operator is defined as(15)Xij=Xij+roundXmax-Xmin×Cauchy0,1,rand≤pbXij,rand>pb(16)Cauchy0,1=tan⁡rand-0.5×π(17)Xij=Xmin,Xij<XminXmax,Xij>XmaxXij,elsewhile,where rand∈[0,1] is a uniformly distributed random value, Cauchy(0,1) is a standard Cauchy-distributed random value, and the adapting mutation rate Pb is computed as(18)Pb=1-tmax⁡It.


### 3.5. Neighborhood Searching Operator

Neighborhood search (NS) algorithm begins with an initial solution and searches the better solution in its neighborhood range [21]. At present, define a solution X=1,2,3,4, and a new vector can be obtained by exchanging any two positions of the solution. A collection of all the new vectors is referred to as exchange neighborhood shown in Figure 3.Figure 3
Neighborhood exchange criterion.NS algorithm is usually carried out at the end of the global search for individual local optimization, which can improve the convergence accuracy of algorithm. In MOPSO, the nondominated solutions in the external archive are regarded as leader particles to guide the swarm to fly. Additionally, the nondominated solutions are evolved into the Pareto optimal solutions set. In view of the particularity of the external archive, NS algorithm is introduced into the external archive. This paper proposed three kinds of operations: NS-1, NS-2, and NS-3. The procedures are described as follows.For NS-1, do the following steps.Step 1.
Take each nondominated solution in the external archive as the initial solution.Step 2.
Exchange randomly two positions of the initial solution to get a new solution. If the new solution is better than the initial solution, then replace it with the new solution; otherwise, keep the initial solution.For NS-2, do the following steps.Step 1.
Take each nondominated solution in the external archive as the initial solution.Step 2.
Select the best solution from the neighborhood range of the initial solution.Step 3.
If the best solution is better than the initial solution, then replace it with the new solution; otherwise, keep the initial solution.For NS-3, do the following steps.Step 1.
Take each nondominated solution in the external archive as the initial solution.Step 2.
Exchange two positions of the initial solution in sequence to get a new solution; if the new solution is better than the initial solution, replace it with the new solution; then serve the new solution as the initial solution and repeat the above operation; otherwise, keep the initial solution.The corresponding pseudocode of the local search operations is shown in Algorithm2.Algorithm 2:The pseudocode of the local search operations.
NS-1
(
1
) for i=1:1:L
(
2
)  NX=crossover(X,X)
(
3
)  evaluatefitness(NX)
(
4
)  X=Compare(X,NX)
(
5
) end for
NS-2
(
1
) for i=1:1:L
(
2
)  for j=1:1:D-1
(
3
)   for k=j+1:1:D
(
4
)   NX=exchange(X(j,k))
(
5
)   evaluatefitness(NX)
(
6
)   end for
(
7
)  end for
(
8
)  X=Compare(X,max(NX))
(
9
) end for
NS-3
(
1
) for i=1:1:L
(
2
)  for j=1:1:D-1
(
3
)   for k=j+1:1:D
(
4
)   NX=exchange(X(j,k))
(
5
)   evaluatefitness(NX)
(
6
)   X=Compare(X,NX)
(
7
)   end for
(
8
)  end for
(
9
) end for


### 3.6. The Procedure of HMODPSO Algorithm

Summarizing the above procedures, we obtain the following pseudocode of the HMODPSO algorithm in Algorithm3, whereCheckBoundaries validates the particles to search in the solution space and max⁡It is the maximum iterations. Furthermore, use the crowding distance sorting strategy [22] to maintain the external archive. Figure 4 illustrates the entire flowchart of HMODPSO.Algorithm 3:The pseudocode of the HMODPSO algorithm.
(
1
) Initialize the swarm with H particles
(
2
) Initialize the external archive
(
3
) for k=1:1:max⁡It
(
4
) for j=1:1:H
(
5
)  Select leader for each particle
(
6
)  Update velocity (Eq. (9))
(
7
)  Update position (Eq. (10))
(
8
)  CheckBoundaries (Eq. (12))
(
9
)  Repairing operator
(
10
)  Evaluate the particle
(
11
)  Update Pbest
(
12
)  Cauchy mutation operator (Eq. (15))
(
13
)  CheckBoundaries (Eq. (17))
(
14
)  Update Pbest
(
15
)  Update external archive
(
16
)  Neighborhood searching operator
(
17
) end for
(
18
) if archive_size>max_archive
(
19
)  Maintain the external archive by crowding distance sorting strategy
(
20
) end if
(
21
) end for
(
22
) Return the external archiveFigure 4
The entire flowchart of HMODPSO.


## 3.1. Particle Position Encoding

Decimal encoding for particles is adopted in this paper. The length of the particle position encoding (denoted byD) is the total number of weapons. Each weapon is treated as a dimension of particle position encoding, and the value of each dimension indicates the number of the targets to which the weapon is assigned. Figure 2 provides an example to explain the encoding method (M=6,N=4,t=1), and X is the corresponding 0-1 decision matrix. Obviously, such a particle position encoding method can guarantee that all the generated solutions satisfy the constraint set (5).Figure 2
Particle position encoding.


## 3.2. The Leader Particle Selecting

In PSO [17], each particle in the swarm corresponds to a potential solution of the optimization problem. In a D-dimensional search space, each particle has a position Xi(k)=[Xi,1(k),Xi,2(k),…,Xi,D(k)] and a velocity Vi(k)=[Vi,1(k),Vi,2(k),…,Vi,D(k)]. During the movement, the personal best position (Pbest) for each particle is recorded as Pi(k)=[Pi,1(k),Pi,2(k),…,Pi,D(k)], and the global best particle (Gbest) of the swarm is referred to as Gk=[G1(k),G2(k),…,GD(k)]. For specific DWTA problem, the velocity and position of the particle are updated by the following equations, respectively:(9)Vijk+1=wVijk+c1r1Pijk-Xijk+c2r2Gjk-Xijk(10)Xijk+1=roundXijk+Vijk+1(11)Vijk+1⟶max⁡min⁡Vijk+1,Vmax,Vmin(12)Xijk+1=max⁡min⁡Xijk+1,Xmax,Xmin,where i represents the i-th particle of the swarm, j represents the j-th dimension in the search space, k is the number of current iterations, w is the inertia weight, c1,c2 are the acceleration coefficients, and r1,r2∈[0,1] are uniformly distributed random variables. Xmax and Xmin are the lower and upper boundaries of the position, respectively; Vmax and Vmin are the lower and upper boundaries of the velocity, respectively. round is said to be an integer operator. In MOPSO, selecting the global best position or the leader (Gbest) randomly from the external archive is a popular way [12]. However, the roulette method may damage evolution direction of the particles. Considering the weakness, a new method of selecting the leader particle is presented. First, the square root distance (SRD) [18] of particles X1 and X2 is defined as follows:(13)SRDX1,X2=∑i=1mfiX1-fiX2,where m represents the number of the objective functions. During the process of iteration, first, calculate all the particles’ SRD with the nondominated solutions in the external archive; then, select the nondominated solution which has the minimum SRD as the leader particle; that is,(14)min∑i=1mfiY1-fiX,…,∑i=1LfiYL-fiX,where Y1,Y2,…,YL represents the L nondominated solutions in the external archive and X represents the current position of particle. Particles can choose closer nondominated solutions as their leader particles through the new method, which can help the algorithm overcome the uncertainty of random selection and improve the convergence accuracy.


## 3.3. Repairing Operator

Clearly, unfeasible solutions which do not satisfy the constraint set (4) may be generated for solving the DWTA. To cope with this issue, a repairing operator is introduced [19], which includes two parts.(
1) Deleting the Redundant Allocation. First, select the unfeasible solutions and mark the targets satisfying the threshold of the damage probability. Then, try to delete the missile attacking the target with minimum damage probability. After deleting, if the target still satisfies the threshold, confirm the deletion; otherwise, retain the missile.(
2) Supplementing the Insufficient Allocation. First, select the targets which do not satisfy threshold of damage probability. Then mark the target with the minimum joint damage probability, and select one missile attacking the target with the maximum damage probability. Repeat the operation until no weapons are left or all the targets satisfy the threshold of damage probability.The corresponding pseudocode of the repairing operator is shown in Algorithm1.Algorithm 1:The pseudocode of the repairing operator.
Input. Unfeasible solutions
Output. New solutions
(
1
)Deleting the redundant allocation
(1)for j=1:1:N
(2) Flagj≔i∣xij=1,i=1,2,…,M
(3) while Pj>Pdj
(4)  k≔argmini∈FlagPij
(5)  Pj′=(Pj-Pkj)/(1-Pkj)
(6)  if Pj′>Pdj
(7)    Flagj:xkj=0
(8)    Pj=Pj′
(9)  end if
(10) end while
(11)end for
(
2
)Supplementing the insufficient allocation
(12)  M′≔i∣∑‍xij=0,i=1,2,…,M
(13)  T′≔j∣Pj<Pdj,j=1,2,…,N
(14)  P′′=argminj∈T′Pj
(15)while M′>0&&P′′>Pdj
(16)   l≔argminj∈T′Pj
(17)   k≔argmaxi∈M′Pil
(18)   xkl=1
(19)   M′=M′-1
(20)  P′′=argminj∈T′Pj
(21)end while


## 3.4. Cauchy Mutation Operator

In order to further maintain the diversity of particles, the Cauchy mutation [20] is introduced into HMODPSO algorithm. When the current particle is dominated by Pbest, the Cauchy mutation is applied to disturb the current particle. The Cauchy mutation operator is defined as(15)Xij=Xij+roundXmax-Xmin×Cauchy0,1,rand≤pbXij,rand>pb(16)Cauchy0,1=tan⁡rand-0.5×π(17)Xij=Xmin,Xij<XminXmax,Xij>XmaxXij,elsewhile,where rand∈[0,1] is a uniformly distributed random value, Cauchy(0,1) is a standard Cauchy-distributed random value, and the adapting mutation rate Pb is computed as(18)Pb=1-tmax⁡It.


## 3.5. Neighborhood Searching Operator

Neighborhood search (NS) algorithm begins with an initial solution and searches the better solution in its neighborhood range [21]. At present, define a solution X=1,2,3,4, and a new vector can be obtained by exchanging any two positions of the solution. A collection of all the new vectors is referred to as exchange neighborhood shown in Figure 3.Figure 3
Neighborhood exchange criterion.NS algorithm is usually carried out at the end of the global search for individual local optimization, which can improve the convergence accuracy of algorithm. In MOPSO, the nondominated solutions in the external archive are regarded as leader particles to guide the swarm to fly. Additionally, the nondominated solutions are evolved into the Pareto optimal solutions set. In view of the particularity of the external archive, NS algorithm is introduced into the external archive. This paper proposed three kinds of operations: NS-1, NS-2, and NS-3. The procedures are described as follows.For NS-1, do the following steps.Step 1.
Take each nondominated solution in the external archive as the initial solution.Step 2.
Exchange randomly two positions of the initial solution to get a new solution. If the new solution is better than the initial solution, then replace it with the new solution; otherwise, keep the initial solution.For NS-2, do the following steps.Step 1.
Take each nondominated solution in the external archive as the initial solution.Step 2.
Select the best solution from the neighborhood range of the initial solution.Step 3.
If the best solution is better than the initial solution, then replace it with the new solution; otherwise, keep the initial solution.For NS-3, do the following steps.Step 1.
Take each nondominated solution in the external archive as the initial solution.Step 2.
Exchange two positions of the initial solution in sequence to get a new solution; if the new solution is better than the initial solution, replace it with the new solution; then serve the new solution as the initial solution and repeat the above operation; otherwise, keep the initial solution.The corresponding pseudocode of the local search operations is shown in Algorithm2.Algorithm 2:The pseudocode of the local search operations.
NS-1
(
1
) for i=1:1:L
(
2
)  NX=crossover(X,X)
(
3
)  evaluatefitness(NX)
(
4
)  X=Compare(X,NX)
(
5
) end for
NS-2
(
1
) for i=1:1:L
(
2
)  for j=1:1:D-1
(
3
)   for k=j+1:1:D
(
4
)   NX=exchange(X(j,k))
(
5
)   evaluatefitness(NX)
(
6
)   end for
(
7
)  end for
(
8
)  X=Compare(X,max(NX))
(
9
) end for
NS-3
(
1
) for i=1:1:L
(
2
)  for j=1:1:D-1
(
3
)   for k=j+1:1:D
(
4
)   NX=exchange(X(j,k))
(
5
)   evaluatefitness(NX)
(
6
)   X=Compare(X,NX)
(
7
)   end for
(
8
)  end for
(
9
) end for


## 3.6. The Procedure of HMODPSO Algorithm

Summarizing the above procedures, we obtain the following pseudocode of the HMODPSO algorithm in Algorithm3, whereCheckBoundaries validates the particles to search in the solution space and max⁡It is the maximum iterations. Furthermore, use the crowding distance sorting strategy [22] to maintain the external archive. Figure 4 illustrates the entire flowchart of HMODPSO.Algorithm 3:The pseudocode of the HMODPSO algorithm.
(
1
) Initialize the swarm with H particles
(
2
) Initialize the external archive
(
3
) for k=1:1:max⁡It
(
4
) for j=1:1:H
(
5
)  Select leader for each particle
(
6
)  Update velocity (Eq. (9))
(
7
)  Update position (Eq. (10))
(
8
)  CheckBoundaries (Eq. (12))
(
9
)  Repairing operator
(
10
)  Evaluate the particle
(
11
)  Update Pbest
(
12
)  Cauchy mutation operator (Eq. (15))
(
13
)  CheckBoundaries (Eq. (17))
(
14
)  Update Pbest
(
15
)  Update external archive
(
16
)  Neighborhood searching operator
(
17
) end for
(
18
) if archive_size>max_archive
(
19
)  Maintain the external archive by crowding distance sorting strategy
(
20
) end if
(
21
) end for
(
22
) Return the external archiveFigure 4
The entire flowchart of HMODPSO.


## 4. Simulations and Results

In the experimental scenario, setF=5, Mf=4(f=1,2,3,4,5), and N=10, and obtain M=20. When t=1, the time window of each target tjT(j=1,2,…,10) is [8.3,11.4,16.2,13.8,20.1,25.4,19.6,13.5,17.2,22.8]; and each target’s threat coefficient wj(j=1,2,…,10) is [0.6,0.7,0.3,0.5,0.6,0.35,0.65,0.55,0.4,0.75]. The damage probability threshold of each target can be set to 0.9. The weapon’s damage probability pij(i=1,2,…,20,j=1,2,…,10) is(19)P1=0.540.830.900.880.750.660.860.810.770.630.730.800.780.650.560.770.820.840.710.850.750.820.900.770.900.790.810.650.590.820.680.730.710.870.760.640.830.640.810.840.810.750.870.530.800.900.610.800.730.810.850.800.910.850.850.840.590.770.640.660.700.720.730.890.770.670.850.680.850.830.810.830.620.750.830.800.750.680.600.760.830.590.550.910.560.730.630.880.770.650.850.840.590.770.860.850.780.660.610.770.800.730.610.900.890.660.780.910.840.570.670.640.850.880.560.730.890.770.920.660.830.850.890.810.720.730.890.770.920.660.840.890.770.880.570.770.820.790.810.650.860.850.880.660.850.880.560.730.630.840.660.650.830.820.600.710.650.840.830.580.650.710.810.560.770.860.850.710.660.610.890.860.780.660.850.800.640.850.850.840.560.730.630.880.590.770.860.850.710.870.900.670.820.680.610.800.890.880.780.66.To verify the efficiency of HMODPSO algorithm, three different kinds of algorithms were proposed: HMODPSO-1 (adopting NS-1 operation), HMODPSO-2 (adopting NS-2 operation), and HMODPSO-3 (adopting NS-3 operation). At the same time, compare the proposed algorithms with NSGA-II [23], MODPSO [12], and MODPSO-GSA [24] to show their potential competences. In the experiment, population size is 60, 100 iterations are carried out, the external archive size is chosen as 30, the crossover probability of NSGA-II is 0.8, and the mutation rate of NSGA-II is 0.1. All the simulations were performed under the same environment (Matlab) on Intel Core i5-4590 3.3 GHz CPU with 4 GB RAM.The Pareto fronts produced by six algorithms are shown in Figure5; it can be seen that all the solutions using HMODPSO are better than the other three algorithms. MODPSO-GSA can also get good solutions. However, MODPSO and NSGA-II easily fall into the local optimum. The different Pareto fronts are compared in Figure 6.Figure 5
The Pareto fronts obtained by six algorithms.
(a)
NSGA-II
(b)
MODPSO
(c)
MODPSO-GSA
(d)
HMODPSO-1
(e)
HMODPSO-2
(f)
HMODPSO-3Figure 6
The compared Pareto fronts.For each algorithm, 30 independent runs were executed. The computational time of algorithm is represented byT. The average results of the Pareto optimal solutions attained by six algorithms are shown in Table 1.Table 1
The average results of the Pareto optimal solutions.
Pareto solutions
NSGA-II
MODPSO
MODPSO-GSA
HMODPSO-1
HMODPSO-2
HMODPSO-3
13
∗
∗
5.0500
5.0524
5.0524
5.0524
14
∗
∗
5.0779
5.0998
5.1052
5.1054
15
5.0551
5.0550
5.1516
5.1510
5.1536
5.1576
16
5.1220
5.0635
5.1971
5.1975
5.1976
5.1980
17
5.1346
5.1306
5.2380
5.2372
5.2382
5.2291
18
5.2158
5.2031
5.2682
5.2639
5.2680
5.2682
19
5.2375
5.2283
5.2894
5.2876
5.2899
5.2904
20
5.2637
5.2614
5.3038
5.3095
5.3135
5.3183
T/s
56.23
3.1364
6.8735
4.4328
8.7426
12.5279Owing to the time window constraint of the DWTA multiobjective optimization model, the computational time of algorithm must satisfy requirements of the targets’ time window. From Table1, NSGA-II, HMODPSO-2, and HMODPSO-3 algorithms do not meet the requirements. MODPSO has the fastest convergence speed, but it easily falls into the local optimum. Compared with MODPSO-GSA, HMODPSO-1 has better comprehensive performance. So, in the t=1 stage, the assignment (f=5.0524,g=13) which was obtained once by HMODPSO-1 can be selected for the engagement. Figure 7 shows the damage probability of each target and the scheme of the assignment (f=5.0524,g=13).Figure 7
The WTA scheme obtained by HMODPSO-1.
(a)
The damage probability
(b)
The schemeAfter the first engagement, the parametersM(t+1), N(t+1), wj(t+1), pij(t+1), and tjT(t+1) need to be updated in the next stage according to the target damage assessment and the new air combat situation. Assume that, in the t=2 stage, M(2)=7 and N(2)=4. The damage probability threshold of each target can be set to 0.9. The time window of each target tjT(2) is [6.5,8.1,5.9,0.2]; and each target’s threat coefficient wj(2) is [0.6,0.8,0.5,0.7]. The weapon’s damage probability pij(2) is(20)P2=0.650.740.850.810.720.770.920.590.860.630.710.780.530.800.730.650.900.620.790.740.820.910.750.840.560.650.730.94.Sincet4T(2)=0.2s is so small that it is difficult to satisfy the time window of the special target, a new priority selecting method is efficiently adopted to deal with several special targets whose time windows do not satisfy the time window constraint of the proposed algorithm. From the existing weapons at a certain stage, select the missile with the maximum damage probability to attack the special target; repeat this operation until the special target’s joint damage probability satisfies the threshold. As seen from P(2), the 7-th missile attacking the target can satisfy the threshold of damage probability. The assignment of the remaining targets can be solved by HMODPSO-1 algorithm. Figure 8 shows the Pareto front.Figure 8
The Pareto front obtained by HMODPSO-1.HMODPSO-1 was independently run for 30 times. The average results of the Pareto optimal solutions attained by HMODPSO-1 are shown in Table2.Table 2
The average results of the Pareto optimal solutions.
Pareto solutions
HMODPSO-1
3
1.7280
4
1.7856
5
1.8372
6
1.8712
T
/
s
1.1155In thet=2 stage, the assignment (f=1.7280,g=3) can be selected for the engagement. Figure 9 shows the damage probability of each target and the scheme of the assignment (f=1.7280,g=3). HMODPSO-1 can efficiently solve the DWTA until no targets or weapons are left in the cooperative air combat.Figure 9
The WTA scheme obtained by HMODPSO-1.
(a)
The damage probability
(b)
The scheme


## 5. Conclusions

This paper has presented a new hybrid multiobjective DPSO called HMODPSO algorithm to solve the cooperative air combat DWTA multiobjective optimization problems. The proposed algorithm has three advantages. (a) The leader particle selecting operator and neighborhood searching operator can improve the search ability and the rate of convergence. (b) The repairing operator can enhance the efficiency of generating feasible solutions. (c) The Cauchy mutation operator can boost the diversity and distribution of Pareto optimal solutions. HMODPSO can find solutions with good accuracy, convergence, and diversity, which can generate obviously better DWTA decisions without the cost of overmuch extra computation time compared to NSGA-II, MODPSO, MODPSO-GSA, and MOEA/D. HMODPSO can also improve the cooperative air combat effectiveness efficiently.Future research will focus on optimizing DWTA instances under larger scales. In addition, more new mechanisms on reducing the time complexity of the proposed algorithm will also be investigated.



---

*Source: 8063767-2017-04-11.xml*