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

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. 1 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, 2 such as particle swarm algorithm 3 and improved ant colony algorithms. 4,5 In addition, Zaman et al. 6 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. 7 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, 8 MWMTA has become one of the key technologies for unmanned combat aircraft 9 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, 10 analytic hierarchy process (AHP), 11 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 12 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) 15 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, 16 ant colony optimization, 17 particle swarm optimization, 18 and genetic algorithm (GA) with greedy eugenics, 19 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) 20 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 offboresight 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 2. Velocity threat factor T v 3. Distance threat factor T r T r ¼ 0:5 r r mr ; r r mt 0:5 À 0:2 r À r mt r mr À r mt 0 @ 1 A r mt < r < r mr 1:0 r mr < r < r mt 0:8 max r mr ; r mt ð Þ< r < r dr 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 Y r k¼1 ð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 Â Y r k¼1 ð1 À p kj Þ x kj .
Accordingly, when we just consider the maximum expected kill value, the function can be derived as and when reducing casualties is taken into account separately, the remaining threat of the targets after a round of engagement can be expressed by 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, 27 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 28 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 .
Definition 2. The grey number with both lower bound a and upper bound a is called the interval grey number, denoted as 2 ½a; a. And when the lower and upper bounds are infinite, it is called the black number; when 2 ½a; a and a ¼ a, we call the white number. In addition, ðgÞ ¼ À þ ð À À þ Þg,ð0 g 1Þ is the standardized form of interval grey number.
Definition 3. The basic arithmetic of interval grey number is as follows.
Subtraction 1 À 2 2 ½a À d; b À c. Multiplication 1 Á 2 2 ½minfac; ad; bc; bdg; maxfac; ad; bc; bdg. Division is defined as 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.

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.

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][30][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¼ 11 After that, the weight vector of the interval grey number is calculated by here j ¼ n À 1; n À 2; Á Á Á; 1 and w n > 0. Thus the weight value of each index can be written as where w 1 is the weight of fighter combat capability, w 2 is the weight of fighter angle threat, w 3 is the weight of fighter distance threat,w 4 is the weight of fighter velocity threat, w 5 is the weight of fighter altitude threat, and w 6 is the weight of fighter stealth threat. Finally, through normalization, we can get the single interval weight

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 here the entry in the jth row and the ith column represents the evaluation value of expert i for indicator j.

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.

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 Afterward, sum the grey evaluation coefficients of target J's every evaluation grey clustering by 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 where K ¼ 1; 2; Á Á Á ; k, J ¼ 1; 2; Á Á Á ; j, relative to the evaluation index A, a grey evaluation matrix of the target can be derived as

Comprehensive evaluation
The synthetic grey evaluation vector 32 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

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, 2 ½ À ; þ is an interval grey number, the point possibility degree function can be defined by when a is less than .
Definition 5. Suppose 1 2 ½ 1 À ; 1 þ and 2 2 ½ 2 À ; 2 þ are two arbitrary interval grey numbers, ðgÞ is the standardization of interval grey number. If g 1 is given and 1 ðg 1 Þ is a fixed point, pð 1 ðg 1 Þ < 2 Þ would be the possibility degree function under the circumstance where 1 ðg 1 Þ is less than 2 . Therefore, the possibility degree function with 1 less than 2 can be expressed as which appears to be complementary, namely, According to the integral solution method given by the literature, 33 the possibility degree when 1 < 2 can be expressed as 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 1 ; 2 ; Á Á Á ; N are N grey intervals and the entry in the ith row and the jth column represents the possibility degree with grey interval i less than j .
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 ¼ ij Â Ã mÂ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 þ ::: þ 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 where F indicates that the minimum threat value of enemy residual targets, and E means destroying enemy targets as many as possible. Constraints s:t: X r k¼1 x kj s j ; j ¼ 1; 2;:::; n X r k¼1 X n j¼1 x kj r X n j¼1 x kj 1; k ¼ 1; 2;:::;r 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 34 is finally expressed by 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 MAX GEN, 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 0 that satisfy the constraint are generated. According to the above objective function model, the fitness value of the particle is calculated by 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 Â 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 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 20 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 35 and cosine function, an improved adaptive genetic operator is presented in this article. The functions of genetic operator are expressed as follows Here, P c max is the maximum crossover probability, P c min is the minimum crossover probability,P m max is the maximum mutation probability, P m min is the minimum mutation probability, A is a constant, f min is the minimum fitness value of the population, f avg is the average fitness value of the population, f 0 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 0 and coefficient g & ð0; 1Þ are defined. When T 0 > T min , that is, SA is not frozen, the optimal individual after mutation is denoted by P 0 and the new individual is denoted by P 1 . 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 ¼ expðÀDF=T Þ where DF ¼ FitðP 0 Þ þ À FitðP 1 Þ À . At last, the current evolution temperature can be calculated by T ¼ T Âg. Compare the calculated T with the baseline temperature T min , if T is less than T 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.
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.

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. Here T c j ; T a 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; 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 ¼ ½0:5; 1:5; 3:57; 5:5; 4:31; 2:62.
Then combined with the indexes given by experts, we can evaluate the sample matrix According to (11), the grey evaluation coefficient can be calculated by   Table 3. Dominant accumulation matrix under league system. Afterward, through (12) and (13), the grey evaluation weight of each indicator can be determined by At last, associated with the grey cluster rating D ¼ 9; 7; 5; 3; 1 ð Þ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 48dimensional 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 ! 1 and ! 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.
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.
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.  In order to compare the three algorithms more directly, their performance indicators mentioned above are listed in Table 6. Scheme of autonomous target allocation 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. Similarly, the experimental data of the three algorithms in scenario 2 are displayed in Table 7.
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    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.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (grant nos 61673209, 71971115, 61973158).