# MOQPSO-D/S for Air and Missile Defense WTA Problem under Uncertainty

Authors: Hao Xu; Qinghua Xing; Zhenhao Tian
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Aiming at the shortcomings of single objective optimization for solving weapon target assignment (WTA) and the existing multiobjective optimization based WTA method having problems being applied in air and missile defense combat under uncertainty, a fuzzy multiobjective programming based WTA method was proposed to enhance the adaptability of WTA decision to the changes of battlefield situation. Firstly, a multiobjective quantum-behaved particle swarm optimization with double/single-well (MOQPSO-D/S) algorithm was proposed by adopting the double/single-well based position update method, the hybrid random mutation method, and the two-stage based guider particles selection method. Secondly, a fuzzy multiobjective programming WTA model was constructed with consideration of air and missile defense combat’s characteristics. And, the uncertain WTA model was equivalently clarified based on the necessity degree principle of uncertainty theory. Thirdly, with particles encoding and illegal particles adjusting, the MOQPSO-D/S algorithm was adopted to solve the fuzzy multiobjective programming based WTA model. Finally, example simulation was conducted, and the result shows that the WTA model constructed is rational and MOQPSO-D/S algorithm is efficient.


---

## Body

## 1. Introduction

In information warfare, air attack is a highly integrated operation form. Defenders need to carry out the task of antiaerodynamic targets and antiballistic missiles simultaneously [1]. As one of the key phases in air and missile defense combat, weapon target assignment (WTA) plays a significant role on operational effectiveness. Therefore, study on air and missile defense WTA problem is of great significance.Given the importance, WTA has been studied with application of Lagrange relaxation algorithm [2], ant colony algorithm [3], genetic algorithm [4], clone selection algorithm [5], and particle swarm optimization algorithm [6]. Most of these researches construct WTA model based on single objective optimization with the objective function of maximizing the kill probability to enemy. They need to allocate all the firepower without consideration of the sequent operation need and cannot fit with the operational practice. In order to avoid allocating all firepower, a WTA model was constructed with the objective function of maximizing effectiveness-cost ratio in literature [6]. However, the adaptability of WTA decision, based on this method, to the changes of battlefield situation is still weak. Because this method is still a single objective optimization method, when conducting this method, the decision-making preferences need to be provided beforehand and only one WTA alternative can be given. Aiming at overcoming these weak points, the multiobjective optimization method was introduced into studying WTA problem [7–9]. With consideration of every objective, the Pareto solution set, obtained by multiobjective optimization method, can provide a series of WTA alternatives for the commander. And the final WTA decision will be made by choosing one from the alternatives. These studies [7–9] could provide useful reference for research on air and missile defense WTA problem. However, they cannot be applied to air and missile defense combat, which has the characteristic of multiweapon composite disposition, directly. Moreover, with the growing complexity of battlefield, the certain WTA method [3–9] can hardly be applied to air and missile defense combat under uncertainty.In order to address the issues above, the WTA model for air and missile defense with uncertain threat information is constructed based on fuzzy multiobjective programming, with consideration of the air and missile defense combat’s characteristics. Due to the superior performance of quantum-behaved particle swarm optimization algorithm, a multiobjective quantum-behaved particle swarm optimization with double/single-well algorithm is proposed for solving the WTA model.


## 2. Basic Theory

### 2.1. Concepts of Multiobjective Optimization

Definition 1 (multiobjective optimization problem).
Taking the minimizing problem as an example, a multiobjective optimization problem can be defined as follows:(1)minfx=f1x,f2x,…,foxs.t.x∈Ω.In (1), fl(x)  (l=1,2,…,o) is objective function, x is the decision vector, and Ω is the decision space defined by a series of constraints.Definition 2 (Pareto dominance).
Letxa and xb be two feasible solutions in Ω; if(2)allflxa≤flxb&&  anyflxa<flxbl=1,2,…,o,then xa Pareto dominates xb, denoted by xa≺xb.Definition 3 (Pareto optimal solution).
Letxa∈Ω; if ¬∃xb∈Ω, meeting xb≺xa, then xa would be defined as a Pareto optimal solution. And the set Ωp={xa∈Ω∣¬∃xb∈Ω,xb≺xa} would be defined as Pareto optimal solution set.Definition 4 (Pareto front).
The region which consisted of the objective values, corresponding to all solutions inΩp, is defined as Pareto front.


### 2.2. QPSO Algorithm

Quantum-behaved particle swarm optimization (QPSO) algorithm [10, 11], compared with the particle swarm optimization algorithm, has many advantages, such as faster convergence rate, fewer control parameters, and better global convergence. It has attracted much attention of scholars [12, 13].QPSO algorithm is a quantum mechanics based optimization algorithm. It is assumed that the motion state of the particles in the optimal space relative to the attractorP can be described by the wave function ψ(Y) in δ potential well. And ψ(Y) is described as follows:(3)ψY=1Le-Y/L.In (3), L is the characteristic length of δ potential well.The corresponding probability density function isQY=1/Le-2Y/L, and the probability distribution function is FY=e-2Y/L.Supposeu~U(0,1), and let u=F(Y); the position of the particle relative to the attractor P can be obtained as Y=±L/2ln1/u.And the absolute position of the particle in the optimal space isx=P+Y. Therefore, the position update formula of particle xk can be obtained as follows:(4)xkt+1=Pkt±Lkt2ln1ukt.In (4),t is the current generation, uk(t)~U(0,1). If uk(t)>0.5, “±” would be chosen as “+”, or it would be chosen as “−”. Pk(t) and Lk(t) are as follows:(5)Pkt=φktpbestkt+1-φktgbestt,(6)Lkt=2αPkt-xkt.In (5), φk(t)~U(0,1), pbestk(t) is the individual best position of particle xk, and gbest(t) is the global best position. In (6), α is the expansion-constriction factor.


### 2.3. QPSO with Double-Well Algorithm

Based on the double-δ potential well quantum model, a QPSO with double-well algorithm is proposed by Xu et al. [14] to resolve the problem of population diversity decline due to the fast convergence of QPSO algorithm. In this algorithm, the motion of particles is simulated as in a double-well space. The double-well based position update formula of particle xk can be defined as follows:(7)xkt+1=121+βktPk1t+1-βktPk2t±Lkt2ln1ukt.In (7), βk(t) is the distance coefficient between two δ potential wells, and βk(t) is a random number uniformly distributing over [0,1]. If uk(t)>0.5, “±” would be chosen as “+”, or it would be chosen as “−”. The two attractors Pk1t and Pk2t are defined as (8), and the characteristic length Lkt is defined as (9).(8)Pk1t=βk1tpbestkt+1-βk1tgbestk1t,Pk2t=βk2tpbestkt+1-βk2tgbestk2t(9)Lkt=α1+βktPk1t+1-βktPk2t-2xkt.In (8), βk1t~U(0,1), βk2t~U(0,1), pbestkt is the individual best position, and gbestk1t and gbestk2t are two guider particles.


## 2.1. Concepts of Multiobjective Optimization

Definition 1 (multiobjective optimization problem).
Taking the minimizing problem as an example, a multiobjective optimization problem can be defined as follows:(1)minfx=f1x,f2x,…,foxs.t.x∈Ω.In (1), fl(x)  (l=1,2,…,o) is objective function, x is the decision vector, and Ω is the decision space defined by a series of constraints.Definition 2 (Pareto dominance).
Letxa and xb be two feasible solutions in Ω; if(2)allflxa≤flxb&&  anyflxa<flxbl=1,2,…,o,then xa Pareto dominates xb, denoted by xa≺xb.Definition 3 (Pareto optimal solution).
Letxa∈Ω; if ¬∃xb∈Ω, meeting xb≺xa, then xa would be defined as a Pareto optimal solution. And the set Ωp={xa∈Ω∣¬∃xb∈Ω,xb≺xa} would be defined as Pareto optimal solution set.Definition 4 (Pareto front).
The region which consisted of the objective values, corresponding to all solutions inΩp, is defined as Pareto front.


## 2.2. QPSO Algorithm

Quantum-behaved particle swarm optimization (QPSO) algorithm [10, 11], compared with the particle swarm optimization algorithm, has many advantages, such as faster convergence rate, fewer control parameters, and better global convergence. It has attracted much attention of scholars [12, 13].QPSO algorithm is a quantum mechanics based optimization algorithm. It is assumed that the motion state of the particles in the optimal space relative to the attractorP can be described by the wave function ψ(Y) in δ potential well. And ψ(Y) is described as follows:(3)ψY=1Le-Y/L.In (3), L is the characteristic length of δ potential well.The corresponding probability density function isQY=1/Le-2Y/L, and the probability distribution function is FY=e-2Y/L.Supposeu~U(0,1), and let u=F(Y); the position of the particle relative to the attractor P can be obtained as Y=±L/2ln1/u.And the absolute position of the particle in the optimal space isx=P+Y. Therefore, the position update formula of particle xk can be obtained as follows:(4)xkt+1=Pkt±Lkt2ln1ukt.In (4),t is the current generation, uk(t)~U(0,1). If uk(t)>0.5, “±” would be chosen as “+”, or it would be chosen as “−”. Pk(t) and Lk(t) are as follows:(5)Pkt=φktpbestkt+1-φktgbestt,(6)Lkt=2αPkt-xkt.In (5), φk(t)~U(0,1), pbestk(t) is the individual best position of particle xk, and gbest(t) is the global best position. In (6), α is the expansion-constriction factor.


## 2.3. QPSO with Double-Well Algorithm

Based on the double-δ potential well quantum model, a QPSO with double-well algorithm is proposed by Xu et al. [14] to resolve the problem of population diversity decline due to the fast convergence of QPSO algorithm. In this algorithm, the motion of particles is simulated as in a double-well space. The double-well based position update formula of particle xk can be defined as follows:(7)xkt+1=121+βktPk1t+1-βktPk2t±Lkt2ln1ukt.In (7), βk(t) is the distance coefficient between two δ potential wells, and βk(t) is a random number uniformly distributing over [0,1]. If uk(t)>0.5, “±” would be chosen as “+”, or it would be chosen as “−”. The two attractors Pk1t and Pk2t are defined as (8), and the characteristic length Lkt is defined as (9).(8)Pk1t=βk1tpbestkt+1-βk1tgbestk1t,Pk2t=βk2tpbestkt+1-βk2tgbestk2t(9)Lkt=α1+βktPk1t+1-βktPk2t-2xkt.In (8), βk1t~U(0,1), βk2t~U(0,1), pbestkt is the individual best position, and gbestk1t and gbestk2t are two guider particles.


## 3. MOQPSO-D/S Algorithm

Given the superior performance, the QPSO algorithm is introduced to solve the multiobjective problems. And the multiobjective quantum-behaved particle swarm optimization algorithm (MOQPSO) [14–18] is put forward. However, the requirement for solving a multiobjective optimization problem is different from solving a single objective optimization problem. It requires that the diversity of solutions should be good, the distribution of solutions should be uniform, and the distance between solution and the true Pareto front should be as close as possible. The fast convergence performance of QPSO would lead to the premature convergence of MOQPSO. And the global convergence to an optimal solution of QPSO would go against the diversity and distribution of solutions obtained by MOQPSO. In order to improve the performance of MOQPSO algorithm to solve the multiobjective programming problem, a multiobjective quantum-behaved particle swarm optimization with double/single-well (MOQPSO-D/S) algorithm was proposed by combining the QPSO algorithm and the QPSO with double-well algorithm. In this algorithm, a double/single-well based position update method is adopted to balance the relationship of solutions’ diversity, convergence precision, and convergence rate. A hybrid random mutation method is also adopted to avoid premature convergence and improve the convergence precision. A two-stage based guider particles selection method is applied to this algorithm to improve the solution distribution uniformity and convergence rate.


### 3.1. Double/Single-Well Based Position Update Method

The so-called double/single-well based position update method includes two aspects. On one hand, in order to improve solutions’ diversity and avoid premature convergence, the particle position is updated based on QPSO with double-well algorithm in the early optimization stage. On the other hand, in order to improve convergence precision and convergence rate, the particle position is updated based on QPSO algorithm in the final optimization stage. The double/single-well based position update formula is as follows:(10)xkt+1=121+βktPk1t+1-βktPk2t±Lkt2ln1ukt,t≤tmax2Pk1t±Lkt2ln1ukt,t>tmax2.In (10), tmax is the maximum generation. If uk(t)>0.5, “±” would be chosen as “+”, or it would be chosen as “−”. The characteristic length Lk(t) is defined as in the following equation.(11)Lkt=α1+βktPk1t+1-βktPk2t-2xkt,t≤tmax22αPk1t-xkt,t>tmax2.


### 3.2. Hybrid Random Mutation Method

In order to avoid premature convergence and improve convergence precision, a hybrid random mutation method is proposed. The so-called hybrid random mutation is to conduct uniform random mutation in the early optimization stage of the algorithm and to conduct Gaussian random mutation in the final optimization stage. Particles can traverse on the optimal space with equal probability completely by conducting uniform random mutation, which can improve the diversity of solution and avoid premature convergence. Particles can move around the original position by conducting Gaussian mutation so as to avoid lowering optimization efficiency caused by the too large mutation range in the final optimization stage and improve convergence precision. The detail of hybrid random mutation is as follows:(12)xkt=xkt,rand>pmRandxmin,xmax,rand≤pm∩t≤tmax2Gaussianxkt,σ,rand≤pm∩t>tmax2.In (12), rand is a random number uniformly distributing on [0,1], [xmin,xmax] is optimization range, Rand(•) is the uniform random mutation operator, and Gaussian(•) is the Gaussian random mutation operator. pm is the mutation probability, which is defined as follows:(13)pm=0.2,t≤tmax4-415ttmax+415,t>tmax4.


### 3.3. Two-Stage Based Guider Particles Selection Method

In MOQPSO-D/S algorithm, the guider particles for each particle should be selected during iteration. The selection of guider particles from the external archive, with application of the roulette method [17], is not conducive to the distribution uniformity of solution, due to its strong randomicity. In order to address the issue above, a two-stage based guider particles selection method is proposed. The so-called two-stage based guider particles selection method includes the following two steps. Firstly, particles in the external archive should be sorted in descending order according to their crowding distance, and the 10% of particles with the biggest crowding distance is preliminary selected. Secondly, the particle with the closest sigma value to particle xk is selected as the guider particle gbestk1 for xk. The calculation method of the crowding distance for each particle in external archive can be referred to literature [19]. The calculation method of sigma value [20] for each particle is shown as in (14), and the detail of guider particle selection by adopting sigma value method is shown in Figure 1. When the particle position is updated based on QPSO with double-well, the second guider particle should be selected. The neighbor particle, with the larger distance to gbestk1, will be selected as the second guider particle gbestk2. The detail of this method for gbestk2 selection is shown in Figure 2.(14)sigmax=f12x-f22x,…,fo2x-fo2x∑l=1ofl2x.Figure 1
Sigma method based guider particle selection.Figure 2
Thegbestk2 selection.


### 3.4. Algorithm’s Complexity

Suppose that the number of objective functions in the multiobjective problem iso; the particle population is popsize; the size of external archive is Narchive. When operating one iteration, the complexity of MOQPSO-D/S for solving the multiobjective problem includes the following parts:(1)
The complexity of calculating the fitness of all particles equalsOo×popsize.(2)
The complexity of calculating the individual best positions of all particles equalsOo×popsize.(3)
After updating the positions of all particles, the complexity of comparing all particles in the population with the particles in the external archive equalsOo×popsize×Narchive.(4)
For updating the external archive, the complexity of conducting the crowding distance sorting equalsOo×popsize+Narchive2.(5)
After updating the external archive, the complexity of calculating the sigma values of all particles in the external archive equalsONarchive.Thus, when operating one iteration, the complexity of MOQPSO-D/S for solving the multiobjective problem equalsOo×popsize+Narchive2, which equals the complexity of MOQPSO-CD (the multiobjective quantum-behaved particle swarm optimization algorithm based on QPSO and crowding distance sorting) [17] and is acceptable.


## 3.1. Double/Single-Well Based Position Update Method

The so-called double/single-well based position update method includes two aspects. On one hand, in order to improve solutions’ diversity and avoid premature convergence, the particle position is updated based on QPSO with double-well algorithm in the early optimization stage. On the other hand, in order to improve convergence precision and convergence rate, the particle position is updated based on QPSO algorithm in the final optimization stage. The double/single-well based position update formula is as follows:(10)xkt+1=121+βktPk1t+1-βktPk2t±Lkt2ln1ukt,t≤tmax2Pk1t±Lkt2ln1ukt,t>tmax2.In (10), tmax is the maximum generation. If uk(t)>0.5, “±” would be chosen as “+”, or it would be chosen as “−”. The characteristic length Lk(t) is defined as in the following equation.(11)Lkt=α1+βktPk1t+1-βktPk2t-2xkt,t≤tmax22αPk1t-xkt,t>tmax2.


## 3.2. Hybrid Random Mutation Method

In order to avoid premature convergence and improve convergence precision, a hybrid random mutation method is proposed. The so-called hybrid random mutation is to conduct uniform random mutation in the early optimization stage of the algorithm and to conduct Gaussian random mutation in the final optimization stage. Particles can traverse on the optimal space with equal probability completely by conducting uniform random mutation, which can improve the diversity of solution and avoid premature convergence. Particles can move around the original position by conducting Gaussian mutation so as to avoid lowering optimization efficiency caused by the too large mutation range in the final optimization stage and improve convergence precision. The detail of hybrid random mutation is as follows:(12)xkt=xkt,rand>pmRandxmin,xmax,rand≤pm∩t≤tmax2Gaussianxkt,σ,rand≤pm∩t>tmax2.In (12), rand is a random number uniformly distributing on [0,1], [xmin,xmax] is optimization range, Rand(•) is the uniform random mutation operator, and Gaussian(•) is the Gaussian random mutation operator. pm is the mutation probability, which is defined as follows:(13)pm=0.2,t≤tmax4-415ttmax+415,t>tmax4.


## 3.3. Two-Stage Based Guider Particles Selection Method

In MOQPSO-D/S algorithm, the guider particles for each particle should be selected during iteration. The selection of guider particles from the external archive, with application of the roulette method [17], is not conducive to the distribution uniformity of solution, due to its strong randomicity. In order to address the issue above, a two-stage based guider particles selection method is proposed. The so-called two-stage based guider particles selection method includes the following two steps. Firstly, particles in the external archive should be sorted in descending order according to their crowding distance, and the 10% of particles with the biggest crowding distance is preliminary selected. Secondly, the particle with the closest sigma value to particle xk is selected as the guider particle gbestk1 for xk. The calculation method of the crowding distance for each particle in external archive can be referred to literature [19]. The calculation method of sigma value [20] for each particle is shown as in (14), and the detail of guider particle selection by adopting sigma value method is shown in Figure 1. When the particle position is updated based on QPSO with double-well, the second guider particle should be selected. The neighbor particle, with the larger distance to gbestk1, will be selected as the second guider particle gbestk2. The detail of this method for gbestk2 selection is shown in Figure 2.(14)sigmax=f12x-f22x,…,fo2x-fo2x∑l=1ofl2x.Figure 1
Sigma method based guider particle selection.Figure 2
Thegbestk2 selection.


## 3.4. Algorithm’s Complexity

Suppose that the number of objective functions in the multiobjective problem iso; the particle population is popsize; the size of external archive is Narchive. When operating one iteration, the complexity of MOQPSO-D/S for solving the multiobjective problem includes the following parts:(1)
The complexity of calculating the fitness of all particles equalsOo×popsize.(2)
The complexity of calculating the individual best positions of all particles equalsOo×popsize.(3)
After updating the positions of all particles, the complexity of comparing all particles in the population with the particles in the external archive equalsOo×popsize×Narchive.(4)
For updating the external archive, the complexity of conducting the crowding distance sorting equalsOo×popsize+Narchive2.(5)
After updating the external archive, the complexity of calculating the sigma values of all particles in the external archive equalsONarchive.Thus, when operating one iteration, the complexity of MOQPSO-D/S for solving the multiobjective problem equalsOo×popsize+Narchive2, which equals the complexity of MOQPSO-CD (the multiobjective quantum-behaved particle swarm optimization algorithm based on QPSO and crowding distance sorting) [17] and is acceptable.


## 4. Modeling WTA Problem Based on Fuzzy Multiobjective Programming

WTA in air and missile defense is to decide the problems of which fire unit should be chosen to intercept air targets, which air target will be intercepted by the chosen fire unit, and how many interceptors should be launched to intercept the chosen air target by the chosen fire unit. The goal of WTA decision-making is to obtain maximum interception efficiency and minimum interception consumption.


### 4.1. Model Assumptions and Parameter Settings

Suppose that in air and missile defense operation,m fire units W={W1,W2,…,Wm} need intercept n incoming targets T={T1,T2,…,Tn}, where T1~Ta belong to Type A (aerodynamic target), the remaining targets belong to Type B (ballistic target). Assuming that NA or more interceptors are needed to intercept a target belonging to Type A and NB or more interceptors are needed to intercept a target belonging to Type B, different types of targets need to be intercepted by different interceptors; for example, a target belonging to Type A can only be intercepted by Type A interceptors. The fire unit Wi   (i=1,2,…,m) can intercept different types of targets but can only intercept a type of targets simultaneously. Due to the sensor error and other issues, the targets threat value is uncertain. The threat value of Tj  (j=1,2,…,n) will be characterized with application of triangular fuzzy number w~j=(wjp,wjm,wjo), where wjp,wjm, and wjo denote the most pessimistic, the most likely, and the most optimistic threat value, respectively. Wi   (i=1,2,…,m) launches yij interceptors to intercept Tj  (j=1,2,…,n), and the single-shot kill probability is Dij.The correlative parameters are defined as follows:miA denotes the number of Type A targets that Wi can intercept simultaneously,m^iA denotes the number of Type A interceptors that Wi can launch simultaneously,NiA denotes the number of Type A interceptors that Wi stores,CiA denotes the value of a Type A interceptor that Wi stores,miB denotes the number of Type B targets that Wi can intercept simultaneously,m^iB denotes the number of Type B interceptors that Wi can launch simultaneously,NiB denotes the number of Type B interceptors that Wi stores,CiB denotes the value of a Type B interceptor that Wi stores.


### 4.2. Modeling Based on Fuzzy Multiobjective Programming

Taking the maximizing interception efficiency and the minimizing interception consumption as goals, a multiobjective optimization based WTA model is established as follows.(15)maxf1=∑j=1nw~j1-∏i=1m1-Dijyij,(16)minf2=∑i=1m∑j=1aCiAyij+∑j=a+1nCiByij.The model needs to meet the resource constraints and the firepower constraint. For the fire unitWi   (i=1,2,…,m), the resource constraints are needed to be met as follows.(1) Interception Constraint. The fire unit Wi can only intercept a type of target simultaneously; thus the model should meet(17)δ∑j=1ayij+δ∑j=a+1nyij≤1.In (17), δ(ξ)=1,ξ>00,ξ≤0.(2) Interception Capability Constraint. The fire unit Wi can intercept less than miA Type A targets or miB Type B targets simultaneously; thus the model should meet(18)∑j=1aδyij≤miA∑j=a+1nδyij≤miB.(3) Number Constraint of Interceptors. The fire unit Wi can launch less than m^iA and NiA to intercept Type A targets or launch less than m^iB and NiB to intercept Type B targets; thus the model should meet(19)∑j=1ayij≤minm^iA,NiA∑j=a+1nyij≤minm^iB,NiB.For intercepting the incoming targetTj  (j=1,2,…,n), the firepower constraint is needed to be met as follows.N
A or more Type A interceptors are needed to intercept a Type A target, and NB or more Type B interceptors are needed to intercept a Type B target; thus the model should meet(20)∑i=1myij≥NA,j=1,2,…,a∑i=1myij≥NB,j=a+1,a+2,…,n.


### 4.3. Clarifying Model Equivalently

This WTA model cannot be resolved by adopting multiobjective optimization algorithm, for its objective function containing fuzzy parameter shown as (15). Therefore, it is necessary to clarify the fuzzy multiobjective programming based WTA model to be a certain model equivalently. The necessity degree principle of uncertainty theory [21] is adopted to turn the objective function with fuzzy parameter into a certain objective function. The detailed process is as follows.Firstly, the objective function shown as (15) is turned into a fuzzy chance constrained programming based model(21)mincs.t.Nec∑j=1nw~j∏i=1m1-Dijyij≤c≥θ.In (21), θ denotes the necessity degree that the whole rest incoming targets’ threat value is equal to or less than c.Secondly, the fuzzy chance constrained programming based model shown as in (21) is turned into certain form with application of the necessity degree principle of uncertainty theory. Because Nec∑j=1nw~j∏i=1m1-Dijyij≤c≥θ can be turned into ∑j=1nwjm+θwjo-wjm∏i=1m1-Dijyij≤c equivalently, (17) can be turned into(22)mincs.t.∑j=1nwjm+θwjo-wjm∏i=1m1-Dijyij≤c.Finally, the certain objective function is obtained as follows:(23)minf1′=∑j=1nwjm+θwjo-wjm∏i=1m1-Dijyij.In this way, the fuzzy multiobjective programming based WTA model is turned into a certain multiobjective optimization model. The model clarified could be solved with application of multiobjective optimization algorithm. The WTA model constructed in this paper will be solved by adopting MOQPSO-D/S algorithm in the next section.


## 4.1. Model Assumptions and Parameter Settings

Suppose that in air and missile defense operation,m fire units W={W1,W2,…,Wm} need intercept n incoming targets T={T1,T2,…,Tn}, where T1~Ta belong to Type A (aerodynamic target), the remaining targets belong to Type B (ballistic target). Assuming that NA or more interceptors are needed to intercept a target belonging to Type A and NB or more interceptors are needed to intercept a target belonging to Type B, different types of targets need to be intercepted by different interceptors; for example, a target belonging to Type A can only be intercepted by Type A interceptors. The fire unit Wi   (i=1,2,…,m) can intercept different types of targets but can only intercept a type of targets simultaneously. Due to the sensor error and other issues, the targets threat value is uncertain. The threat value of Tj  (j=1,2,…,n) will be characterized with application of triangular fuzzy number w~j=(wjp,wjm,wjo), where wjp,wjm, and wjo denote the most pessimistic, the most likely, and the most optimistic threat value, respectively. Wi   (i=1,2,…,m) launches yij interceptors to intercept Tj  (j=1,2,…,n), and the single-shot kill probability is Dij.The correlative parameters are defined as follows:miA denotes the number of Type A targets that Wi can intercept simultaneously,m^iA denotes the number of Type A interceptors that Wi can launch simultaneously,NiA denotes the number of Type A interceptors that Wi stores,CiA denotes the value of a Type A interceptor that Wi stores,miB denotes the number of Type B targets that Wi can intercept simultaneously,m^iB denotes the number of Type B interceptors that Wi can launch simultaneously,NiB denotes the number of Type B interceptors that Wi stores,CiB denotes the value of a Type B interceptor that Wi stores.


## 4.2. Modeling Based on Fuzzy Multiobjective Programming

Taking the maximizing interception efficiency and the minimizing interception consumption as goals, a multiobjective optimization based WTA model is established as follows.(15)maxf1=∑j=1nw~j1-∏i=1m1-Dijyij,(16)minf2=∑i=1m∑j=1aCiAyij+∑j=a+1nCiByij.The model needs to meet the resource constraints and the firepower constraint. For the fire unitWi   (i=1,2,…,m), the resource constraints are needed to be met as follows.(1) Interception Constraint. The fire unit Wi can only intercept a type of target simultaneously; thus the model should meet(17)δ∑j=1ayij+δ∑j=a+1nyij≤1.In (17), δ(ξ)=1,ξ>00,ξ≤0.(2) Interception Capability Constraint. The fire unit Wi can intercept less than miA Type A targets or miB Type B targets simultaneously; thus the model should meet(18)∑j=1aδyij≤miA∑j=a+1nδyij≤miB.(3) Number Constraint of Interceptors. The fire unit Wi can launch less than m^iA and NiA to intercept Type A targets or launch less than m^iB and NiB to intercept Type B targets; thus the model should meet(19)∑j=1ayij≤minm^iA,NiA∑j=a+1nyij≤minm^iB,NiB.For intercepting the incoming targetTj  (j=1,2,…,n), the firepower constraint is needed to be met as follows.N
A or more Type A interceptors are needed to intercept a Type A target, and NB or more Type B interceptors are needed to intercept a Type B target; thus the model should meet(20)∑i=1myij≥NA,j=1,2,…,a∑i=1myij≥NB,j=a+1,a+2,…,n.


## 4.3. Clarifying Model Equivalently

This WTA model cannot be resolved by adopting multiobjective optimization algorithm, for its objective function containing fuzzy parameter shown as (15). Therefore, it is necessary to clarify the fuzzy multiobjective programming based WTA model to be a certain model equivalently. The necessity degree principle of uncertainty theory [21] is adopted to turn the objective function with fuzzy parameter into a certain objective function. The detailed process is as follows.Firstly, the objective function shown as (15) is turned into a fuzzy chance constrained programming based model(21)mincs.t.Nec∑j=1nw~j∏i=1m1-Dijyij≤c≥θ.In (21), θ denotes the necessity degree that the whole rest incoming targets’ threat value is equal to or less than c.Secondly, the fuzzy chance constrained programming based model shown as in (21) is turned into certain form with application of the necessity degree principle of uncertainty theory. Because Nec∑j=1nw~j∏i=1m1-Dijyij≤c≥θ can be turned into ∑j=1nwjm+θwjo-wjm∏i=1m1-Dijyij≤c equivalently, (17) can be turned into(22)mincs.t.∑j=1nwjm+θwjo-wjm∏i=1m1-Dijyij≤c.Finally, the certain objective function is obtained as follows:(23)minf1′=∑j=1nwjm+θwjo-wjm∏i=1m1-Dijyij.In this way, the fuzzy multiobjective programming based WTA model is turned into a certain multiobjective optimization model. The model clarified could be solved with application of multiobjective optimization algorithm. The WTA model constructed in this paper will be solved by adopting MOQPSO-D/S algorithm in the next section.


## 5. Solving WTA Model Based on MOQPSO-D/S

### 5.1. Particle Coding

The MOQPSO-D/S algorithm cannot solve the discrete problem directly. Therefore, it is necessary to code the position of particles when the WTA model is resolved by adopting MOQPSO-D/S algorithm. The position of particles will be coded with application of matrix encoding.y denotes the position of a particle, and it can be coded as follows:(24)y=y11y12⋯y1ny21y22⋯y2n⋮⋮⋱⋮ym1ym2⋯ymn.In (24), yij denotes the number of interceptors that Wi launches to intercept Tj.


### 5.2. Illegal Particle Code Adjusting

When initializing and updating the particle position, the particle position may not satisfy the constraints of the WTA model; namely, the infeasible solution comes out. Therefore, it is necessary to check the particle coding legitimacy and to adjust the illegal particle code. The process of particle code checking and illegal code adjusting is shown in Figure3.Figure 3
The process of particle code checking and illegal code adjusting.For example, when a particle code does not meet the interception capability constraint for fire unitWi   (i=1,2,…,m), Operation 2 should be conducted as follows:Firstly, judge whether∑j=1ayij=0 or ∑j=a+1nyij=0.Case 1.
One has(25)∑j=1ayij=0
C1S1 (Case1 Step 1). Let count_A=∑j=1aδyij; if δyij=1, then store the subscript j into the set ΩA.
C1S2. Let Nexc_A=count_A-miA.
C1S3. Select randomly Nexc_A subscript from the set ΩA, and let the corresponding yij=0.Case 2.
One has(26)∑j=a+1nyij=0
C2S2. Let count_B=∑j=a+1nδyij; if δyij=1, then store the subscript j into the set ΩB.
C1S2. Let Nexc_B=count_B-miB.
C1S3. Select randomly Nexc_B subscript from the set ΩB, and let the corresponding yij=0.Finally, output the particle code that meets the interception capability constraint for fire unitWi  (i=1,2,…,m).Due to limitations on space, Operation 1, Operation 3, and Operation 4 will not be introduced here.


### 5.3. Solving Steps

After the particles are encoded and the illegal particles are adjusted, the MOQPSO-D/S algorithm can be adopted to solve the fuzzy multiobjective programming based WTA model clarified equivalently. The detailed steps are as follows.Step 1.
Set the algorithm parameters, and initialize the particle swarm randomly.Step 2.
Check out whether every particle code meets the constraints. If a particle code does not meet the constraints, then adjust it. Otherwise, turn to Step3.Step 3.
Calculate the fitness of every particlef(y)=(f1′(y),f2(y)) according to (16) and (23); let pbest=y,f(pbest)=f(y), where pbest denotes the individual best position of y.Step 4.
Initialize the external storage file ARCHIVE. Calculate the dominated relationship between particles according to their fitness, and store the nondominated particles into ARCHIVE, then calculate the crowding distance and sigma value of particles in ARCHIVE.Step 5.
Initialize iteration, and lett=1.Step 6.
Select guider particlesgbest1 and gbest2 with application of two-stage method. If t>tmax/2, the guider particle gbest2 would not be selected. Update the position of every particle according to (10).Step 7.
Judge whether the mutation condition is satisfied or not. If it is satisfied, conduct mutation for the particle by adopting the hybrid random mutation method. Otherwise, turn to Step8.Step 8.
Check out whether every particle code meets the constraints. If a particle code does not meet the constraints, then adjust it. Otherwise, turn to Step9.Step 9.
Calculate the fitness of every particlef(y(t))=(f1′(y(t)),f2(y(t))), and calculate the dominated relationship between the current position and the individual best position of every particle. If the current position dominates the individual best position of the particle yk, then let pbestk=yk(t),f(pbestk)=f(yk(t)). If the current position and the individual best position do not dominate each other, then conduct the operation pbestk=yk(t), f(pbestk)=f(yk(t)) with the probability of 0.5. Otherwise, turn to Step 10.Step 10.
Calculate the dominated relationship between every particle and particles in ARCHIVE. If the particleyk dominates the particle yg in ARCHIVE, then store the particle yk into ARCHIVE, delete the particle yg, and turn to Step 11. If the particle yk and the particle yg do not dominate each other, then store the particle yk into ARCHIVE, and turn to Step 11. Otherwise, turn to Step 12.Step 11.
Calculate the crowding distance of particles in ARCHIVE updated, and calculate the sigma value of the particles newly stored. Judge whether the number of particles in ARCHIVE exceeds its capacity. If it exceeds the capacity, then delete the particle with the smallest crowning distance. Otherwise, turn to Step12.Step 12.
Lett=t+1. If t>tmax, then end the iteration, and output all the particles in ARCHIVE as the WTA alternatives for air and missile defense combat. Otherwise, turn to Step 6.


## 5.1. Particle Coding

The MOQPSO-D/S algorithm cannot solve the discrete problem directly. Therefore, it is necessary to code the position of particles when the WTA model is resolved by adopting MOQPSO-D/S algorithm. The position of particles will be coded with application of matrix encoding.y denotes the position of a particle, and it can be coded as follows:(24)y=y11y12⋯y1ny21y22⋯y2n⋮⋮⋱⋮ym1ym2⋯ymn.In (24), yij denotes the number of interceptors that Wi launches to intercept Tj.


## 5.2. Illegal Particle Code Adjusting

When initializing and updating the particle position, the particle position may not satisfy the constraints of the WTA model; namely, the infeasible solution comes out. Therefore, it is necessary to check the particle coding legitimacy and to adjust the illegal particle code. The process of particle code checking and illegal code adjusting is shown in Figure3.Figure 3
The process of particle code checking and illegal code adjusting.For example, when a particle code does not meet the interception capability constraint for fire unitWi   (i=1,2,…,m), Operation 2 should be conducted as follows:Firstly, judge whether∑j=1ayij=0 or ∑j=a+1nyij=0.Case 1.
One has(25)∑j=1ayij=0
C1S1 (Case1 Step 1). Let count_A=∑j=1aδyij; if δyij=1, then store the subscript j into the set ΩA.
C1S2. Let Nexc_A=count_A-miA.
C1S3. Select randomly Nexc_A subscript from the set ΩA, and let the corresponding yij=0.Case 2.
One has(26)∑j=a+1nyij=0
C2S2. Let count_B=∑j=a+1nδyij; if δyij=1, then store the subscript j into the set ΩB.
C1S2. Let Nexc_B=count_B-miB.
C1S3. Select randomly Nexc_B subscript from the set ΩB, and let the corresponding yij=0.Finally, output the particle code that meets the interception capability constraint for fire unitWi  (i=1,2,…,m).Due to limitations on space, Operation 1, Operation 3, and Operation 4 will not be introduced here.


## 5.3. Solving Steps

After the particles are encoded and the illegal particles are adjusted, the MOQPSO-D/S algorithm can be adopted to solve the fuzzy multiobjective programming based WTA model clarified equivalently. The detailed steps are as follows.Step 1.
Set the algorithm parameters, and initialize the particle swarm randomly.Step 2.
Check out whether every particle code meets the constraints. If a particle code does not meet the constraints, then adjust it. Otherwise, turn to Step3.Step 3.
Calculate the fitness of every particlef(y)=(f1′(y),f2(y)) according to (16) and (23); let pbest=y,f(pbest)=f(y), where pbest denotes the individual best position of y.Step 4.
Initialize the external storage file ARCHIVE. Calculate the dominated relationship between particles according to their fitness, and store the nondominated particles into ARCHIVE, then calculate the crowding distance and sigma value of particles in ARCHIVE.Step 5.
Initialize iteration, and lett=1.Step 6.
Select guider particlesgbest1 and gbest2 with application of two-stage method. If t>tmax/2, the guider particle gbest2 would not be selected. Update the position of every particle according to (10).Step 7.
Judge whether the mutation condition is satisfied or not. If it is satisfied, conduct mutation for the particle by adopting the hybrid random mutation method. Otherwise, turn to Step8.Step 8.
Check out whether every particle code meets the constraints. If a particle code does not meet the constraints, then adjust it. Otherwise, turn to Step9.Step 9.
Calculate the fitness of every particlef(y(t))=(f1′(y(t)),f2(y(t))), and calculate the dominated relationship between the current position and the individual best position of every particle. If the current position dominates the individual best position of the particle yk, then let pbestk=yk(t),f(pbestk)=f(yk(t)). If the current position and the individual best position do not dominate each other, then conduct the operation pbestk=yk(t), f(pbestk)=f(yk(t)) with the probability of 0.5. Otherwise, turn to Step 10.Step 10.
Calculate the dominated relationship between every particle and particles in ARCHIVE. If the particleyk dominates the particle yg in ARCHIVE, then store the particle yk into ARCHIVE, delete the particle yg, and turn to Step 11. If the particle yk and the particle yg do not dominate each other, then store the particle yk into ARCHIVE, and turn to Step 11. Otherwise, turn to Step 12.Step 11.
Calculate the crowding distance of particles in ARCHIVE updated, and calculate the sigma value of the particles newly stored. Judge whether the number of particles in ARCHIVE exceeds its capacity. If it exceeds the capacity, then delete the particle with the smallest crowning distance. Otherwise, turn to Step12.Step 12.
Lett=t+1. If t>tmax, then end the iteration, and output all the particles in ARCHIVE as the WTA alternatives for air and missile defense combat. Otherwise, turn to Step 6.


## 6. Simulation Experiment and Analysis

In order to verify the performance of MOQPSO-D/S algorithm for solving multiobjective optimization problems and WTA problem in air and missile defense operation, two simulation experiments are designed. In experiment 1, the convergence and solution distribution uniformity of the MOQPSO-D/S algorithm for solving multiobjective optimization problem are validated. In experiment 2, after the simulation scenarios are set, the rationality of the fuzzy multiobjective programming based WTA model and the feasibility of the MOQPSO-D/S algorithm for solving the model constructed in this paper are validated by solving the WTA model with application of MOQPSO-D/S algorithm coded. The experiments above are conducted on the computer with INTEL CORE i5-4590, 3.3 GHz CUP, and 4 G RAM after programming with MATLAB 2013a. The capability of ARCHIVE for all multiobjective optimization algorithms is 100.


### 6.1. Verify the Performance of MOQPSO-D/S for Solving Multiobjective Optimization Problems

The test functions ZDT1~ZDT4 [22] are chosen to verify the performance of MOQPSO-D/S algorithm for solving multiobjective optimization problems. The multiobjective quantum-behaved particle swarm optimization algorithm with double-potential well and share-learning (MOQPSO-DPS) [14], the MOQPSO-CD [17], and the multiobjective particle swarm optimization algorithm based on crowding distance sorting (MOPSO-CD) [23] are chosen as the comparative algorithms. Set the particle population to be 100, and set the iteration number to be 200. The variables number of the test functions ZDT1~ZDT3 is 30, and the variables number of the test function ZDT4 is 10. The algorithms run 30 times independently. Select the solution of one time from the 30 times of every algorithm for solving test functions ZDT1~ZDT4 randomly, and the results are shown in Figure 4.Figure 4
Pareto front obtained by solving ZDT1~ZDT4 with application of four algorithms.
(a)
ZDT1
(b)
ZDT2
(c)
ZDT3
(d)
ZDT4Error ration (ER) [14] and spacing metric (SP) [24] are chosen to weigh the convergence of four algorithms for solving the test functions and the distribution uniformity of the Pareto optimal solution, respectively. And, the mean running time is chosen to weigh the operating rate. The calculations of ER and SP are shown in (27) and (28). ER index represents the ratio of the Pareto solutions obtained by four algorithms in the true Pareto optimal solution set. The smaller the ER value is, the better the convergence is. The smaller the SP value is, the more uniformly the Pareto optimal solutions distribute.(27)ER=∑i=1nPeinP.In (27), ei denotes the probability that a Pareto optimal solution is not in the true Pareto optimal solution set, if it is in, ei=0; or, ei=1. nP denotes the number of elements in a Pareto optimal solution set.(28)SP=1nP-1∑k=1nPd--dk2.In (28), dk denotes the Hamming distance between the dot corresponding to solution xk and its neighbor dot in the Pareto front. The calculation of dk is shown as in (29). d- denotes the mean distance, shown in (30).(29)dk=minr≠k∑ξ=1ofξxk-fξxr,(30)d-=1nP∑k=1nPdk.The comparative results are shown in Tables 1–3, where the ER values, SP values, and the mean running time values are the average values of four algorithms for solving the test functions 30 times.Table 1
ER of four algorithms for solving the test functions.
Algorithms
Test functions
ZDT1
ZDT2
ZDT3
ZDT4
MOQPSO-D/S
0.0097
0.0087
0.0083
0.0290
MOQPSO-DPS
0.0150
0.0127
0.0187
0.0370
MOQPSO-CD
0.0113
0.0253
0.0343
0.0337
MOPSO-CD
0.0283
0.0363
0.1143
0.2783Table 2
SP of four algorithms for solving the test functions.
Algorithms
Test functions
ZDT1
ZDT2
ZDT3
ZDT4
MOQPSO-D/S
0.0043
0.0075
0.0089
0.0074
MOQPSO-DPS
0.0117
0.0073
0.0114
0.0131
MOQPSO-CD
0.0071
0.0101
0.0121
0.0087
MOPSO-CD
0.0079
0.0084
0.0143
0.0241Table 3
Mean running time of four algorithms for solving the test functions.
Algorithms
Test functions
ZDT1
ZDT2
ZDT3
ZDT4
MOQPSO-D/S
5.1134
6.7634
4.0197
4.5163
MOQPSO-DPS
7.2471
7.3521
4.4110
6.0217
MOQPSO-CD
8.3283
8.1126
4.7853
6.4313
MOPSO-CD
10.3590
11.8781
6.7596
8.8357From Table1, we can see that the ER indexes of MOQPSO-D/S for solving the test functions are smaller than the remaining algorithms’. Combining with Figure 4, it can be concluded that the convergence of MOQPSO-D/S is better than the remaining algorithms’. And, it is proved that the double/single-well based position update method and the hybrid random mutation method can improve convergence precision.From Table2, we can see that the SP indexes of MOQPSO-D/S for solving the test functions are smaller than the remaining algorithms’, except for solving ZDT2. Combining with Figure 4, it can be concluded that the distribution uniformity of the Pareto optimal solution of MOQPSO-D/S is better than the remaining algorithms’. And, it is proved that the two-stage based guider particles selection method can improve solution distribution uniformity.From Table3, we can see that the mean running time of MOQPSO-D/S for solving the test functions is smaller than the remaining algorithms’. It can be concluded that the operating rate of MOQPSO-D/S is better than the remaining algorithms’. And, it is proved that the two-stage based guider particles selection method can improve operating rate.


### 6.2. Verify the Performance of MOQPSO-D/S for Solving Uncertain WTA Problem

In order to validate the rationality of the fuzzy multiobjective programming based WTA model and to verify the performance of the MOQPSO-D/S algorithm for solving the model simultaneously, the following case for air and missile defense combat is given.Suppose that, in air and missile defense operation, 6 fire unitsWi    (i=1,2,…,6) need intercept 5 incoming targets Tj   (j=1,2,…,5). The number of targets that Wi    (i=1,2,…,6) can intercept simultaneously, the number of interceptors that Wi can launch simultaneously, the number of interceptors that Wi stores, and the dimensionless value of an interceptor that Wi stores are listed in Table 4. The incoming targets T1~T3 belong to Type A, and the remaining targets belong to Type B. The target threat values are weighed with application of triangular fuzzy number as follows: w~1=(0.41,0.43,0.45), w~2=(0.48,0.51,0.54), w~3=(0.33,0.35,0.37), w~4=(0.87,0.89,0.91), w~5=(0.90,0.92,0.94). The least interceptors to intercept a Type A target and a Type B target are NA=1 and NB=2, respectively. The necessity degree θ is set to be 0.8. The single-shot kill probability of the fire unit Wi   (i=1,2,…,6) that launches an interceptor to intercept the target Tj   (j=1,2,…,5) is listed in Table 5. The particle encoding method and the illegal particle code adjusting method are adopted in the subsequent algorithms for solving the example.Table 4
Parameters value.
W
1
W
2
W
3
W
4
W
5
W
6
m
i
a
1
3
2
5
4
0
m
^
i
a
3
6
4
10
8
0
m
i
b
0
0
0
1
1
1
m
^
i
b
0
0
0
2
2
2
N
i
a
6
6
8
9
8
0
N
i
b
0
0
0
2
2
2
C
i
a
0.1
0.2
0.15
0.25
0.3
-
C
i
b
-
-
-
0.45
0.5
0.6Table 5
Single-shot kill probability.
W
1
W
2
W
3
W
4
W
5
W
6
T
1
0.45
0.32
0.54
0.47
0.75
0
T
2
0.5
0.27
0.34
0.73
0.78
0
T
3
0.28
0.62
0.58
0.56
0.57
0
T
4
0
0
0
0.63
0.55
0.7
T
5
0
0
0
0.45
0.5
0.85


#### 6.2.1. Analyze the Result of MOQPSO-D/S for Solving Uncertain WTA Problem

With different population size and different iterations, the WTA case is solved with application of MOQPSO-D/S algorithm. The algorithm for solving the WTA case under every condition runs 30 times independently. The mean running time is listed in Table6. It can be seen from Table 6 that the real-time performance of MOQPSO-D/S algorithm can meet the need of WTA in air and missile defense combat.Table 6
Mean running time/s.
Population size
Iterations
50
100
200
50
0.7504
1.4830
2.9733
100
1.5179
2.9313
5.9883
200
3.0576
6.0737
11.8876With the population size and iterations both being 50, the solution of the WTA case solved by MOQPSO-D/S algorithm is shown in Figure5. It can be seen from Figure 5 that the Pareto optimal solution, obtained by MOQPSO-D/S algorithm for solving the fuzzy multiobjective programming based WTA model in air and missile defense combat, distributes uniformly, where every dot refers to a WTA alternative. Therefore, it can provide a series of alternatives for commanders to make decision. The commanders could make the final decision by choosing a suitable WTA alternative from the solution with consideration of the battlefield situation and risk preference. Comparing with the single objective optimization based WTA method, which can only give a WTA alternative, this method could preferably adapt to the change of the battlefield situation and combine the tactical thought of the commander to make decision.Figure 5
Solution of MOQPSO-D/S for solving the case.


#### 6.2.2. Compare with the Single Objective Optimization Based WTA Method

In order to validate the feasibility of the solution obtained by MOQPSO-D/S algorithm for solving WTA model, in comparison, the example is solved by adopting the method proposed in literature [6], taking maximizing effectiveness-cost ratio as the objective function and adopting IDPSO algorithm to solve the WTA model. With the population size and iterations being both 50, the single objective optimization method for solving the WTA case runs 30 times. The mean running time is 0.7193 s. The best fitness over generations based on IDPSO is shown as in Figure 6. The optimal WTA result is (31)T1T2T3T4T5W1W2W3W4W5W6001000200000001000000200000011.Figure 6
Best fitness over generations based on IDPSO algorithm.The optimal result is just the WTA alternative corresponding to the dot marked by a red five-angle star shown in Figure5. And the mean running time is a little shorter than the MOQPSO-D/S algorithm with the same population size and iterations. It is proved that MOQPSO-D/S algorithm, when solving the WTA model, would disperse computation resource but can still obtain the satisfactory solution containing the optimal result obtained by single optimization method. Meanwhile, the fuzzy multiobjective programming WTA model for air and missile defense is proved to be reasonable and adopting MOQPSO-D/S algorithm to solve the model is feasible.


#### 6.2.3. Compare with Other Multiobjective Optimization Algorithm for Solving WTA Method

In order to further verify the performance of MOQPSO-D/S algorithm for solving WTA problem, in comparison, the example is solved by adopting MOQPSO-DPS, MOQPSO-CD, and MOPSO-CD algorithm with the population size and iterations being both 100. Every algorithm for solving the WTA case runs 30 times independently. The statistical results for running time of every algorithm are shown in Figure7. It can be seen from Figure 7 that the real-time performance of MOQPSO-D/S algorithm is obviously superior to the comparable algorithms’.Figure 7
Statistical results for running time of four algorithms.The convergence of every algorithm for solving the WTA case is weighed by dominated ratioC(X1,X2) [25]. The calculation of C(X1,X2) is shown as in (32). X1 and X2 denote the Pareto optimal solution sets obtained by two algorithms, respectively. If C(X1,X2) is bigger than C(X2,X1), the convergence of the algorithm corresponding to X1 is superior to the algorithm corresponding to X2. The convergence comparison among four algorithms is listed in Table 7. A1,A2,A3, and A4 listed in Table 7 denote the Pareto optimal solution sets obtained by MOQPSO-D/S, MOQPSO-DPS, MOQPSO-CD, and MOPSO-CD algorithm, respectively. It can be concluded from Table 7 that the convergence of MOQPSO-D/S algorithm for solving the WTA case is the best of four algorithms.(32)CX1,X2=numelx2∈X2∣∃x1∈X1:x1≺x2numelX2.In (32), the function numel(•) is used to count the number of elements in a set.Table 7
Dominated ratio among four algorithms’ Pareto optimal sets.
X
2
X
1
A1
A2
A3
A4
A1
0
0
0.0185
0.0370
A2
0.0588
0
0.1176
0.1765
A3
0.0556
0
0
0
A4
0.0755
0.0755
0.0755
0The distribution uniformity of the Pareto optimal solution set obtained by every algorithm for solving the WTA case is weighed by SP. The statistical results for SP value of the Pareto optimal solution set obtained by every algorithm are shown in Figure8. It can be concluded from Figure 8 that the Pareto optimal solutions obtained by MOQPSO-D/S algorithm for solving the WTA case distribute the most uniformly, and its distribution uniformity is the most stable.Figure 8
Statistical results for SP value of four algorithms.


## 6.1. Verify the Performance of MOQPSO-D/S for Solving Multiobjective Optimization Problems

The test functions ZDT1~ZDT4 [22] are chosen to verify the performance of MOQPSO-D/S algorithm for solving multiobjective optimization problems. The multiobjective quantum-behaved particle swarm optimization algorithm with double-potential well and share-learning (MOQPSO-DPS) [14], the MOQPSO-CD [17], and the multiobjective particle swarm optimization algorithm based on crowding distance sorting (MOPSO-CD) [23] are chosen as the comparative algorithms. Set the particle population to be 100, and set the iteration number to be 200. The variables number of the test functions ZDT1~ZDT3 is 30, and the variables number of the test function ZDT4 is 10. The algorithms run 30 times independently. Select the solution of one time from the 30 times of every algorithm for solving test functions ZDT1~ZDT4 randomly, and the results are shown in Figure 4.Figure 4
Pareto front obtained by solving ZDT1~ZDT4 with application of four algorithms.
(a)
ZDT1
(b)
ZDT2
(c)
ZDT3
(d)
ZDT4Error ration (ER) [14] and spacing metric (SP) [24] are chosen to weigh the convergence of four algorithms for solving the test functions and the distribution uniformity of the Pareto optimal solution, respectively. And, the mean running time is chosen to weigh the operating rate. The calculations of ER and SP are shown in (27) and (28). ER index represents the ratio of the Pareto solutions obtained by four algorithms in the true Pareto optimal solution set. The smaller the ER value is, the better the convergence is. The smaller the SP value is, the more uniformly the Pareto optimal solutions distribute.(27)ER=∑i=1nPeinP.In (27), ei denotes the probability that a Pareto optimal solution is not in the true Pareto optimal solution set, if it is in, ei=0; or, ei=1. nP denotes the number of elements in a Pareto optimal solution set.(28)SP=1nP-1∑k=1nPd--dk2.In (28), dk denotes the Hamming distance between the dot corresponding to solution xk and its neighbor dot in the Pareto front. The calculation of dk is shown as in (29). d- denotes the mean distance, shown in (30).(29)dk=minr≠k∑ξ=1ofξxk-fξxr,(30)d-=1nP∑k=1nPdk.The comparative results are shown in Tables 1–3, where the ER values, SP values, and the mean running time values are the average values of four algorithms for solving the test functions 30 times.Table 1
ER of four algorithms for solving the test functions.
Algorithms
Test functions
ZDT1
ZDT2
ZDT3
ZDT4
MOQPSO-D/S
0.0097
0.0087
0.0083
0.0290
MOQPSO-DPS
0.0150
0.0127
0.0187
0.0370
MOQPSO-CD
0.0113
0.0253
0.0343
0.0337
MOPSO-CD
0.0283
0.0363
0.1143
0.2783Table 2
SP of four algorithms for solving the test functions.
Algorithms
Test functions
ZDT1
ZDT2
ZDT3
ZDT4
MOQPSO-D/S
0.0043
0.0075
0.0089
0.0074
MOQPSO-DPS
0.0117
0.0073
0.0114
0.0131
MOQPSO-CD
0.0071
0.0101
0.0121
0.0087
MOPSO-CD
0.0079
0.0084
0.0143
0.0241Table 3
Mean running time of four algorithms for solving the test functions.
Algorithms
Test functions
ZDT1
ZDT2
ZDT3
ZDT4
MOQPSO-D/S
5.1134
6.7634
4.0197
4.5163
MOQPSO-DPS
7.2471
7.3521
4.4110
6.0217
MOQPSO-CD
8.3283
8.1126
4.7853
6.4313
MOPSO-CD
10.3590
11.8781
6.7596
8.8357From Table1, we can see that the ER indexes of MOQPSO-D/S for solving the test functions are smaller than the remaining algorithms’. Combining with Figure 4, it can be concluded that the convergence of MOQPSO-D/S is better than the remaining algorithms’. And, it is proved that the double/single-well based position update method and the hybrid random mutation method can improve convergence precision.From Table2, we can see that the SP indexes of MOQPSO-D/S for solving the test functions are smaller than the remaining algorithms’, except for solving ZDT2. Combining with Figure 4, it can be concluded that the distribution uniformity of the Pareto optimal solution of MOQPSO-D/S is better than the remaining algorithms’. And, it is proved that the two-stage based guider particles selection method can improve solution distribution uniformity.From Table3, we can see that the mean running time of MOQPSO-D/S for solving the test functions is smaller than the remaining algorithms’. It can be concluded that the operating rate of MOQPSO-D/S is better than the remaining algorithms’. And, it is proved that the two-stage based guider particles selection method can improve operating rate.


## 6.2. Verify the Performance of MOQPSO-D/S for Solving Uncertain WTA Problem

In order to validate the rationality of the fuzzy multiobjective programming based WTA model and to verify the performance of the MOQPSO-D/S algorithm for solving the model simultaneously, the following case for air and missile defense combat is given.Suppose that, in air and missile defense operation, 6 fire unitsWi    (i=1,2,…,6) need intercept 5 incoming targets Tj   (j=1,2,…,5). The number of targets that Wi    (i=1,2,…,6) can intercept simultaneously, the number of interceptors that Wi can launch simultaneously, the number of interceptors that Wi stores, and the dimensionless value of an interceptor that Wi stores are listed in Table 4. The incoming targets T1~T3 belong to Type A, and the remaining targets belong to Type B. The target threat values are weighed with application of triangular fuzzy number as follows: w~1=(0.41,0.43,0.45), w~2=(0.48,0.51,0.54), w~3=(0.33,0.35,0.37), w~4=(0.87,0.89,0.91), w~5=(0.90,0.92,0.94). The least interceptors to intercept a Type A target and a Type B target are NA=1 and NB=2, respectively. The necessity degree θ is set to be 0.8. The single-shot kill probability of the fire unit Wi   (i=1,2,…,6) that launches an interceptor to intercept the target Tj   (j=1,2,…,5) is listed in Table 5. The particle encoding method and the illegal particle code adjusting method are adopted in the subsequent algorithms for solving the example.Table 4
Parameters value.
W
1
W
2
W
3
W
4
W
5
W
6
m
i
a
1
3
2
5
4
0
m
^
i
a
3
6
4
10
8
0
m
i
b
0
0
0
1
1
1
m
^
i
b
0
0
0
2
2
2
N
i
a
6
6
8
9
8
0
N
i
b
0
0
0
2
2
2
C
i
a
0.1
0.2
0.15
0.25
0.3
-
C
i
b
-
-
-
0.45
0.5
0.6Table 5
Single-shot kill probability.
W
1
W
2
W
3
W
4
W
5
W
6
T
1
0.45
0.32
0.54
0.47
0.75
0
T
2
0.5
0.27
0.34
0.73
0.78
0
T
3
0.28
0.62
0.58
0.56
0.57
0
T
4
0
0
0
0.63
0.55
0.7
T
5
0
0
0
0.45
0.5
0.85


### 6.2.1. Analyze the Result of MOQPSO-D/S for Solving Uncertain WTA Problem

With different population size and different iterations, the WTA case is solved with application of MOQPSO-D/S algorithm. The algorithm for solving the WTA case under every condition runs 30 times independently. The mean running time is listed in Table6. It can be seen from Table 6 that the real-time performance of MOQPSO-D/S algorithm can meet the need of WTA in air and missile defense combat.Table 6
Mean running time/s.
Population size
Iterations
50
100
200
50
0.7504
1.4830
2.9733
100
1.5179
2.9313
5.9883
200
3.0576
6.0737
11.8876With the population size and iterations both being 50, the solution of the WTA case solved by MOQPSO-D/S algorithm is shown in Figure5. It can be seen from Figure 5 that the Pareto optimal solution, obtained by MOQPSO-D/S algorithm for solving the fuzzy multiobjective programming based WTA model in air and missile defense combat, distributes uniformly, where every dot refers to a WTA alternative. Therefore, it can provide a series of alternatives for commanders to make decision. The commanders could make the final decision by choosing a suitable WTA alternative from the solution with consideration of the battlefield situation and risk preference. Comparing with the single objective optimization based WTA method, which can only give a WTA alternative, this method could preferably adapt to the change of the battlefield situation and combine the tactical thought of the commander to make decision.Figure 5
Solution of MOQPSO-D/S for solving the case.


### 6.2.2. Compare with the Single Objective Optimization Based WTA Method

In order to validate the feasibility of the solution obtained by MOQPSO-D/S algorithm for solving WTA model, in comparison, the example is solved by adopting the method proposed in literature [6], taking maximizing effectiveness-cost ratio as the objective function and adopting IDPSO algorithm to solve the WTA model. With the population size and iterations being both 50, the single objective optimization method for solving the WTA case runs 30 times. The mean running time is 0.7193 s. The best fitness over generations based on IDPSO is shown as in Figure 6. The optimal WTA result is (31)T1T2T3T4T5W1W2W3W4W5W6001000200000001000000200000011.Figure 6
Best fitness over generations based on IDPSO algorithm.The optimal result is just the WTA alternative corresponding to the dot marked by a red five-angle star shown in Figure5. And the mean running time is a little shorter than the MOQPSO-D/S algorithm with the same population size and iterations. It is proved that MOQPSO-D/S algorithm, when solving the WTA model, would disperse computation resource but can still obtain the satisfactory solution containing the optimal result obtained by single optimization method. Meanwhile, the fuzzy multiobjective programming WTA model for air and missile defense is proved to be reasonable and adopting MOQPSO-D/S algorithm to solve the model is feasible.


### 6.2.3. Compare with Other Multiobjective Optimization Algorithm for Solving WTA Method

In order to further verify the performance of MOQPSO-D/S algorithm for solving WTA problem, in comparison, the example is solved by adopting MOQPSO-DPS, MOQPSO-CD, and MOPSO-CD algorithm with the population size and iterations being both 100. Every algorithm for solving the WTA case runs 30 times independently. The statistical results for running time of every algorithm are shown in Figure7. It can be seen from Figure 7 that the real-time performance of MOQPSO-D/S algorithm is obviously superior to the comparable algorithms’.Figure 7
Statistical results for running time of four algorithms.The convergence of every algorithm for solving the WTA case is weighed by dominated ratioC(X1,X2) [25]. The calculation of C(X1,X2) is shown as in (32). X1 and X2 denote the Pareto optimal solution sets obtained by two algorithms, respectively. If C(X1,X2) is bigger than C(X2,X1), the convergence of the algorithm corresponding to X1 is superior to the algorithm corresponding to X2. The convergence comparison among four algorithms is listed in Table 7. A1,A2,A3, and A4 listed in Table 7 denote the Pareto optimal solution sets obtained by MOQPSO-D/S, MOQPSO-DPS, MOQPSO-CD, and MOPSO-CD algorithm, respectively. It can be concluded from Table 7 that the convergence of MOQPSO-D/S algorithm for solving the WTA case is the best of four algorithms.(32)CX1,X2=numelx2∈X2∣∃x1∈X1:x1≺x2numelX2.In (32), the function numel(•) is used to count the number of elements in a set.Table 7
Dominated ratio among four algorithms’ Pareto optimal sets.
X
2
X
1
A1
A2
A3
A4
A1
0
0
0.0185
0.0370
A2
0.0588
0
0.1176
0.1765
A3
0.0556
0
0
0
A4
0.0755
0.0755
0.0755
0The distribution uniformity of the Pareto optimal solution set obtained by every algorithm for solving the WTA case is weighed by SP. The statistical results for SP value of the Pareto optimal solution set obtained by every algorithm are shown in Figure8. It can be concluded from Figure 8 that the Pareto optimal solutions obtained by MOQPSO-D/S algorithm for solving the WTA case distribute the most uniformly, and its distribution uniformity is the most stable.Figure 8
Statistical results for SP value of four algorithms.


## 6.2.1. Analyze the Result of MOQPSO-D/S for Solving Uncertain WTA Problem

With different population size and different iterations, the WTA case is solved with application of MOQPSO-D/S algorithm. The algorithm for solving the WTA case under every condition runs 30 times independently. The mean running time is listed in Table6. It can be seen from Table 6 that the real-time performance of MOQPSO-D/S algorithm can meet the need of WTA in air and missile defense combat.Table 6
Mean running time/s.
Population size
Iterations
50
100
200
50
0.7504
1.4830
2.9733
100
1.5179
2.9313
5.9883
200
3.0576
6.0737
11.8876With the population size and iterations both being 50, the solution of the WTA case solved by MOQPSO-D/S algorithm is shown in Figure5. It can be seen from Figure 5 that the Pareto optimal solution, obtained by MOQPSO-D/S algorithm for solving the fuzzy multiobjective programming based WTA model in air and missile defense combat, distributes uniformly, where every dot refers to a WTA alternative. Therefore, it can provide a series of alternatives for commanders to make decision. The commanders could make the final decision by choosing a suitable WTA alternative from the solution with consideration of the battlefield situation and risk preference. Comparing with the single objective optimization based WTA method, which can only give a WTA alternative, this method could preferably adapt to the change of the battlefield situation and combine the tactical thought of the commander to make decision.Figure 5
Solution of MOQPSO-D/S for solving the case.


## 6.2.2. Compare with the Single Objective Optimization Based WTA Method

In order to validate the feasibility of the solution obtained by MOQPSO-D/S algorithm for solving WTA model, in comparison, the example is solved by adopting the method proposed in literature [6], taking maximizing effectiveness-cost ratio as the objective function and adopting IDPSO algorithm to solve the WTA model. With the population size and iterations being both 50, the single objective optimization method for solving the WTA case runs 30 times. The mean running time is 0.7193 s. The best fitness over generations based on IDPSO is shown as in Figure 6. The optimal WTA result is (31)T1T2T3T4T5W1W2W3W4W5W6001000200000001000000200000011.Figure 6
Best fitness over generations based on IDPSO algorithm.The optimal result is just the WTA alternative corresponding to the dot marked by a red five-angle star shown in Figure5. And the mean running time is a little shorter than the MOQPSO-D/S algorithm with the same population size and iterations. It is proved that MOQPSO-D/S algorithm, when solving the WTA model, would disperse computation resource but can still obtain the satisfactory solution containing the optimal result obtained by single optimization method. Meanwhile, the fuzzy multiobjective programming WTA model for air and missile defense is proved to be reasonable and adopting MOQPSO-D/S algorithm to solve the model is feasible.


## 6.2.3. Compare with Other Multiobjective Optimization Algorithm for Solving WTA Method

In order to further verify the performance of MOQPSO-D/S algorithm for solving WTA problem, in comparison, the example is solved by adopting MOQPSO-DPS, MOQPSO-CD, and MOPSO-CD algorithm with the population size and iterations being both 100. Every algorithm for solving the WTA case runs 30 times independently. The statistical results for running time of every algorithm are shown in Figure7. It can be seen from Figure 7 that the real-time performance of MOQPSO-D/S algorithm is obviously superior to the comparable algorithms’.Figure 7
Statistical results for running time of four algorithms.The convergence of every algorithm for solving the WTA case is weighed by dominated ratioC(X1,X2) [25]. The calculation of C(X1,X2) is shown as in (32). X1 and X2 denote the Pareto optimal solution sets obtained by two algorithms, respectively. If C(X1,X2) is bigger than C(X2,X1), the convergence of the algorithm corresponding to X1 is superior to the algorithm corresponding to X2. The convergence comparison among four algorithms is listed in Table 7. A1,A2,A3, and A4 listed in Table 7 denote the Pareto optimal solution sets obtained by MOQPSO-D/S, MOQPSO-DPS, MOQPSO-CD, and MOPSO-CD algorithm, respectively. It can be concluded from Table 7 that the convergence of MOQPSO-D/S algorithm for solving the WTA case is the best of four algorithms.(32)CX1,X2=numelx2∈X2∣∃x1∈X1:x1≺x2numelX2.In (32), the function numel(•) is used to count the number of elements in a set.Table 7
Dominated ratio among four algorithms’ Pareto optimal sets.
X
2
X
1
A1
A2
A3
A4
A1
0
0
0.0185
0.0370
A2
0.0588
0
0.1176
0.1765
A3
0.0556
0
0
0
A4
0.0755
0.0755
0.0755
0The distribution uniformity of the Pareto optimal solution set obtained by every algorithm for solving the WTA case is weighed by SP. The statistical results for SP value of the Pareto optimal solution set obtained by every algorithm are shown in Figure8. It can be concluded from Figure 8 that the Pareto optimal solutions obtained by MOQPSO-D/S algorithm for solving the WTA case distribute the most uniformly, and its distribution uniformity is the most stable.Figure 8
Statistical results for SP value of four algorithms.


## 7. Conclusions

The goal of this article is to find a method to aid decision-making of WTA in air and missile defense combat under uncertainty. A fuzzy multiobjective programming based WTA model has been constructed. And multiobjective optimal algorithm, called MOQPSO-D/S, has been put forward to solve the fuzzy multiobjective programming based WTA model. The simulation results shows that the performance of MOQPSO-D/S algorithm is superior, the WTA model constructed in this paper is rational and the MOQPSO-D/S algorithm can solve the WTA model efficiently.In future research, we may study WTA problem under uncertainty based on other uncertain methods, for example, the Taguchi method and the Bayesian method. Yet, we may study multiobjective optimization algorithms, for example, the game theory based methodologies, to solve the multiobjective programming based WTA model.



---

*Source: 9897153-2017-12-14.xml*