Multi-weapon multi-target assignment based on hybrid genetic algorithm in uncertain environment

Abstract

The multi-weapon multi-target assignment is always an unavoidable problem in military field. It does make sense to find a proper assignment of weapons to targets which may help maximize the attack effect. In this article, as the information achieved from the battlefield is becoming more and more uncertain, a novel threat assessment method and target assignment algorithm are proposed against the background of unmanned aerial vehicles intelligent air combat. Specifically, with regard to the threat assessment issue, a possibility degree function based on grey theory is structured to further improve the grey analytic hierarchy process. It can transform the interval weight of threat factors into scalar-valued weight, with which the accuracy of threat assessment can be improved. Regarding the target assignment problem, combining with interval grey number, an improved hybrid genetic algorithm is developed. The improvements are mainly consisting of adaptive crossover and mutation operators which can help to find an approximate solution within certain time constraints. Meanwhile, the simulated annealing operation is incorporated to avoid local optimum and premature phenomenon. In addition, the selection operation and fitness function are also redesigned to handle the interval numbers. Simulation results demonstrate the effectiveness of our algorithm in completing the multi-objective weapon-target assignment under uncertain environment.

Introduction

With the progress in the information technology, modern air combat technology is developing rapidly in the track of autonomous and intelligent combat. Throughout this process, the weapon-target assignment (WTA) problem has always been a fundamental problem in the battlefield decision for firepower strike. A general WTA problem attempts to allocate and schedule various weapon resources with the objective of minimizing surviving target value. It should take into account the real-time states of the warring parties, including the damage ability of our weapon systems against hostile threats and the different economic value of each weapon. Essentially, the WTA problem is a typical combinatorial optimization problem with complex constraints and has been proved to be NP-complete. ¹ In contrast to traditional WTA problem, multi-weapon multi-target assignment (MWMTA) is a higher dimensional problem involving more criterions. It is obviously more complicated and more in line with the real operational decision-making. In previous research, academies endeavor to solve WTA problem from different aspects theoretically by proposing several effective algorithms, ² such as particle swarm algorithm ³ and improved ant colony algorithms. ^4,5 In addition, Zaman et al. ⁶ proposed a heuristic genetic approach to address such optimization problem, and simulated annealing genetic algorithm (SAGA) is also used as a solution to it in the study by Tokgöz and Bulkan. ⁷

Nowadays, WTA has been applied in a variety of combat modes. Among these combat modes, air to air combat undoubtedly plays the most decisive role, especially in the advent of the unmanned aerial vehicle (UAV) era. UAV can make autonomous intelligent decision quickly and accomplish cooperative air combat through aircraft formation. In the context of beyond-visual-range air combat, ⁸ MWMTA has become one of the key technologies for unmanned combat aircraft ⁹ fire control system. Thus, this article takes UAVs as an example platform and unites the reality of UAVs air combat to analyze the MWMTA problem. Besides, another two problems also attract considerable attention: threat assessment and uncertain environment.

Threat assessment is an estimation of the lethality or threat of the enemy to our forces. For multi-UAV coordinated air combat, based on the situation assessment, threat assessment integrates the prior knowledge of enemy’s destructive capability, maneuverability, motion patterns, and operational intention to evaluate the severity of an engagement. It can be regarded as a prerequisite for WTA, and a number of methods have been put forward to solve this problem. For example, TOPSIS, ¹⁰ analytic hierarchy process (AHP), ¹¹ and so on have focused on the matter of battlefield environment assessment. There is one more point that needs to be stated in advance. For the sake of convenience, the process composed of both threat assessment and WTA, hereafter this article will be defined as the “decision-making.”

Although the above methods can deal with such optimization problems, they did not consider the actual situation filled with uncertainty ¹² no matter in the threat assessment or the allocation stage. The uncertainty discussed here can be divided into two aspects. On the one hand, due to the increasing interference factors and the error of sensors, the information obtained from the air combat process must have considerable uncertainties. On the other hand, the stealth of fighters is constantly improved, and we may not grasp all the performance parameters of enemy aircraft. Therefore, during the evaluation process, some information is incomplete and it is difficult to characterize them with an accurate value. In this study, it is essential to find ways to describe and deal with the flexible information. The traditional intelligent algorithm also needs to be improved so that it can apply to the actual situation and overcome its own drawbacks.

This article mainly focuses on the issue of threat assessment and MWMTA with air combat as the background. Nevertheless, with the continuous improvement of aircraft stealth, the increase of interference and the measurement error of the sensor itself, the uncertainty of the information acquired also keeps increasing. The main contributions in this article are summarized as follows.

Grey number-based threat assessment: For threat assessment, the interval grey number-based grey analytic hierarchy process (GAHP), which could characterize uncertain and incomplete information in air combat, is developed to comprehensively evaluate the threat of the target. Although the fuzzy set ^13,14 and interval analytic hierarchy process (IAHP) ¹⁵ can handle ambiguous information, they would lead to high computational complexity of the whole decision-making process. For this part, we propose a possibility degree function to quantize the interval grey number, which basically convert the grey interval weights into scalar-valued weights. This operation will greatly reduce the complexity of all the assignment algorithms which are not compatible with the interval threat values.

Interval threat value compatible solution: We combine the uncertainty into the target assignment problem, while the conventional intelligent optimization algorithm often requires a certain value. For the WTA problem, although there have been artificial bee colony algorithm, ¹⁶ ant colony optimization, ¹⁷ particle swarm optimization, ¹⁸ and genetic algorithm (GA) with greedy eugenics, ¹⁹ they fail to be applied in air combat WTA under uncertain environments and do not apply to the case where the solution of the fitness function is in interval form. To solve this problem, the possibility degree function is combined with the GA and used as the basis for selection operation. We design an appropriate fitness function and selection mechanism, thereby making the algorithm become suitable for the optimization whose individual fitness value is grey interval.

A novel hybrid GA-based target assignment: For the target assignment algorithm, to get rid of the shortcomings of the algorithm, for example, slow convergence speed and challenge to achieve optimal solution, a novel hybrid GA is presented in this article. The adaptive crossover and mutation operators are introduced to help find an approximate solution within certain time constraints. Meanwhile, the simulated annealing (SA) operation is incorporated to avoid premature phenomenon. Especially, since common adaptive genetic algorithm (AGA) ²⁰ also has some drawbacks such as poor stability and easy to fall into local optimum, a new modified adaptive genetic operator is presented to address these problems. Compared with the traditional GA ^21,22 or the genetic approach only with SA, ^23,24 the proposed method does have better performance.

The rest of this article is organized as follows. In the second section, the air combat threat assessment and MWMTA problems are described. In the third section, improved GAHP methods are developed to determine the interval weights of threat factors and the possibility degree function is introduced. In the fourth section, a target assignment model is established. In the fifth section, an improved AGA which involves SA mechanism is presented. In the sixth section, simulations are conducted to verify the effectiveness of our proposed algorithm. Finally, conclusions are given in the seventh section.

Problem description

Before formulating the air combat decision-making problem under uncertain environment, the relations between the various parts of the problem are introduced. Firstly, the synthetic threat assessment values of the enemy aircraft against ours are obtained through the proposed threat evaluation method. Only when the threat values are obtained can we allocate the weapons in the next step. Secondly, a target assignment model is set up. Thirdly, an optimization algorithm is introduced to solve the previous model. Finally, according to the threat evaluation values and the target assignment model, the MWMTA problem is simulated with the proposed intelligent algorithm. During this process, some following key issues should be addressed.

Air combat threat index assessment model

Since the multi-objective cooperative air combat will absolutely be the trend in the future, the command system can evaluate and rank threats through battlefield information, and then fully deploy the aircraft to develop the best cooperative attack or defense scheme. The air combat situation has a great impact on threat assessment.

The relationship between both sides in air combat is shown in Figure 1, here v_R and v_B , respectively, denote the velocity vectors of our fighter and the target, q_R is the off-boresight angle of our fighter relative to the target, q_B is the target approaching angle, and r represents the relative distance between our aircraft and the target. Then, the threat factors of angle, velocity, and distance can be expressed by

Figure 1.

Sketch of air combat situation.

• 1. Angle threat factor $T_{\alpha }$

$T_{\alpha }=\frac{\left|q_B\right|+\left|q_R\right|}{{360}^{\circ }}$

1

• 2. Velocity threat factor T_v

$T_v=\left\{\begin{array}{cc}0.1& v_B<0.6v_R\\ -0.5+\frac{v_B}{v_R}& 0.6v_R\le v_B\le 1.5v_R\\ 1.0& v_B>1.5v_R\end{array}\right.$

2

• 3. Distance threat factor T_r

$\begin{array}{cc}{T}_r=\hfill & \left\{\begin{array}{cc}0.5\hfill & r\le r_{mr},r\le r_{mt}\hfill \\ 0.5-0.2\left(\frac{r-r_{mt}}{r_{mr}-r_{mt}}\right)\hfill & r_{mt}<r<r_{mr}\hfill \\ 1.0\hfill & r_{mr}<r<r_{mt}\hfill \\ 0.8\hfill & \text{max}\left(r_{mr},r_{mt}\right)<r<r_{dr}\hfill \end{array}\right.\hfill \end{array}$

3

where $r_{mr}$ is the maximum range of our aircraft missiles, $r_{mt}$ is the target maximum attack range, $r_{dr}$ is the maximum tracking distance of our radar.

In addition to the above threat factors, altitude and stealth capability also affect the threat assessment. Besides, when the aircraft is assembled and put into the battle, its maneuverability and weapon performance can be considered to remain unchanged, accordingly the targets air combat effectiveness can be taken into account as a static threat factor.

Based on different air combat environments and modes, the threat assessment uses established models to evaluate various types of threat factors and calculates the weights of them according to the relative situation information on the battlefield. The threat assessment array makes a comparative analysis on the threaten degree of individuals participating in the war and establishes the final operational scheme.

MWMTA problem in cooperative air combat

Coordinated target assignment is another important part of the decision-making problem. Before we establish a mathematical model, some assumptions should be defined at first.

Assumption 1

This air battle is composed of m fighters and n enemy targets while our aircraft carry r missiles. Each fighter in the multi-UAV swarm is equivalent to an independent unit, which possesses different quantities of missiles.

Assumption 2

Each fighter can use any number of missiles it carries to attack a target. Every missile must be assigned to targets and each missile can only attack one target.

Assumption 3

The kill probability between the kth missile and the jth target is already known and is labeled as $p_{kj}$ . The threat value of the jth target against our ith aircraft is also calculated beforehand through the threat assessment.

The air combat target assignment is to cause more damage to the targets on the premise of preserving ourselves as much as possible. Let $x_{kj}=1$ denote that missile k is allocated to attack target j, and $x_{kj}=0$ denote that missile k does not attack target j. If the missile k is assigned to the target j, the survival probability of target j is $(1-p_{kj})$ . After performing a coordinated attack, survival probability of target j becomes ${\displaystyle \prod _{k=1}^r}{(1-p_{kj})}^{x_{kj}}$ . Moreover, let $T_{ji}$ be the threat value of the jth target to our ith aircraft, the remaining threat ^25,26 can be written as $T_{ji}\times {\displaystyle \prod _{k=1}^r}{(1-p_{kj})}^{x_{kj}}$ .

Accordingly, when we just consider the maximum expected kill value, the function can be derived as

$E\left(x\right)={\displaystyle \sum _{j=1}^n}\left[1-{\displaystyle \prod _{k=1}^r}{\left(1-p_{kj}\right)}^{x_{kj}}\right]$

4

and when reducing casualties is taken into account separately, the remaining threat of the targets after a round of engagement can be expressed by

$F\left(x\right)={\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{i=1}^m}\left[T_{ji}\times {\displaystyle \prod _{k=1}^r}{\left(1-p_{kj}\right)}^{x_{kj}}\right]$

5

Uncertainty in decision-making problems

Although the threat assessment and target assignment have been described separately, among them, it is not clear which problems caused by uncertainty need to be resolved.

Description of uncertain information in air combat

The command system requires the assistance of airborne equipment for battlefield environment analysis and information management to get a reasonable allocation scheme. However, with the continuous improvement of aircraft stealth, the increase of interference and measurement error of the airborne equipment, the uncertainties of information acquired from the scene also keeps increasing. How to accurately describe this uncertain information has important practical significance.

The rapidity and accuracy of decision-making in uncertain environment

Due to the rapid changes in air combat, the command system must make correct judgments within a very short time. It is of great significance to study how to accurately and timely distribute the targets to be attacked according to the current situation. In order to improve the overall combat capability and survivability of the air fleet, we can start from two aspects. First, since the information we collect is of great uncertainty and incompleteness, the data derived through a series of calculation are unlikely to be given in the form of real number. It does make sense to develop a means which can convert interval numbers into real numbers, thus reducing the computational complexity and saving computation time. Moreover, as traditional intelligent algorithm might have some shortcomings, we have to modify the algorithm so as to guarantee the rapidity and accuracy of weapon-target allocation under uncertain environment.

The connection between uncertain information and target assignment algorithm

As for the intelligent optimization algorithm, besides the drawbacks of the algorithm itself, its optimization process often requires the data to be determined. This article aims to study the issue of target assignment under uncertain environment and needs to deal with flexible information in the instantaneous air combat situation. Since common optimization algorithms are unable to cope with uncertain information, how to establish a connection between uncertainties and the algorithm is another challenge we face.

Interval grey number-based air combat threat assessment

In actual air combat, the performance parameters of enemy aircraft cannot be clearly displayed. This section mainly studies the target threat assessment method under uncertain environment. Based on the modified AHP, ²⁷ we introduce the concept of interval grey number and possibility degree function. Thus, the uncertainty of information can be accurately represented and the rationality of threat assessment results can be guaranteed.

Related conception of interval grey number and GAHP

To help understand, some definitions of interval grey number from grey system theory ²⁸ should be introduced at first.

Definition 1

A number that only knows the approximate range while does not know its exact value is called a grey number. It is the element or cell that makes up the grey system and is usually denoted by the symbol $\otimes$ .

Definition 2

The grey number with both lower bound $\stackrel{\_}{a}$ and upper bound $\overline{a}$ is called the interval grey number, denoted as $\otimes \in [\stackrel{\_}{a},\overline{a}]$ . And when the lower and upper bounds are infinite, it is called the black number; when $\otimes \in [\stackrel{\_}{a},\overline{a}]$ and $\stackrel{\_}{a}=\overline{a}$ , we call $\otimes$ the white number. In addition, $\otimes (\gamma )={\otimes }^{-}+({\otimes }^{-}-{\otimes }^{+})\gamma$ , $(0\le \gamma \le 1)$ is the standardized form of interval grey number.

Definition 3

The basic arithmetic of interval grey number is as follows.

Suppose that ${\otimes }_1\in [a,b]$ , $a<b$ , ${\otimes }_2\in [c,d]$ , $c<d$ .

The addition operation are defined as ${\otimes }_1+{\otimes }_2\in [a+c,b+d]$ .

Subtraction ${\otimes }_1-{\otimes }_2\in [a-d,b-c]$ .

Multiplication ${\otimes }_1\cdot {\otimes }_2\in [\text{min}\{ac,ad,bc,bd\},\text{max}\{ac,ad,bc,bd\}]$ .

Division is defined as $\frac{{\otimes }_1}{{\otimes }_2}\in [\text{min}\{\frac{a}{c},\frac{a}{c},\frac{b}{c},\frac{b}{d}\},\text{max}\{\frac{a}{c},\frac{a}{c},\frac{b}{c},\frac{b}{d}\}]$ , where $c\ne 0$ , $d\ne 0$ , $cd>0$ .

GAHP is the product of the mix of grey theory and AHP. In cybernetics, people often use the depth of color to describe the degree of clarity of information, and the grey indicates whether the information is clear or not. Due to the complexity of the air combat environment, there are often uncertainties in small samples and poor information. We need to adopt GAHP to evaluate the target threat. On the basis of determining the weights of threat indicators by IAHP, the grey theory is used to convert the decentralized information from experts into weight vectors of different grey degrees. Then the comprehensive evaluation value can be obtained through scalarization.

Compared with the traditional AHP, GAHP can deal with incomplete information under the interaction of various factors. In addition, compared with the IAHP, GAHP tends to study objects with clear extension and unclear connotation.

Procedure of the GAHP

• 1. Set up hierarchical structure of threat assessment indicators

First, the object to be evaluated is decomposed into several layers, and each layer consists of several elements. Among them, the top layer is the target layer, which represents the ultimate goal of the entire system. The bottom layer is the solution layer, which provides various solutions for implementing the target layer. This top-down dominance creates a ladder hierarchy.

• 2. Determine the weight of evaluation indicators

Since each evaluation index has different importance to the evaluation target, it is necessary to determine the weight of each indicator first. Here, not only factors like angle, speed, distance, altitude, and air combat efficiency exponent are considered but also the strong stealth and super maneuverability of advanced fourth-generation fighters. ^{29 –31} Due to the uncertainty of some information, the IAHP is adopted in this article. First, according to the judgment matrix A given by the decision maker, where $A\text{=}\left[\begin{array}{cccc}{\otimes }_{11}& {\otimes }_{12}& \cdot \cdot \cdot & {\otimes }_{1n}\\ {\otimes }_{21}& {\otimes }_{22}& \cdot \cdot \cdot & {\otimes }_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {\otimes }_{n1}& {\otimes }_{n2}& \cdot \cdot \cdot & {\otimes }_{nn}\end{array}\right]$ , the length of each grey interval ${\otimes }_{ij}$ can be obtained by calculating $d({\otimes }_{ij})=\overline{a_{ij}}-\underset{\_}{a_{ij}},i,j\in N$ . Then normalize the length of each interval by

${\lambda }_{ij}=\frac{\frac{1}{d({\otimes }_{ij})}}{{\displaystyle \sum _{j=i+1}^n}\frac{1}{d({\otimes }_{ij})}},i,j\in N$

6

After that, the weight vector of the interval grey number is calculated by

$w_i={\displaystyle \sum _{j=i+1}^n}{\otimes }_{ij}{\lambda }_{ij}{w}_j$

7

here $j=n-1,n-2,\cdot \cdot \cdot ,1$ and $w_n>0$ . Thus the weight value of each index can be written as

$W=(w_1,w_2,\cdots ,w_6{)}^T$

8

where w ₁ is the weight of fighter combat capability, w ₂ is the weight of fighter angle threat, w ₃ is the weight of fighter distance threat,w ₄ is the weight of fighter velocity threat, w ₅ is the weight of fighter altitude threat, and w ₆ is the weight of fighter stealth threat. Finally, through normalization, we can get the single interval weight

$w_i=\frac{w_i}{{\displaystyle \sum _{i=1}^n}{w}_i}$

9

• 3. Determine evaluation sample matrix

The threat level of enemy fighter can be divided into five levels: very high, high, normal, low, very low, with their values are 9, 7, 5, 3, 1, respectively. Each indicator is graded by different specialists and given different weights. Then the weighted average method is used to obtain the following evaluation sample matrix

$D_{JI}{\hspace{0.17em}}^A=\left[\begin{array}{cccc}{d}_{11}{\hspace{0.17em}}^{(A)}& d_{12}{\hspace{0.17em}}^{(A)}& \cdots & d_{1i}{\hspace{0.17em}}^{(A)}\\ d_{21}{\hspace{0.17em}}^{(A)}& d_{22}{\hspace{0.17em}}^{(A)}& \cdots & d_{2i}{\hspace{0.17em}}^{(A)}\\ ⋮& ⋮& \ddots & ⋮\\ d_{j1}{\hspace{0.17em}}^{(A)}& d_{j2}{\hspace{0.17em}}^{(A)}& \cdots & d_{ji}{\hspace{0.17em}}^{(A)}\end{array}\right]$

10

here the entry in the jth row and the ith column represents the evaluation value of expert i for indicator j.

• 4. Determine evaluation grey clustering

To determine evaluation grey clustering is to determine the rating, grey number, and definite weighted function of grey clustering. According to the grading criteria, five corresponding grey clusters are set to indicate very high, high, normal, low, very low with e = 1, 2, 3, 4, 5 as their numbers. We can get the concrete form of the definite weighted function according to Figure 2.

Figure 2.

Definite weighted function diagram.

• 5. Calculate grey evaluation coefficients

According to the evaluation sample matrix (10) and the definite weighted function f_k (d_ij ), the grey evaluation coefficient of target J relative to evaluation index A that belongs to grey clustering K can be calculated by

$n_{JK}{\hspace{0.17em}}^{(A)}={\displaystyle \sum _{i=1}^i}{f}_k(d_{JI}{\hspace{0.17em}}^{(A)})$

11

Afterward, sum the grey evaluation coefficients of target J's every evaluation grey clustering by

$n_J{\hspace{0.17em}}^{(A)}={\displaystyle \sum _{i=1}^i}{f}_k(n_{JK}{\hspace{0.17em}}^{(A)})$

12

• 6. Calculate the weight matrix of grey evaluation weight vector

Divide (11) by (12), and the grey evaluation weight of target J that belongs to evaluation grey clustering K is given by

$r_{JK}^{(A)}=\frac{n_{JK}^{(A)}}{n_J^{(A)}}$

13

where $K=1,2,\cdots ,k$ , $J=1,2,\cdots ,j$ , relative to the evaluation index A, a grey evaluation matrix of the target can be derived as

$R^A=\left[\begin{array}{cccc}{r}_{11}{\hspace{0.17em}}^{(A)}& r_{12}{\hspace{0.17em}}^{(A)}& \cdots & r_{1k}{\hspace{0.17em}}^{(A)}\\ r_{21}{\hspace{0.17em}}^{(A)}& r_{22}{\hspace{0.17em}}^{(A)}& \cdots & r_{2k}{\hspace{0.17em}}^{(A)}\\ ⋮& ⋮& \ddots & ⋮\\ r_{j1}{\hspace{0.17em}}^{(A)}& r_{j2}{\hspace{0.17em}}^{(A)}& \cdots & r_{jk}{\hspace{0.17em}}^{(A)}\end{array}\right]$

14

• 7. Comprehensive evaluation

The synthetic grey evaluation vector ³² of index A is calculated by $U^A=W{R}^A$ , we assign $D=(9,7,5,3,1)$ separately to each rating of grey clusters. Finally, the comprehensive evaluation value can be expressed as $S=U^A{D}^T$ .

Description of the possibility degree function

According to the evaluation method in the previous section, the weights of the threat factors are given in the form of grey intervals. It brings a lot of computational trouble to the final decision-making. The interval weights of threat factors are required to be converted into scalar-valued weights accordingly. The possibility degree can be used to measure how much an interval grey number is greater or smaller than another one. First, the definition of point possibility degree function is introduced.

Definition 4

Assume that a is an arbitrary real number, $\otimes \in [{\otimes }^{-},{\otimes }^{+}]$ is an interval grey number, the point possibility degree function can be defined by

$p\left(a<\otimes \right)=\left\{\begin{array}{ll}1,\hfill & a<{\otimes }^{-}\hfill \\ \frac{{\otimes }^{+}-a}{{\otimes }^{+}-{\otimes }^{-}},\hfill & {\otimes }^{-}\le a\le {\otimes }^{+}\hfill \\ 0,\hfill & a>{\otimes }^{+}\hfill \end{array}\right.$

15

when a is less than $\otimes$ .

Definition 5

Suppose ${\otimes }_1\in [{\otimes }_1{\hspace{0.17em}}^{-},{\otimes }_1{\hspace{0.17em}}^{+}]$ and ${\otimes }_2\in [{\otimes }_2{\hspace{0.17em}}^{-},{\otimes }_2{\hspace{0.17em}}^{+}]$ are two arbitrary interval grey numbers, $\otimes (\gamma )$ is the standardization of interval grey number. If ${\gamma }_1$ is given and ${\otimes }_1({\gamma }_1)$ is a fixed point, $p({\otimes }_1({\gamma }_1)<{\otimes }_2)$ would be the possibility degree function under the circumstance where ${\otimes }_1({\gamma }_1)$ is less than ${\otimes }_2$ . Therefore, the possibility degree function with ${\otimes }_1$ less than ${\otimes }_2$ can be expressed as

$p({\otimes }_1<{\otimes }_2)={\displaystyle \underset{0}{\overset{1}{\int }}}p({\otimes }_1({\gamma }_1)<{\otimes }_2)d{\gamma }_1$

16

which appears to be complementary, namely, $p({\otimes }_1\ge {\otimes }_2)=1-p({\otimes }_1<{\otimes }_2)$ .

According to the integral solution method given by the literature, ³³ the possibility degree when ${\otimes }_1<{\otimes }_2$ can be expressed as

$p\left({\otimes }_1<{\otimes }_2\right)=\left\{\begin{array}{ll}1,\hfill & {\otimes }_1^{+}\le {\otimes }_2^{-}\hfill \\ 1-\frac{{\left({\otimes }_1^{+}-{\otimes }_2^{-}\right)}^2}{2\left({\otimes }_1^{+}-{\otimes }_1^{-}\right)\left({\otimes }_2^{+}-{\otimes }_2^{-}\right)},\hfill & {\otimes }_1^{-}\le {\otimes }_2^{-}\le {\otimes }_1^{+}\le {\otimes }_2^{+}\hfill \\ \frac{{\otimes }_2^{+}+{\otimes }_2^{-}-2{\otimes }_1^{-}}{2\left({\otimes }_1^{+}-{\otimes }_1^{-}\right)},\hfill & {\otimes }_1^{-}\le {\otimes }_2^{-}\le {\otimes }_2^{+}\le {\otimes }_1^{+}\hfill \\ \frac{{\left({\otimes }_1^{+}-{\otimes }_2^{-}\right)}^2}{2\left({\otimes }_1^{+}-{\otimes }_1^{-}\right)\left({\otimes }_2^{+}-{\otimes }_2^{-}\right)},\hfill & {\otimes }_2^{-}\le {\otimes }_1^{-}\le {\otimes }_2^{+}\le {\otimes }_1^{+}\hfill \\ 0\hfill & {\otimes }_2^{+}\le {\otimes }_1^{-}\hfill \\ 1-\frac{{\otimes }_2^{+}+{\otimes }_2^{-}-2{\otimes }_1^{-}}{2\left({\otimes }_1^{+}-{\otimes }_1^{-}\right)},\hfill & {\otimes }_2^{-}\le {\otimes }_1^{-}\le {\otimes }_1^{+}\le {\otimes }_2^{+}\hfill \end{array}\right.$

17

Furthermore, according to the competition mechanism of a league, the values of N grey intervals can be obtained after the accumulation of each round of comparison. It is also the theoretical basis for the final sorting of multiple interval grey numbers. Then each round of the comparison process can be described by the following dominant accumulation matrix shown in Table 1. In which ${\otimes }_1,{\otimes }_2,\cdots ,{\otimes }_N$ are N grey intervals and the entry in the ith row and the jth column represents the possibility degree with grey interval ${\otimes }_i$ less than ${\otimes }_j$ .

Table 1.

Weight evaluation matrix.

COMP	${\otimes }_1$	${\otimes }_2$	$\cdots$	${\otimes }_N$
${\otimes }_1$	$CO{M}_{11}$	$CO{M}_{12}$	$\cdots$	$CO{M}_{1N}$
${\otimes }_2$	$CO{M}_{21}$	$CO{M}_{22}$	$\cdots$	$CO{M}_{2N}$
		$\cdots$
${\otimes }_N$	$CO{M}_{N1}$	$CO{M}_{N2}$	$\cdots$	$CO{M}_{NN}$

Obviously, by defining the possibility degree function, the extent of a certain interval relative to another one can be described. Hence, this function is applicable to the comparison of multiple grey intervals.

Air combat MWMTA model

In this section, a target assignment model is set up to prepare for the next heuristic optimization algorithm.

According to the assessment of air combat threat in the last section, we can get the threat grey matrix $T={\left[{\otimes }_{ij}\right]}_{m\times n}$ of enemy aircraft against our aircraft at a certain moment. Assume that the number of aircraft which constitute our UAV formation is m, each aircraft is equipped with r_i missiles, and the quantity of aircraft in enemy formation is m. Numbering the missiles of the m aircraft, and the m aircraft should possess $r_1+r_2+\mathrm{...}+r_m=r$ missiles. Let $p_{kj}$ be the damage probability of the kth missile attacking the jth target, meanwhile, let $x_{kj}$ be the decision variable, so $x_{kj}=1$ denote the allocation of the kth missile to attack the jth target. According to (4) and (5), the mathematical model of air combat target decision is established as follows

$\left\{\begin{array}{l}\text{min}F={\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{i=1}^m}\left[T_{ji}\times {\displaystyle \prod _{k=1}^r}{(1-p_{kj})}^{x_{kj}}\right]\hfill \\ \text{max}E={\displaystyle \sum _{j=1}^n}\left[1-{\displaystyle \prod _{k=1}^r}{(1-p_{kj})}^{x_{kj}}\right]\hfill \end{array}\right.$

18

where F indicates that the minimum threat value of enemy residual targets, and E means destroying enemy targets as many as possible.

Constraints

$\text{s.t.}\left\{\begin{array}{l}{\displaystyle \sum _{k=1}^r}{x}_{kj}\le s_j, \text{j=1,2,}\mathrm{...}\text{, }n\hfill \\ {\displaystyle \sum _{k=1}^r}{\displaystyle \sum _{j=1}^n}{x}_{kj}\le r\hfill \\ {\displaystyle \sum _{j=1}^n}{x}_{kj}\le 1, \text{k=1,2,}\mathrm{...}\text{,}r\hfill \end{array}\right.$

19

here, constraint 1 denotes that up to s_j missiles can be assigned to attack the jth target, constraint 2 indicates that no more than r missiles can be used, and condition 3 denotes that one missile can only attack one target. Therefore, the simplified objective function ³⁴ is finally expressed by

$\begin{array}{ll}\text{min}\hfill & {\omega }_1\times {\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{i=1}^m}\left[T_{ji}\times {\displaystyle \prod _{k=1}^r}{(1-p_{kj})}^{x_{kj}}\right]\hfill \\ \hspace{0.17em}\hfill & + {\omega }_2\times (n-{\displaystyle \sum _{j=1}^n}\left[1-{\displaystyle \prod _{k=1}^r}{(1-p_{kj})}^{x_{kj}}\right])/n\hfill \\ \hspace{0.17em}\hfill & + N\times {\displaystyle \sum _{j=1}^n}{\left[\text{min}(0,s_j-{\displaystyle \sum _{k=1}^r}{x}_{kj})\right]}^2\hfill \end{array}$

20

where N is a penalty factor of the penalty function. It is mainly used to punish the occurrence of one missile attacking more than one target simultaneously in the solution space.

Improved hybrid GA for MWMTA problem

In this section, a new hybrid GA is proposed to solve the target assignment model and realize coordinated target assignment under uncertain environment. Here, improved genetic operators are designed in conjunction with SA to avoid local minima.

Conception of hybrid GA

GA simulates the selection, crossover, and mutation process of biological evolution. It gradually optimizes the search results until a global optimal solution is found or the preset number of iterations is reached. During the iterative process, this rule of survival of the fittest improves the population quality and overcomes the search randomness.

The hybrid GA adopted in this article is made up of two parts: AGA and SA algorithm. It combines the advantages of both, and in particular, adaptive genetic operators have also been redesigned to further improve the performance of the algorithm.

Suppose that the number of individuals is $NIND$ , the frequency of iterations is $MAXGEN$ , and the probability of crossover and mutation are p_c and p_m , respectively. In order to improve the search efficiency, $NIND$ initial solutions X ₀ that satisfy the constraint are generated. According to the above objective function model, the fitness value of the particle is calculated by $f_0=\left\{\begin{array}{cccc}f({\otimes }_{01})& f({\otimes }_{02})& \cdot \cdot \cdot & f({\otimes }_{0NIND})\end{array}\right\}$ .

Encoding operation

The chromosomes, to some degree, can reflect the nature and characteristics of an individual. The quality of chromosome encoding will directly affect the efficiency and search ability of the algorithm. Hence it is necessary for us to encode individual chromosomes based on the specific problem. In this section, each individual represents a missile-target distribution scheme X. As can be seen from the mathematical model of air combat target assignment in the previous section, our UAV formation has r missiles, which need to be allocated to m enemy aircraft. Besides, since $x_{kj}=1$ denotes that our kth missile is assigned to attack the jth target, the individual can be encoded as a binary matrix of $r\times n$ dimensions, which is composed of the decision variable $x_{kj}$ . In addition, there are some constraints on the variables of the target assignment model. For example, a missile can only attack one target. Therefore, in order to simplify the procedures and improve the search efficiency, NIND individuals which satisfy the constraints are generated randomly in the initial population. Finally, through such encoding means, the optimal solution we obtain is already a missile distribution scheme, and it does not even need to be decoded.

Selection

In this article, the possibility degree function is introduced as the basis of the selection. There is situation where the dominance of some individuals is very outstanding. However, if only one round of comparison was performed, for the fitness interval, the possibility degree of its dominant value could only be 1 at most, and the difference in the dominant value between individuals would not be too great. Hence we compare the grey interval ${\otimes }_i$ with all other grey intervals. Its cumulative dominant value can be calculated by summing the ith row or column of the cumulative matrix in Table 1. The cumulative dominant values of all grey intervals are derived in the same way. Then the index weight of each grey interval can be obtained through normalization processing. By doing so, individual fitness with greater dominance is more likely to be selected, and the algorithm can be qualified for the optimization whose individual fitness value is grey interval.

Adaptive genetic operators

The crossover probability p_c and the mutation probability p_m are the critical factors affecting the performance of the algorithm. Among them, p_c is the rate that affects individual population renewal, providing that the value of p_c is reasonable, the search rate of solution space could be improved, whereas an excessively large p_c will lead to premature convergence. Besides, the mutation operation is to maintain the diversity of the population and enhance the local search ability. In case p_m is too small, the population will lose some valuable genes; on the contrary, too large p_m makes the algorithm close to random search.

For different optimization problems, repeated experiments are needed to determine p_c and p_m . In the early stage of evolution, the population overall fitness value is low. It is necessary to use larger p_c and p_m to exchange information between different individuals. At the later stage of evolution, the chromosomes of individuals in the population tend to be stable. We prefer a smaller p_c , which can reduce the crossover operation and accelerate algorithm convergence.

Srinivas and Patnaik ²⁰ proposed an AGA with linear adaptive adjustment of crossover and mutation probabilities. However, common AGA also has problems such as poor stability and easy to fall into local optimal solution. Mixed with logistic function ³⁵ and cosine function, an improved adaptive genetic operator is presented in this article. The functions of genetic operator are expressed as follows

$P_c=\left\{\begin{array}{ll}{P}_{c\text{ max}}-\frac{P_{c\text{ max}}-P_{c\text{ min}}}{1+\text{exp}\left(A\text{cos}\left(\frac{f_{avg}-f^{\text{'}}}{f_{avg}-f_{\text{min}}}\pi \right)\right)}\hfill & f^{\text{'}}\le f_{\text{avg}}\hfill \\ P_{c\text{ max}}\hfill & f^{\text{'}}>f_{\text{avg}}\hfill \end{array}\right.$

21

$P_m=\left\{\begin{array}{ll}{P}_{m\text{ max}}-\frac{P_{m\text{ max}}-P_{m\text{ min}}}{1+\text{exp}\left(A\text{cos}\left(\frac{f_{avg}-f}{f_{avg}-f_{\text{min}}}\pi \right)\right)}\hfill & f\le f_{\text{avg}}\hfill \\ P_{m\text{ max}}\hfill & f>f_{\text{avg}}\hfill \end{array}\right.$

22

Here, $P_{c\text{ max}}$ is the maximum crossover probability, $P_{c\text{ min}}$ is the minimum crossover probability, $P_{m\text{ max}}$ is the maximum mutation probability, $P_{m\text{ min}}$ is the minimum mutation probability, A is a constant, $f_{\text{min}}$ is the minimum fitness value of the population, $f_{\text{avg}}$ is the average fitness value of the population, $f^{\text{'}}$ is the larger fitness value of two individuals involved in the crossover, and f represents the fitness value of the mutant individual.

In order to solve the minimum optimization problem, for individuals whose fitness is less than the average fitness value of the population, their crossover and mutation probabilities are reduced, so as to retain the excellent individuals and make the algorithm jump out of the local optimal solution. And if the individual fitness is higher than the average fitness value, the probabilities of crossover and mutation will be increased, which helps to maintain the diversity of the population and make the algorithm converge faster.

Simulated annealing

The SA starts from a high initial temperature. As the temperature decreases, the solution of the algorithm tends to be stable. However, the feasible solution might be a local optimal solution. The SA has to combine the probabilistic jump characteristics to find the global optimal solution.

In this article, the specific mechanism of SA is designed as follows. First, the initial temperature T ₀ and coefficient $\gamma \subset (0,1)$ are defined. When $T_0>T_{\text{min}}$ , that is, SA is not frozen, the optimal individual after mutation is denoted by P ₀ and the new individual is denoted by P ₁. Then, we calculate fitness values $Fit(P_0)$ and $Fit(P_1)$ . When $p(Fit(P_1)<Fit(P_0))>0.5$ , it can be considered that the new individual is better than the prior one. The individual in original population should be directly replaced by the new individual. Otherwise, when $p(Fit(P_1)<Fit(P_0))<0.5$ , it shows that the new individual is inferior to the prior one, and thus we should accept the inferior individual with a probability $p=\text{exp}(-\Delta F/T)$ where $\Delta F=Fit{(P_0)}^{+}-Fit{(P_1)}^{-}$ . At last, the current evolution temperature can be calculated by $T=T\text{×}\gamma$ . Compare the calculated T with the baseline temperature $T_{\text{min}}$ , if T is less than $T_{\text{min}}$ , go back to the beginning and continue the above operation. When the GA terminal condition is met, jump out of the loop.

Implementation of hybrid GA

From Figure 3, it clearly illustrates the implementation of the hybrid GA that can be outlined as follows.

Figure 3.

Hybrid genetic algorithm flow chart.

• Step 1: Coding

   Adopting binary coding, the solution is in matrix form. The number of our missiles r as well as the number of enemy aircraft n should be given at first. Then the distribution scheme can be expressed in the form of a binary matrix.

Step 2: Selection operation

   In the form of league system, the possibility degree of dominant value obtained by mutual comparison between individuals is accumulated.

Step 3: Adaptive crossover

   The crossover operation adopts single-point crossover. Before each crossover, according to (21), we can work out its probability p_c . By generating a random number r, if $r<p_c$ , a random cross point would be generated, and we can exchanged the substring between the tangent points to obtain a new offspring individual.

Step 4: Adaptive mutation

   Mutation is realized by simple substitution. Same as crossover, by generating a random number r, if $r<p_m$ , we have to reverse the original element.

Step 5: SA operation

   On the basis of the best chromosome which is found among the initial population after crossover and mutation, two rows and two columns are selected for exchange. The solution of the new individual is a feasible solution.

Simulation results and analysis

Target threat determination

In a air combat, take into account the performance of fighters, airborne equipment and weapons, as well as the characteristics of beyond visual range air combat and the coordination of aircraft, we establish an objective and comprehensive threat assessment model. Then, with the parameters given by database and experts, the grey matrix A for threat factor evaluation in Table 2 is selected.

Table 2.

Weight evaluation grey matrix A.

	$T_{c_j}$	$T_{{\alpha }_{ji}}$	$T_{r_{ji}}$	$T_{v_{ji}}$	$T_{h_{ji}}$	$T_{y_{ji}}$
$T_{c_j}^1$	[1.0, 1.0]	[0.2, 0.8]	[0.2, 0.3]	[0.1, 0.2]	[0.25, 0.4]	[0.2, 0.4]
$T_{{\alpha }_{ji}}$	[1.25, 5.0]	[1.0, 1.0]	[1.4, 1.6]	[0.8, 1.2]	[2.0, 2.0]	[1.0, 1.25]
$T_{r_{ji}}$	[3.3, 5.0]	[0.63, 0.71]	[1.0, 1.0]	[0.6, 0.8]	[1.33, 1.40]	[0.33, 0.50]
$T_{v_{ji}}$	[5.0, 10.0]	[0.83, 1.25]	[1.25, 1.67]	[1.0, 1.0]	[2.0, 2.4]	[1.0, 1.0]
$T_{h_{ji}}$	[2.4, 4.0]	[0.5, 0.5]	[0.71, 0.75]	[0.42, 0.50]	[1.0, 1.0]	[0.77, 0.80]
$T_{y_{ji}}$	[2.5, 5.0]	[0.8, 1.0]	[2.0, 3.0]	[1.0, 1.0]	[1.25, 1.3]	[1.0, 1.0]

Here $T_{c_j},T_{{\alpha }_{ji}},T_{r_{ji}},T_{v_{ji}},T_{h_{ji}},T_{y_{ji}}$ respectively represent the interval values of target air combat effectiveness, angle, distance, velocity, altitude, and stealth. And according to the method mentioned above, the interval weights of six threat factors can be obtained, $w_1=[0.18,0.21],w_2=[0.24,0.28],w_3=[0.45,0.53],w_4=[0.61,0.72],w_5=[0.32,0.44],w_6=[0.23,0.41]$ .

Thus the final dominant value of six grey intervals can be obtained by combining the possibility degree function with the interval weights of six threat factors we have calculated. The dominant accumulation matrix of the six grey interval weights derived through league competition is shown in Table 3. After doing a sum over matrix columns, the cumulative dominant value matrix of each grey interval can be calculated as $W=[\mathrm{0.5,1.5,3.57,5.5,4.31,2.62}]$ .

Table 3.

Dominant accumulation matrix under league system.

	${\otimes }_1$	${\otimes }_2$	${\otimes }_3$	${\otimes }_4$	${\otimes }_5$	${\otimes }_6$
${\otimes }_1$	0.5	1	1	1	1	1
${\otimes }_2$	0	0.5	1	1	1	1
${\otimes }_3$	0	0	0.5	1	0.81	0.12
${\otimes }_4$	0	0	0	0.5	0	0
${\otimes }_5$	0	0	0.19	1	0.5	0
${\otimes }_6$	0	0	0.88	1	1	0.5

Then combined with the indexes given by experts, we can evaluate the sample matrix

$D=\left[\begin{array}{cccccc}8& 7& 6.5& 6& 8& 9\\ 7& 8& 6.5& 6& 5.5& 6\\ 6.5& 7& 6.5& 8& 7& 6\\ 7& 6& 5.5& 6& 8& 6\\ 8& 7& 5.5& 6.5& 6& 7\\ 6& 7& 7& 8& 6.5& 7\end{array}\right]$

According to (11), the grey evaluation coefficient can be calculated by

$\begin{array}{l}E=1, x_{11}=f_1(8)+f_1(7)+f_1(6.5)+f_1(6)+f_1(8)+f_1(9)=4.944\\ E=2, x_{12}=f_2(8)+f_2(7)+f_2(6.5)+f_2(6)+f_2(8)+f_2(9)=4.3571\\ E=3, x_{13}=f_3(8)+f_3(7)+f_3(6.5)+f_3(6)+f_3(8)+f_3(9)=3.1\\ E=4, x_{14}=f_4(8)+f_4(7)+f_4(6.5)+f_4(6)+f_4(8)+f_4(9)=0\\ E=5, x_{15}=f_5(8)+f_5(7)+f_5(6.5)+f_5(6)+f_5(8)+f_5(9)=0\end{array}$

Afterward, through (12) and (13), the grey evaluation weight of each indicator can be determined by

$x_1=12.4015,r_1=(0.3987,0.3513,0.25,0,0)$

Similarly

$\begin{array}{l}{r}_2=(0.2711,0.3307,0.3878,0.0104,0)\\ r_3=(0.2714,0.3830,0.2264,0.1192,0)\\ r_4=(0.2338,0.4369,0.2769,0.0107,0)\\ r_5=(0.2927,0.4328,0.2635,0.0110,0)\\ r_6=(0.3144,0.4334,0.2522,0,0)\\ U^A=W{R}^A=(0.3153,0.3707,0.2865,0.0249,0)\\ S=U^A{D}^T=3.25\end{array}$

At last, associated with the grey cluster rating $D=\left(9,7,5,3,1\right)$ defined before, the comprehensive threat value of 3.25 is between low and normal, which lays a foundation for the next step of firepower distribution.

Simulation of target assignment

Scenario 1: Our UAV formation consists of four aircrafts while the enemy formation consists of six aircrafts. Each of our aircraft carries two missiles and is capable of attacking two enemy aircraft. Furthermore, to reduce the consumption of missiles, each enemy aircraft is attacked by at most two missiles. Let the target allocation scheme be a 48-dimensional discrete 0-1 space. A is set to 9.903438. The number of initial population individuals is set to 150, and the iteration times are set to 100. In addition, the weights ${\omega }_1$ and ${\omega }_2$ of the target decision mathematical model from the fourth section are set to 0.6 and 0.4, and the penalty factor N is set to 100. Then, the hit rate of our eight missiles against six enemy targets and the threat grey intervals of six enemy aircraft attacking our four targets are given in Tables 4 and 5.

Table 4.

Hit rate of our missiles against enemy targets.

	Hit rate
Missiles	1	2	3	4	5	6
1	0.29	0.92	0.23	0.89	0.14	0.72
2	0.49	0.82	0.41	0.10	0.51	0.37
3	0.33	0.46	0.39	0.90	0.43	0.30
4	0.12	0.38	0.52	0.26	0.71	0.41
5	0.22	0.15	0.44	0.29	0.21	0.13
6	0.61	0.95	0.76	0.32	0.24	0.19
7	0.44	0.56	0.16	0.22	0.88	0.17
8	0.81	0.42	0.35	0.43	0.77	0.90

Table 5.

Grey threat value matrix of enemy targets against us.

	Target threat
Targets	1	2	3	4
1	[0.50, 0.60]	[0.25, 0.28]	[0.19, 0.25]	[0.61, 0.61]
2	[0.41, 0.52]	[0.20, 0.29]	[0.34, 0.38]	[0.10, 0.14]
3	[0.33, 0.39]	[0.25, 0.25]	[0.23, 0.31]	[0.63, 0.71]
4	[0.28, 0.32]	[0.60, 0.63]	[0.37, 0.37]	[0.48, 0.51]
5	[0.58, 0.64]	[0.27, 0.30]	[0.16, 0.18]	[0.29, 0.29]
6	[0.25, 0.32]	[0.48, 0.54]	[0.33, 0.35]	[0.19, 0.19]

Comparing the simulation results of standard GA, SAGA, and the adaptive genetic simulated annealing algorithm, the curves of the comprehensive optimal fitness expectation, residual target threat value F, and the target killing expected value E are given in Figures 4 to 6.

Figure 4.

Comprehensive optimal fitness expectation (Scenario 1). (a) Standard genetic algorithm. (b) Simulated annealing genetic algorithm. (c) Improved hybrid genetic algorithm.

Figure 5.

Residual target threat value F (Scenario 1). (a) Standard genetic algorithm. (b) Simulated annealing genetic algorithm. (c) Improved hybrid genetic algorithm.

Figure 6.

Target killing expected value E (Scenario 1). (a) Standard genetic algorithm. (b) Simulated annealing genetic algorithm. (c) Improved hybrid genetic algorithm.

From Figures 4 to 6, we can see that:

Due to the introduction of the probability interval, the optimal fitness value F and the residual target evaluation value E are presented in interval form, which respectively stand for the worst and best case of the two values.

In terms of stability, Figures 4 to 6(a) illustrate that, during the evolution of the standard GA, the curves have some abnormal fluctuations. The same phenomenon does not appear on the other two methods. This shows that the proposed algorithm and the SAGA are more stable than the standard GA.

With regard to the convergence rate, we can see that the standard GA converges around the 22nd generation. The GA which only contains SA obtains the optimal solution in the 40th generation. And the idea value can be achieved near the 9th generation with our proposed algorithm. Besides, if the runtime is taken into account, the standard GA takes 0.0633 s per iteration. The SAGA takes 0.2630 s and our proposed algorithm takes 0.2536 s for each iteration step. By synthesizing the number of iterations above, we can see that it takes 1.39 s, 10.52 s, and 2.28 s, respectively, for the three algorithms to reach convergence. All these methods were realized using MATLAB code and ran on a LENOVO laptop (Win 10, Intel Core i5-3230M 2.6Hz CPU, 8 GB RAM).

In our experiment, the target killing expectation E is required to be as large as possible, while enemy residual threat F should be as small as possible. As shown in figures, for the standard GA, its optimal fitness interval is [1.08, 1.21]. The interval value of residual enemy threat F is [1.66, 1.87] and the target killing expectation E is 4.76. As for the proposed hybrid algorithm and SAGA, their values seem to be the same. Their optimal fitness and residual enemy threat interval values are [0.80, 0.89] and [1.23, 1.39], respectively. Both of them can achieve a target killing expected value of 5.14.

In order to compare the three algorithms more directly, their performance indicators mentioned above are listed in Table 6.

Scheme of autonomous target allocation

Table 6.

Performance comparison between the three algorithms.

	Standard genetic algorithm	Simulated annealing genetic algorithm	Improved hybrid genetic algorithm
Number of iterations of convergence	22	40	9
Single iteration time (s)	0.0633	0.2630	0.2536
Convergence time (s)	1.39	10.52	2.28
Optimal fitness	[1.08, 1.21]	[0.80, 0.89]	[0.80, 0.89]
Target residual threat (F)	[1.66, 1.87]	[1.23, 1.39]	[1.23, 1.39]
Target killing expectation (E)	4.76	5.14	5.14

$X=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$

In conclusion, compared with the standard GA, except for the better stability, the proposed hybrid algorithm obviously cost less number of iterations to reach convergence. Moreover, the proposed method can help achieve a higher target killing expected value (E) and make the residual target threat (F) smaller. This shows that our algorithm can effectively avoid local minima with a better search ability. Due to the properties of SA, although each iteration of our method takes longer, there is little difference in convergence time between the two. Besides, compared with the SAGA, both of them can achieve the ideal fitness value. However, the number of iterations of the SAGA reaching convergence is much more than the proposed algorithm. As their single iteration time is similar, the convergence time of the SAGA is much longer too. Therefore, the proposed method can better deal with MWMTA problem under uncertain environment.

According to the result matrix X, the first missile attacks enemy aircraft 2, the second and sixth missiles attack enemy aircraft 1, the third missile attacks enemy aircraft 4, the fourth and fifth missiles attack enemy aircraft 3, the seventh missile attacks enemy aircraft 5, and the eighth missile attacks enemy aircraft 6. Through the above simulation experiments, the interval grey number-based target allocation designed in this article can better adapt to the battlefield uncertainty without affecting the operational effectiveness, and obtain a reasonable attack plan at the same time.

Scenario 2: In order to verify the effectiveness of our proposed algorithm, another case where the number of missiles is less than the number of enemy targets is designed. Let our UAV formation contains three aircrafts and enemy formation has eight aircrafts. Each of our aircraft carries two missiles. Other parameters remain almost unchanged. Through the three algorithms, the results are shown in Figures 7 to 9.

Figure 7.

Optimal fitness expectation (Scenario 2). (a) Standard genetic algorithm. (b) Simulated annealing genetic algorithm. (c) Improved hybrid genetic algorithm.

Figure 8.

Residual target threat value F (Scenario 2). (a) Standard genetic algorithm. (b) Simulated annealing genetic algorithm. (c) Improved hybrid genetic algorithm.

Figure 9.

Target killing expected value E (Scenario 2). (a) Standard genetic algorithm. (b) Simulated annealing genetic algorithm. (c) Improved hybrid genetic algorithm.

Similarly, the experimental data of the three algorithms in scenario 2 are displayed in Table 7.

Table 7.

Experimental results in scenario 2.

	Standard genetic algorithm	Simulated annealing genetic algorithm	Improved hybrid genetic algorithm
Number of iterations of convergence	24	39	8
Single iteration time (s)	0.0515	0.2873	0.2717
Convergence time (s)	1.236	11.20	2.17
Optimal fitness	[2.14, 2.44]	[2.22, 2.55]	[1.89, 2.17]
Target residual threat (F)	[3.29, 3.78]	[3.41, 3.95]	[2.89, 3.36]
Target killing expectation (E)	4.66	4.46	4.93

From the table, we can see that only the proposed algorithm can obtain ideal values and it converges with the least number of iterations. Therefore, our method proves to have better performance than the other two methods. The target assignment scheme can be represented as

$X=\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$

The matrix indicates that missile 1 is allocated to attack target 6, missile 2 is to attack target 2, missile 3 is to attack target 4, missile 4 is assigned to target 5, missile 5 is assigned to target 7, and missile 6 is allocated to attack target 8.

Conclusion

This article uses UAVs as an example platform to study the threat assessment and MWMTA problems in uncertain environment. Based on the interval grey number, a modified GAHP is developed to determine the interval weight of each threat factor. According to grey theory, the possibility degree function is constructed to address the grey interval comparison problem, so that the interval weights of threat factors can be transformed into scalar-valued weights. On the premise of maximizing killing and minimizing casualties, a target assignment model is established. After that, a hybrid GA incorporating adaptive genetic operators and SA is proposed. Combined with the possibility degree function, appropriate fitness function and selection operation are designed to make the algorithm suitable for the optimization of the interval form fitness value. Moreover, improved adaptive crossover and mutation operators and the SA mechanism are presented to overcome the shortcomings of traditional GA. Finally, a reasonable attack scheme and the corresponding intelligent decision-making solution are derived. The simulation results prove satisfactory.