Solving the Dynamic Weapon Target Assignment Problem by an Improved Multiobjective Particle Swarm Optimization Algorithm

: Dynamic weapon target assignment (DWTA) is an effective method to solve the multistage battleﬁeld ﬁre optimization problem, which can reﬂect the actual combat scenario better than static weapon target assignment (SWTA). In this paper, a meaningful and effective DWTA model is established, which contains two practical and conﬂicting objectives, namely, maximizing combat beneﬁts and minimizing weapon costs. Moreover, the model contains limited resource constraints, feasibility constraints and ﬁre transfer constraints. The existence of multi-objective and multi-constraint makes DWTA more complicated. To solve this problem, an improved multiobjective particle swarm optimization algorithm (IMOPSO) is proposed in this paper. Various learning strategies are adopted for the dominated and non-dominated solutions of the algorithm, so that the algorithm can learn and evolve in a targeted manner. In order to solve the problem that the algorithm is easy to fall into local optimum, this paper proposes a search strategy based on simulated binary crossover (SBX) and polynomial mutation (PM), which enables elitist information to be shared among external archive and enhances the exploratory capabilities of IMOPSO. In addition, a dynamic archive maintenance strategy is applied to improve the diversity of non-dominated solutions. Finally, this algorithm is compared with three state-of-the-art multiobjective optimization algorithms, including solving benchmark functions and DWTA model in this article. Experimental results show that IMOPSO has better convergence and distribution than the other three multiobjective optimization algorithms. IMOPSO has obvious advantages in solving multiobjective DWTA problems.


Introduction
Weapon target assignment (WTA) is a key issue of Command & Control (C2) [1] in modern warfare. Its goal is to strike enemy targets by rationally allocating weapon units, optimize the firepower strike system, and achieve the best strike effect. As an important subject of national defense-related applications, WTA has become a hot focus at present [2][3][4][5][6], attracting a large number of scholars to study, which has important military significance [7].
WTA has been proved to be NP-complete [8], and the amount of calculation for solving WTA increases exponentially with the increase of dimension. Hosein et al. divided WTA problems into two versions, including static weapon-target assignment (SWTA) and dynamic weapon-target assignment (DWTA) [9]. SWTA is the process of distributing all weapons to targets in a single stage. DWTA is a multi-stage decision-making problem, which incorporates the concept of time window. Time window refers to the time interval from when a target arrives at the weapon combat zone to when the target escapes from the weapon combat zone. As the battlefield situation is constantly changing, each weapontarget-stage pair will be assigned a corresponding time window. DWTA requires that the change of battlefield situation information should be considered when assigning tasks at different stages. The DWTA algorithm must give the result of allocation scheme in a limited time, which can provide the basis for subsequent decision-making.
In recent years, there have been many methods to solve WTA problems, including traditional methods and intelligent algorithms. Traditional methods include goal programming [10], game theoretic framework [11], Lagrangian relaxation method [12], etc. These algorithms can solve WTA problems with smaller dimensions, but their global optimization ability is poor, which are unable to meet the actual needs. Intelligent algorithm is widely used to solve WTA problems because of its high optimization precision and fast convergence speed. Hongtao et al. [13] put forward a new selection algorithm, which used adaptive chaos optimization and clonal selection method to solve WTA. Sounc et al. [14] used a parallel simulated annealing algorithm to solve the WTA problem. Chen et al. [15] proposed an evolutionary decision algorithm including genetic algorithm and memetic algorithm to solve the DWTA problem. Li et al. [16] introduced the population initialization method for specific problems into evolutionary algorithm, improved the selection operation mechanism, and improved the efficiency of solving the SWTA problem. Hu et al. [17] proposed an ant colony optimization algorithm based on elite strategy, which improved the strategies of path selection and pheromone update and solved four WTA models. Although these methods have a certain improvement in optimization accuracy compared with traditional methods, it is necessary to further improve the optimization ability and solution quality of these algorithms to meet the operational requirements of the battlefield.
WTA is a typical combinatorial constraint optimization problem. Generally, the evaluation index of WTA problem adopted is the overall damage probability of the target. An excellent WTA scheme should not only achieve the purpose of combat, but also minimize the cost of combat. Because the smaller the ammunition cost and the better the combat effect are two conflicting goals, the actual WTA is a multiobjective optimization problem (MOP). Although some algorithms are weighted multiple evaluation indexes into one cost function, these algorithms still aim to optimize a single objective. The defect of this method is that the weight of the objective is usually generated through experience, and improper weight selection will seriously affect the optimization effect of the algorithm and may produce a decision scheme that the commander does not really expect. Many multiobjective optimization techniques have been studied to solve multiobjective WTA optimization problems over the last decade. Li et al. [18] solved the target-based multiobjective optimization DWTA problem by using NSGA-II. Li et al. [19] proposed an adaptive mechanism to solve the multi-stage weapon target allocation problem and compared the adaptive NSGA-II and adaptive MOEA/D with other algorithms [20]. Li et al. [21] proposed an improved Pareto ant colony optimization method to solve the target-based multiobjective SWTA optimization problem. All these researches show the advantage of solving WTA with multiobjective optimization algorithm.
The purpose of this paper is to use multiobjective optimization algorithm to solve DWTA problem. At present, the well-known multiobjective optimization algorithms are SPEA-II [22], PESA-II [23], NSGA-II [24], MOEA/D [25], MOPSO [26], etc. Among these algorithms, MOPSO is simple in principle and easy to implement, so it is used to solve multiobjective optimization problems in engineering application. The MOPSO algorithm was proposed by C.A.Coello et al. in the early 21st century. Since then, many scholars have done a lot of research on MOPSO and put forward many excellent MOPSO algorithms [27][28][29][30][31]. The purpose of these algorithms is to obtain more accurate Pareto front with good distribution characteristics when solving MOPs. MOPSO has many advantages that traditional optimization algorithms cannot achieve. However, many researchers find that MOPSO has some defects in solving problems of engineering applications, the most notable of which is that MOPSO is easy to converge in advance, leading to falling into local optimal solution.
In accordance with the characteristics of DWTA and the weakness of MOPSO in solving DWTA problems, we propose an improved multiobjective particle swarm optimization algorithm (IMOPSO). The main contributions of the algorithm are described as follows: (1) A novel particle learning strategy is proposed. According to the evolutionary characteristics of dominated and non-dominated solutions in MOPSO, various learning strategies are adopted, which makes the algorithm can learn and evolve in a targeted way. (2) A novel search strategy based on simulated binary crossover and polynomial mutation is proposed, which enables elitist information to be shared among external archive and enhances the exploratory capabilities of IMOPSO. (3) A dynamic external archive maintenance strategy is proposed. Compared with earlier archive maintenance methods, a dynamic external archive maintenance strategy can improve the diversity of non-dominated solutions.
The rest of this paper is organized as follows. In Section 2, some basic definitions of the MOP, the DWTA problem and its mathematical model are introduced and analysed. In Section 3, MOPSO and the improved method are described. In Section 4, the experimental research is introduced. In Section 5, the experimental results are analyzed. The experimental results are discussed in Section 6. The conclusion and future work are described in Section 7.

Related Concepts of MOP
The mathematical model of multiobjective optimization problem (MOP) can be expressed as follows: where Ω represents the decision space, x is the decision vector, x = [x 1 , x 2 , · · · , x n ], n is the total number of decision variables. f m (x) denotes the mth objective function.

F(x)
: Ω → R m is composed by m objective functions. Given two vectors u = (u 1 , u 2 , · · · u m ) and v = (v 1 , v 2 , · · · v m ) ∈ Ω, we say that u dominate v (denoted by u ≺ v), if and only if: u i ≤ v i for all i ∈ {1, 2, · · · m} and u j ≤ v j for at least one component. This characteristic is called as Pareto Dominance.
Also, only if there is no vector x ∈ Ω, which makes f(x) ≺ f(x * ), then the decision vector x * Ω can be called Pareto optimal solution or non-dominated solution.
The Pareto optimal solution set within Ω is called Pareto optimal set, which is denoted as P * . Furthermore, the Pareto front can be defined as PF * = {f(x) : x ∈ P * }.

Analysis of DWTA Problem
DWTA is a multistage problem, and it requires decision makers assess available weapons and targets to be attacked for subsequent decisions. The most prominent feature of DWTA is that the continuously changing battlefield situation information should be considered in the process of target assignment in different stages. DWTA should be adapted to the dynamic solving conditions to obtain the optimized attack scheme. The actual DWTA is composed of several stages. The results of the current stage need to be evaluated after the decision is completed, and the results of the current stage will have an important impact on the decision-making of the next stage. This is a cyclic process: Decision Making → Damage Assessment → Decision Making. In this paper, the DWTA problem is discretized into a multi-stage SWTA process, and the time window is considered at the same time, so DWTA can be expressed as follows.
Based on the above analysis, the combat scenario considered in this paper is described as follows. Suppose that in a certain period of time, the defender detects T hostile targets, and the defender has W weapons to intercept targets. Before the hostile targets break through the defense or the targets are completely destroyed by firepower, the defender has at most S stages of shooting targets with weapons.

Expected Damage of Targets
The first objective is very common in WTA-that is, maximizing the damage to hostile targets, which is also the most basic operational intention of WTA. The expected damage formula at stage t is as follows: where t and s are parameters related to the time window, W(t) and T(t) represent the remaining number of weapons and targets in stage t, (W(1) = W, T(1) = T) respectively. v j is the threat value of target j. p ij (s) means the kill probability of weapon i destroys target j in stage s, and the kill probability matrix is P(s) = p ij (s) T×W . X t = [X t , X t+1 , · · · , X s ] with X t = x ij (t) W×T is the decision matrix at stage t, and x ij (t) is a binary variable, which is represented by

Ammunition Consumption
In addition to meeting tactical requirements, a WTA decision should also consider operational costs. The second objective function is to minimize the ammunition consumption of weapons at all stages. The corresponding formulation is expressed as: where u ij (s) is the ammunition consumption allocated by weapon i to target j at stage s, β i is the ammunition unit cost of weapon i. We assume that all weapons are the same, and β i is assumed to be constant. The ammunition consumption matrix is U(s) = u ij (s) T×W .
In Constraint set (6), m j represents the maximum quantity of weapons consumed when hitting target j at stage t. This constraint limits the maximum amount of ammunition used by each target to be destroy at each stage. Constraint set (7) indicates that weapon i can shoot at most n i targets at a time, which is related to the property of the weapon itself. Generally, most weapons can fire only one target at one time; that is, n i = 1 with ∀i ∈ {1, 2, · · · , W}. Constraint set (8) indicates the number of ammunition available for weapon i at all stages.

Feasibility Constraints
x ij (t) ≤ f ij (t), ∀i ∈ {1, 2, · · · , W}, ∀j ∈ {1, 2, · · · , T}, ∀t ∈ {1, 2, · · · , S} where f ij (t) is a 0 or 1 variable. If the target j cannot be shot because it exceeds the maximum range of weapon i at stage t, then f ij (t) is 0; otherwise, it is 1. The feasibility constraint has the time attribute of DWTA, the influence of time window on the engagement feasibility is considered, so the feasibility matrix F(s) = f ij (s) W×(T×S) needs to be updated after each stage.

Fire Transfer Constraints
Usually, a weapon needs to adjust its shooting orientation and reload its ammunition after one shooting and before the next shooting, which is widely exist in WTA problems. The fire transferring time is defined as the time interval between when a weapon finished shooting from one direction and started shooting from another direction. If the time interval between two consecutive shots is less than the corresponding fire transfer time, then the weapon can only complete one shot.
In order to introduce the constraints of fire transfer, we set a combat scenario including W weapons, T targets and S stages, so the fire transferring feasibility can be expressed as: where f t f w i , S m , S n , t j , t k is a binary variable. f t f w i , S m , S n , t j , t k = 0 indicates that if weapon w i is able to shoot target t j at stage S m (i.e., x w i t j (S m ) = 1), it cannot complete the shooting of target t k at stage S n (i.e., x w i t k (S n ) = 0) and vice versa. f t f w i , S m , S n , t j , t k = 1 means no such restrictions. According to the fire transfer feasibility matrix, the fire transfer constraint is as follows: which indicates that weapon w i has a fire transfer constraint for targets t j at stage S m and target t k at stage S n . Based on the above description, the multiobjective DWTA model is constructed as follows: (7), (8), (9), (10).

Multiobjective Particle Swarm Optimization Algorithm
The MOPSO algorithm is the name of the PSO algorithm used to solve multiobjective optimization problems. Both update the velocity and position of the particle in the same way. The PSO algorithm updates the velocity and position of a particle through two "guides": the optimal solution found by the particle (i.e., the individual guide p best ) and the optimal solution found by the entire population thus far (i.e., the global guide g best ). In the PSO algorithm, the velocity and position update with the following conditions: x id (t + 1) = x id (t) + v id (t + 1) In the d-dimension direction of search space, v id (t) and x id (t) mean the velocity and position of the i-th particle in the t-th iteration respectively, c 1 and c 2 are acceleration factors of individual and global optimal particles, respectively, r 1 and r 2 are random numbers between [0, 1], p ibestd is the individual optimal position found by the i-th particle until the t-th iteration, g bestd is the global optimal position found by the entire population until the t-th iteration, v id (t + 1) and x id (t + 1) indicate the updated velocity and updated position of the i-th particle obtained after the t + 1 iteration, respectively.
Although there is no difference between MOPSO and PSO in formula expression, MOPSO is different from PSO in the following ways. Firstly, MOPSO is suitable for solving MOPs, so the updating and selection methods guided by individuals and globally are quite different from PSO. In addition, the solution obtained by MOPSO is a non-inferior solution. Moreover, MOPSO algorithm has an archive dedicated to storing the global optimal solution currently found by particles. Researchers have studied and improved MOPSO from many aspects and made some meaningful achievements [32][33][34]. Even though MOPSO has been widely used in practical problems, there are still problems such as premature convergence and poor accuracy in searching global optimal solution.

Improvement of Learning Strategies
When MOPSO is used to solve the MOPs, its method of selecting the global optimal solution is to select it randomly from the external archive. However, due to the lack of pertinence, the guiding role of g best cannot be effectively played. In this paper, an improved learning strategy is proposed to deal with g best in different situations, so that it can be better applied to solve MOPs.
(1) If x i (t) is not included in the external archive, then the individuals that dominate x i (t) in the archive can be used as its global optimal solution. Therefore, it can be considered that g best of the solution is a linear combination of all the archive individuals that dominate x i (t).
It is assumed that there are k solutions in the external archive that dominate x i (t): a 1 (t), a 2 (t), · · · a k (t). We obtain a set of weights w j by random generation, so that ∑ k j=1 w j = 1. Thus, g best can be expressed as follows: The velocity updating formula of the particle is expressed as follows: (2) If x i (t) is the solution in the external archive, then there is no individual that dominates x i (t), and it is a non-inferior solution in itself. It can be considered that there is no global optimal solution, that is, g best does not exist. The velocity updating formula of the particle is expressed as follows: In the above two cases, only the velocity is defined differently, and the position updating method of particles remains unchanged.

The Search Strategy Based on Crossover and Mutation
The genetic algorithm optimizes the objective by simulating the law of survival of the fittest in the biological evolution in nature [35]. In the genetic algorithm, two individuals are randomly selected from the population as parents, crossover operation is performed with a higher probability, mutation operation is performed with a lower probability for each individual, then the better individual is selected to enter the next generation, and the optimal solution is found through such iterative process [36]. Inspired by the above method, this paper proposes a search strategy based on Simulated Binary Crossover (SBX) and Polynomial Mutation (PM).
In each iteration, in order to make better use of the dimensional information of p best , the external archive is used to store the generated p best , and then the p best in the external archive is applied to generate new particles through the SBX criterion. The new particles merge all dimensional information of p best generated by each iteration, which enhances the communication of all dimensional information among particles and effectively improves the ability of particles to jump out of the local optimal solution.
If a population stops evolving (i.e., the velocity of all particles in the population is almost equal to zero), the global optimal solution will not be obtained. In this article, polynomial mutation is used to generate a new global optimal particle (g best ), which makes the current population jump out of the local optimal state. Different from the usual mutation strategy (i.e., one particle mutates with a certain probability), the strategy we propose is as follows: if the g best of the particle has not been improved for several consecutive generations, the newly generated g best is mutated by polynomial.
SBX operator allows elite solutions to exchange useful gene segments, while PM operation adds a small turbulence to search local region. In order to execute SBX and PM operators on the external archive A, an elite subset E is firstly selected from A, which contains a number of non-dominated solutions with bigger crowding distance in A. Generally, the size of E is smaller than that of A, which is set as half of |A| in this paper. For each solution A i (i = 1, 2, · · · , |A|), a random integer j is generated in [1, |E|]. Then, A i and E j are used as the parent solutions for performing SBX operator. A child solution is randomly selected from SBX operator, and then the PM operator is proceeded to perform.

Maintenance Strategy of External Archive
The non-inferior solution obtained by MOPSO algorithm is stored in the external archive. Due to the consideration of storage capacity and computational efficiency, the archive cannot be infinite. In order to avoid the maximum size of external archive, it is necessary to delete the most crowded individuals in the external archive. Most optimization algorithms use the crowding distance in NSGA-II to indicate the crowding degree [24].
Although NSGA-II can maintain the external archive to a certain extent, it also has the defect of reducing the diversity of non-dominated solutions. Let's illustrate with an example. If there is a given set of non-dominated solutions as shown in Figure 1a, two of these solutions are redundant. Since solution C and solution D are the two solutions with the minimal crowding distance in the set of solutions, according to the processing method of NSGA-II, these two solutions need to be deleted at the same time. The result is shown in Figure 1b. We can intuitively see from the figure that the solution obtained by this method has poor distribution. In view of the above problems, this paper adopts a dynamic maintenance strategy to maintain external archive. The way to deal with crowding distance in dynamic maintenance strategy is to delete only the solution with the smallest crowding distance in the external archive every time, then identify those individuals whose order has been changed in the solution set and recalculate their crowding distance values. Repeat this process several times to ensure that the number of solutions in the external archive does not exceed the maximum threshold. Based on the dynamic maintenance strategy, solution C is deleted first, then calculate the crowding distance between B and D, and then solution E is deleted instead of solution D. The maintenance result is shown in Figure 1c. Comparing Figure 1c with Figure 1b, it is obvious that the solution set generated by dynamic maintenance strategy is more evenly distributed.

Updating Method of External Archive
The way to update external archive by MOPSO is to compare the current population with the original archive by using the archive update rules to generate a new external archive. In this paper, we use the two sets to update the original archive together, which can further improve the performance of the external archive. The specific updating methods are shown in Figure 2.

Elite Individuals Increasing Strategy
Fundamentally, both the parent population and the archive will have an impact on the generation of the offspring population. In order to increase the rate of the population convergence, this paper presents a strategy of elitist incremental, which forms a new offspring population by adding archive individuals in the offspring population. The corresponding strategy is defined as follows: A new offspring population = randomly selected (N − L) offspring solutions +randomly selectedLexternalarchivesolutions (16) N in the formula represents the population size. T is the number of iterations. |A| is the archive size of the algorithm for T iterations.
The above method can further improve the search performance of the algorithm. First, the population has excellent exploration ability at the early stage and can find a more complete range of non-inferior solutions. Secondly, at the later stage, the population can also search near the known non-inferior solution to make it closer to the true PF.

The IMOPSO Algorithm Flow
Based on the analysis of improvement strategies, the process of IMOPSO is described in Algrorithm 1 as follows: Algorithm 1. IMOPSO Input: r 1 , r 2 , c 1 , c 2 , w, T max , the population size N, the maximum archive size A max Output: The Pareto front solutions Step 1. Particle position and velocity are initialized by random generation, and related parameters r 1 , r 2 , c 1 , c 2 and w of IMOPSO are set. Set the maximum number of iterations T max , population size N, and the maximum size of external archive A max , and set T = 0.
Step 2. According to the set parameters, the fitness value of each individual in the initial population is calculated, and the optimal position p ibest of the initial individual is obtained.
Step 3. Update the external archive A.
Step 4. For T = 1: (1) Use Equations (14) or (15) to update the particle velocity and use Equation (12) to update the particle position to form the middle population.
(2) The update procedure of archive A is executed to gather the new non-dominated solutions with bigger crowding-distance values.
(3) Use the method in Section 3.2.2 to perform the evolutionary search process.

Experimental Studies
To verify the performance of the IMOPSO algorithm, the following two experiments are conducted in this paper: (1) Performance experiment of IMOPSO in solving benchmark functions Eight benchmark functions are solved by using the IMOPSO algorithm, and the optimization effect of IMOPSO is verified by comparing with algorithms of MOPSO, NSGA-II and MOEA/D.
(2) Performance experiment of IMOPSO in solving DWTA problem Using the IMOPSO algorithm to solve the established DWTA model, and comparing with the MOPSO, NSGA-II, MOEA/D algorithms to verify the superiority of the IMOPSO algorithm.
The above experimental process was implemented on a personal computer with 64-bit Windows 10, an Intel Core i7-8700K @ 3.70GHz processor, and 16 GB of RAM.

Parameters Setting of Algorithms
To verify the performance of IMOPSO, we selected eight unconstrained benchmark functions; namely, ZDT1-ZDT6, Schaffer (SCH), Kursawe (KUR) and Fonseca (FON). Considering that the DWTA model established in this paper is a bi-objective optimization problem, the benchmark functions selected here are also bi-objective. These test functions are characterized with bi-objective, discontinuous, concave and convex, etc. They provide enough difficulty for testing multiobjective optimization algorithm in searching Pareto optimal solution, challenging algorithm convergence and diversity of solution set distribution. The related information of benchmark functions is shown in Table 1.
Parameters of each algorithm are set as follows. The population size of the four algorithms is set to 100. The crossover probability of NSGA-II is 0.9, tournament selection is adopted, and the mutation probability is 1/n (where n is the number of decision variables). The external archive size of IMOPSO and MOPSO is 100, the inertia weight w decreases from 0.9 to 0.4, and c 1 = c 2 = 2.0. The crossover probability and mutation probability of MOEA/D and NSGA-II are the same, which are 0.9 and 1/n, respectively. The neighbour size of MOEA/D is 30, and the decomposition method is tchebucheff. The number of iterations of ZDT1 and ZDT2 are 60. Since ZDT3, ZDT4 and ZDT6 functions are relatively difficult to optimize, the number of iterations is set to 100. The iteration number of SCH is 60, and that of FON and KUR is 100.

Performance Metrics
To comprehensively evaluate the performance of multiobjective optimization algorithms, it is necessary to analyze them from three aspects: convergence performance, diversity and width of solution set distribution [37]. Therefore, the following three measurement metrics are used to compare the multiobjective optimization methods: Closeness to PF * Generational distance (GD) [38] metric is used to measure the distance between the non-dominated solution set (PF know ) searched by the optimization algorithm and the true Pareto front (PF * ), which is described as follows: where n represents the number of individuals in PF know , and D i refers to the Euclidean distance between the i-th individual in PF know and the nearest individual in PF * . The smaller GD value indicates the better convergence performance. Diversity In 2002, Deb [39] introduced an index to evaluate the diversity of distribution of non-inferior solution sets, and its formula is defined as follows: Parameters D f and D l in the formula are Euclidean distances of extreme solutions in true Pareto front set PF * and boundary solutions of non-dominated solution set PF know obtained by multiobjective optimization algorithm. D i is the distance between two adjacent solutions in the PF know , and D is the average of all the D i . A smaller value of ∆ means that the better the performance of the algorithm.
Distribution width Zitzler et al. [37] proposed a metric to measure the distribution width of non-inferior solution sets, and its calculation formula is as follows: where NP is the set of non-inferior solutions and M is the dimension of solutions. The larger the value of M * 3 , the wider the distribution of non-inferior solutions. However, if the Pareto front is widely distributed, the dimension of the solutions is too large, or the non-dominated solution does not converge, the value of M * 3 will also increase, which is not convenient for comparison and collation of algorithm results. Therefore, this paper adopts the following measurement method of distribution width: where RP is the Pareto optimal set or the reference point set. Notably, a small value of DW reflects a better distribution uniformity for the non-dominated solution set. Figure 3 shows the decision-making process of the DWTA problem. We can see from the figure that DWTA decision-making needs to determine the targets to be attacked and available weapons at the current stage on the basis of the results of battlefield situation analysis. Then, the determined targets and weapons are input into the IMOPSO algorithm to perform optimization iteration to obtain WTA decision. DWTA decision-making process also needs to evaluate and monitor the allocation results and battlefield situation, and determine whether iteration needs to be continued in accordance with the effect of the allocation scheme and battlefield situation.

Integer Encoding Method
The allocation scheme obtained by IMOPSO algorithm in solving multiobjective DWTA can be obtained by decision matrix. To complete this process, we need to encode the particles. In order to obtain an efficient WTA scheme, the coding method can not only meet the solution of the model, but also meet the constraint requirements. Moreover, the encoding length should be made as small as possible. Based on the above conditions, this paper adopts integer encoding method. Take the following scenario (W = 5, T = 6, S = 1) as an example to illustrate the coding method shown in Figure 4, where X is the corresponding decision matrix.

Setting of Experimental Parameters
To verify the optimization performance of the IMOPSO algorithm in solving multiobjective DWTA problem, this section also uses IMOPSO, MOPSO, NSGA-II and MOEA/D for comparative analysis. Three different combat scenarios, i.e., small-, medium-and large-scaled cases are randomly generated in line with the number of weapons and targets. For each case, we generate the killing probability matrix P(s), ammunition consumption matrix U(s), feasibility matrix F(s) and fire transferring feasibility matrix FTF(s), and the elements in the matrix are randomly generated. The related parameters of DWTA problems are as follows: the number of weapons/targets/stages are 5/8/3, 15/24/5 and 30/50/8 in small/medium/large-scale cases, respectively. In all three cases, n i is 1, N i is 2. m j is 1 in small-scale case, and m j is 2 in medium-scale and large-scale case. The value vector of targets comes from Ref [40]. The number of iterations of all algorithms in the three scales are 100,150,200, respectively. The crossover probability and mutation probability of NSGA-II and MOEA/D are 0.9 and 0.1, respectively. Other parameters of the four algorithms are the same as those of solving the benchmark functions.

Performance Metrics
In this paper, the inverted generational distance (IGD) is used to measure the performance of multiobjective optimization algorithms in solving DWTA problem [41]. IGD is a variant of GD, which can not only evaluate the convergence of multiobjective solutions, but also evaluate the distribution characteristics of solution set to a certain extent. The calculation formula is as follows: where P * is set of uniformly distributed points in the objective space along the true Pareto front (PF) or nearly true PF when it is hard or impossible to get the true PF, P is an approximate PF of the true PF, |P * | represents the number of points in P * , and dist P * i , P is the minimum Euclidean distance between P * i and elements in P. The smaller the value of IGD, the better the convergence and distribution of the searched solution.
As to multiobjective DWTA problems, it is hard to obtain the true PF. To address this issue, we run each instance for 30 times, execute the non-dominated sorting on these solutions obtained earlier and select 500 evenly distributed points in Pareto front as the true Pareto front P * .

Simulation Results of Benchmark Functions
In order to eliminate the influence of the randomness of all algorithms on the simulation results, the four algorithms were independently run 30 times on all benchmark functions, and the results were statistically analyzed. The experimental results are shown in Tables 2-4, and the corresponding results of the maximum, minimum, mean and standard deviation are listed in above tables. The best results are bolded to facilitate the visualization.     Table 2 shows the statistical results of the GD metric. It can be seen that the IMOPSO algorithm has obtained 28 optimal values out of 32 index statistics. For the function KUR, the GD of IMOPSO has the largest standard deviation, which indicates that IMOPSO does not perform well in convergence stability. With respect to the other three algorithms, MOPSO has better convergence performance in functions ZDT1, ZDT2 and ZDT3, but its convergence stability is not as good as NSGA-II. MOEA/D has better performance in the GD metric in functions ZDT4, ZDT6, SCH and FON than MOPSO and NSGA-II. Table 3 shows the results of ∆ indicator. We can see that IMOPSO has achieved the best average value in all eight test problems, which shows that the non-dominated solution generated by IMOPSO has the best distribution. However, it can be seen from the standard deviation of some functional indicators that IMOPSO is inferior to other algorithms in the ∆ metric. The distribution performance of NSGA-II in the functions ZDT1, ZDT2 and ZDT3 is second only to that of IMOPSO. MOEA/D has obtained better ∆ values in the functions ZDT4, ZDT6 and FON, while MOPSO has the worst ∆ performance and has not gained any optimal value. Table 4 is the statistical data of distribution width index DW. It can be seen from the data that IMOPSO is the best among the four algorithms except the minimum values of ZDT3, ZDT4, ZDT6 and KUR, and the data stability is also the best. The index performance of NSGA-II is second only to IMOPSO. Compared with MOPSO and MOEA/D, it has obtained the best values in functions ZDT1, ZDT2, ZDT3, ZDT4 and SCH, which implies that the solutions obtained by NSGA-II has the wider distribution range. MOEA/D performs well on functions ZDT6 and FON. The distribution width metric DW of MOPSO is still the worst among the four algorithms.
In order to show the ability of convergence to Pareto front more clearly, the comparison results of four algorithms for any one run of eight test problems are shown in Figures 5-12.          It can be seen from Figure 5 that the non-dominated solutions generated by IMOPSO are superior to the other three algorithms in both convergence and distribution. The solutions obtained by the other three algorithms can converge to the true PF, but the distribution of the solutions is poor, especially MOPSO and MOEA/D. It can be seen from the figure that many solutions cannot be searched, resulting in the solutions generated by the algorithms not covering the entire Pareto front. The performance of NSGA-II is second only to that of IMOPSO. Figure 6 shows that the non-dominated solutions generated by the four algorithms can converge to the vicinity of the true PF. However, we can also find that the solution set generated by MOPSO and NSGA-II is not uniform enough and does not cover the entire Pareto front, which is also the reason for the poor statistical value of its distribution index. The distribution of the solution set obtained by MOEA/D is second only to that of IMOPSO. The solution set generated by IMOPSO has the best convergence and distribution among the four algorithms.
We can intuitively see from Figure 7 that the effect of MOPSO in solving KUR function is very poor. The non-dominated solutions obtained by MOPSO fail to converge to the true PF, and the solution set is not evenly distributed. The non-dominated solutions generated by NSGA-II and MOEA/D converge to the true PF. However, the solutions obtained by MOEA/D fail to cover the entire Pareto front, and the distribution of solutions is not as good as that of NSGA-II. IMOPSO is the best of the four algorithms.
It can be seen from Figure 8 that the solutions obtained by IMOPSO and MOPSO converge to the true PF. The solution set generated by MOPSO is not evenly distributed and its distribution width is poor, and the non-dominated solutions solved by NSGA-II also exists the problem of uneven distribution. The convergence of MOEA/D is poor, and many solutions deviate from the true PF. IMOPSO is superior to the other three algorithms in both convergence and distribution. Figure 9 indicates that the performance of the three metrics of the non-dominated solutions set generated by IMOPSO is the best. The non-dominated solution set generated by MOPSO can converge to the true PF, but the solutions distribution is poor. Compared with IMOPSO, NSGA-II also has the weakness of poor solutions distribution. The nondominated solutions obtained by MOEA/D deviate from the true PF, especially in the later stage. Figure 10 shows that the solution set generated by IMOPSO is superior to the other three algorithms in both convergence and distribution. The non-dominated solutions obtained by MOPSO can converge to the true PF, but the distribution of the solutions is uneven. There are cases where MOPSO fails to find the non-dominant solutions. The solution set generated by NSGA-II deviates slightly from the true PF, but it can be seen from the distribution of solutions that NSGA-II is better than MOPSO and MOEA/D. MOEA/D fails to search the complete Pareto front effectively, which is one of the reasons why the solution set generated by MOEA/D is poorly distributed. Figure 11 indicates the distribution of solution sets generated by the four algorithms. IMOPSO is the best among the four algorithms in terms of convergence and distribution. The non-dominated solutions generated by MOPSO is far away from the true PF, and its performance in solving this problem is very poor. Non-dominated solutions generated by NSGA-II and MOEA/D can converge to the true PF, and the distribution characteristics of solutions obtained by MOEA/D are better than those of NSGA-II, which is mainly manifested in the relatively uniform distribution of solutions.
We can intuitively see from Figure 12 that both IMOPSO and MOEA/D have obtained solutions with excellent convergence and distribution, but IMOPSO performs better than MOEA/D. Non-dominated solutions solved by NSGA-II all converge to the true PF, but the distribution of solutions is not as uniform as MOEA/D. The solutions generated by MOPSO also converge to the true PF, but the distribution of the solution is poor, and many solutions are not found.

Simulation Results of the DWTA Problem
For three DWTA problems of different scales (small, medium and large), four algorithms, i.e., IMOPSO, MOPSO, NSGA-II and MOEA/D, are used to solve them respectively. Each algorithm runs independently for 30 times. The experimental results are shown in the Table 5. It can be seen from the results in Table 5 that the IMOPSO algorithm has achieved the best IGD values in all three instances, and it indicates that IMOPSO has better performance than the other three algorithms in solving multiobjective DWTA problems. With the increase of the scale of DWTA problems, the non-dominated solutions generated by IMOPSO also obtain excellent convergence and distribution. The performance of NSGA-II in solving multiobjective DWTA issues is second only to that of IMOPSO, especially in medium scale, and the statistical values of NSGA-II are very close to that of IMOPSO. MOPSO has little difference with NSGA-II in small-scale and large-scale instance, and its overall performance is second only to NSGA-II. The effect of MOEA/D in solving DWTA problems is very poor. Its index values are far behind the other three algorithms, and this gap is becoming more and more obvious with the increase of the problem scale. This is due to the strictness of the five constraints in the multiobjective DWTA problem discussed in this paper, which leads to the increased difficulty of solving the model. Moreover, MOEA/D needs to normalize the objectives when dealing with multiobjective DWTA problems, so the computational complexity is greatly increased.
To visually compare the abilities of the four algorithms in solving multiobjective DWTA problems of three scales, we give the comparison results between the solution sets generated by any one run of the four algorithms at each scale and the true PF, and the comparison results are shown in Figures 13-18.      It can be seen by Figures 13 and 14 that IMOPSO performs very well in solving small-scale multiobjective DWTA problem, and the solutions obtained by IMOPSO can converge to the true PF. The solutions obtained by MOPSO failed to converge to the true PF, and the distribution of the solution set is poor, which failed to completely cover the entire Pareto front. NSGA-II is second only to IMOPSO, and its non-dominated solutions mostly converge to Pareto front. MOEA/D is the worst one, which almost obtain no useful Pareto solutions. Figures 15 and 16 show the distribution comparison of solution sets of four algorithms in medium scale. It can be seen that most of the solutions obtained by IMOPSO converge to the true PF, and only a few solutions fail to fall on the Pareto front. The non-dominated solutions generated by NSGA-II have a distance from true PF, and its convergence is far less than that of IMOPSO. Compared with NSGA-II, the solution set obtained by MOPSO is farther from the true PF. MOEA/D is the worst among the four algorithms, which generates very few solutions, and the solution set cannot completely cover the entire Pareto front. Figures 17 and 18 are the comparison results of four algorithms for solving large-scale multiobjective DWTA problem. It can be seen that the set of non-dominated solutions obtained by IMOPSO for solving large-scale multiobjective DWTA problem can fully converge to the Pareto front, and most of the non-dominated solutions have the same values as the solutions of true PF. Some solutions generated by NSGA-II can converge to the true PF, but there is still a distance from the Pareto front in the middle location. The solution set obtained by MOPSO fails to converge to the true PF, but the distribution of solutions is relatively good. MOEA/D only finds three solutions. It can be seen that with the increase of DWTA problem scale, MOEA/D is getting worse and worse.

Discussion
In order to evaluate the improved multiobjective particle swarm optimization algorithm proposed in this paper, we conducted two sets of comparative experiments with three multiobjective optimization algorithms, i.e., IMOPSO, NSGA-II and MOEA/D. According to the experimental results of benchmark functions, the IMOPSO algorithm proposed in this paper has achieved the best values in all three evaluation indexes, which shows that the solutions obtained by IMOPSO in solving MOPs have better convergence and distribution. Through the solution set distribution graphs, we can also intuitively see the superiority of IMOPSO in solving MOPs compared with the other three algorithms. The results in the figure just correspond to the statistical data of the evaluation index, which shows the rationality of the data. Comparative experiments of multiobjective DWTA problems show that IMOPSO algorithm can also solve satisfactory solutions with the increasing scale of DWTA problems, thus obtaining an effective weapon-target allocation scheme. The results of solving the multiobjective DWTA problem reflect the ability of IMOPSO in dealing with practical MOPs, which shows that IMOPSO algorithm is more suitable for solving large-scale mixed integer programming problems than other algorithms.
We use the average value of each performance metric in the experiment to measure the improvement degree of IMOPSO compared with other algorithms of dealing with MOPs. The specific method is as follows: We sum the metric average values of each benchmark function, and then find the average value as the comprehensive index value. The DWTA problem also adopts the same processing method. Through the calculation, we can conclude that IMOPSO is superior to MOPSO, NSGA-II and MOEA/D in GD/∆/DW/IGD index, which are 99.61%/74.46%/99.68%/42.99%, 29.97%/39.90%/87.67%/22.50%, 35.33%/ 57.22%/96.96%/93.04%, respectively. It can be seen from the above values that the improved effect of the proposed IMOPSO algorithm on the basis of MOPSO algorithm is particularly remarkable, and the performance of IMOPSO in solving MOPs is also significantly ahead of NSGA-II and MOEA/D. IMOPSO can be regarded as the best of the four algorithms, which successfully overcomes the shortcoming that MOPSO is easy to fall into local optimum, and can also balance the contradiction between exploration and exploitation. Experimental results show that compared with other algorithms, IMOPSO is not only more robust and accurate, but also has a wide range of applications, which is more suitable for dealing with complex MOPs.
Although the IMOPSO algorithm has many excellent performances, it is not completely flawless when dealing with MOPs. At first, the number of times that the IMOPSO algorithm needs to call the objective function in each iteration is very large. Since the benchmark functions have simple and clear mathematical expressions, it is not an insurmountable problem for the optimization system itself to call the objective function more times. However, many problems that need to be solved in engineering application are "black box" types, which cannot be directly expressed by mathematical formulas. This requires huge computer resources as support, so it is necessary to properly handle them in the optimization program. Secondly, in this paper, it is difficult to determine the specific number of unimproved generations when the g best proceeds mutation operation, so we conducted many experiments to determine it.

Conclusions and Future Works
To solve the multiobjective DWTA problem effectively, a novel MOPSO algorithm (IMOPSO) is proposed in this paper. The algorithm constructs different learning strategies for dominated and non-dominated solutions in the current population, which makes the learning pertinence of particles clearer. IMOPSO adopts a search strategy based on SBX and PM. The exchange of information among particles is enhanced by the SBX operator of p best in external archive, and the PM operation is performed on newly generated g best , thereby generating new global optimal particles. SBX and PM address the drawback that MOPSO algorithm is easy to fall into local optimal solution. The algorithm also uses a dynamic external archive maintenance strategy, which improves the diversity of non-dominated solutions. Two sets of experiments based on benchmark functions and multiobjective DWTA problem are used to verify the performance of IMOPSO. In order to illustrate the effectiveness of IMOPSO in solving MOPs, we compared three state-of-the-art multiobjective optimization algorithms: MOPSO, NSGA-II and MOEA/D. The experimental results show that the non-dominated solution obtained by the proposed IMOPSO algorithm is superior to other algorithms in convergence and distribution. Furthermore, IMOPSO has advantages in solving large-scale mixed integer programming problems.
The benchmark functions and multiobjective DWTA model selected in this paper only have two objectives. In reality, there are often more optimization objectives in optimization problems. In the future, we will use more types of benchmark functions which include more objectives to verify the ability of the IMOPSO algorithm. In addition, IMOPSO will be used to solve more multiobjective optimization problem models. Moreover, we will compare IMOPSO with more intelligent algorithms. At present, there are many solvers for multiobjective optimization models, and we will use other solvers (e.g., Gurobi and CPLEX) to solve DWTA problems in future research.