# Air Target Threat Assessment Based on Improved Moth Flame Optimization-Gray Neural Network Model

Authors: Longfei Yue; Rennong Yang; Jialiang Zuo; Hao Luo; Qiuliang Li
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Air target threat assessment is a key issue in air defense operations. Aiming at the shortcomings of traditional threat assessment methods, such as one-sided, subjective, and low-accuracy, a new method of air target threat assessment based on gray neural network model (GNNM) optimized by improved moth flame optimization (IMFO) algorithm is proposed. The model fully combines with excellent optimization performance of IMFO with powerful learning performance of GNNM. Finally, the model is trained and evaluated using the target threat database data. The simulation results show that compared with the GNNM model and the MFO-GNNM model, the proposed model has a mean square error of only 0.0012 when conducting threat assessment, which has higher accuracy and evaluates 25 groups of targets in 10 milliseconds, which meets real-time requirements. Therefore, the model can be effectively used for air target threat assessment.


---

## Body

## 1. Introduction

Air target threat assessment refers to comprehensively considering various factors affecting the target threat value, establishing a reasonable indicator system, and quantifying it, and then establishing a threat assessment model to evaluate the target threat value. As battlefield environment becomes more and more complex and incoming targets are characterized by high speed, high maneuverability, stealth, antijamming, and remote precision guidance, the pressure on air defense system is increasing. Therefore, establishing a reasonable threat assessment model and conducting a rapid and accurate threat assessment of incoming targets is a prerequisite for target allocation in air defense operations and an important support for efficient command and control.At present, a lot of research studies have been carried out on air target threat assessment at home and abroad. Two methods are mainly used. One is the method of reasoning and the other is the method of machine learning.The method of reasoning is to analyze the relationship between index values and threat values and combine expert experience and prior probability to infer the threat values of different targets. This method mainly includes the Bayesian network [1–3] and intuitionistic fuzzy logic [4–6]. The Bayesian network method can intuitively express the relationship between indicators and threats, but the determination of prior probability depends too much on expert experience. Therefore, subjectivity is too strong. An intuitionistic fuzzy logic method can better deal with qualitative indicators, but the rules of reasoning are complicated and the coupling relationship between the indicators is not considered enough.The machine learning method is to establish a threat assessment index system, quantify each index, and then learn the nonlinear relationship between the index values and the threat values through machine learning. This method mainly includes neural network [7, 8] and support vector machine [9]. Neural network has strong nonlinear fitting ability, but it needs large sample training. Support vector machine has advantages for small sample prediction, but it has disadvantages of unreasonable index selection and low evaluation accuracy.Therefore, how to choose a reasonable index system and scientific quantitative methods, while overcoming the strong subjective deficiencies of traditional methods and improving the accuracy of assessment have always been the focus of research on threat assessment.GNNM is a predictive model that combines the advantages of small sample prediction with gray system theory and the advantages of self-learning of neural network. Among them, gray system theory [10] is a systematic scientific theory proposed by Prof. Deng to predict the eigenvalues of uncertain system behaviors. Gray system modeling requires few samples, and the requirement for distribution law of data is not high. Neural network is a mathematical model mainly used to learn the nonlinear function relationship between input data and output data, with the advantages of parallel computing and self-learning, self-organization, and robustness. Combining gray theory with neural network to establish a gray neural network model can make full use of the advantages of both and improve the prediction accuracy. The performance of GNNM depends mainly on the choice of initial weights, so how to choose better initial weight parameters is particularly important.Moth flame optimization (MFO) algorithm [11] is a new intelligent optimization algorithm [12–14] proposed by Prof. Mirjalili in 2015. It has the advantages of less adjustment parameters, faster convergence, and simple implementation. Hence, it is often used in engineering optimization problems [15, 16]. However, MFO also has some shortcomings in the process of optimization. For example, the selection of current optimal solution in each iteration leads to trap in local optima and premature convergence, and the global optimization performance is poor due to low population diversity. Therefore, after improving the optimization performance of MFO, the optimal initial weight parameter is obtained by IMFO and the learning performance of GNNM is improved.Therefore, this paper combines the actual situation of air target threat assessment, selects a reasonable indicator system, establishes a threat assessment framework, and scientifically quantifies each indicator. Then, this paper improves the global optimization performance of the original MFO by introducing Tent chaos, Lévy flight, and Metropolis criterion. Then this paper combined IMFO with GNNM, a threat assessment model based on IMFO-GNNM is established. Finally, the proposed model is evaluated from accuracy and real-time performance. The results show that the proposed model has higher accuracy and better real-time performance in threat assessment. Therefore, this study has certain advantages and important practical significance.


## 2. Threat Assessment Problem Modeling

Air target threat assessment needs to analyze the air target and our asset to obtain the threat value of air target, thus providing a basis for our command and control. Due to the particularity of combat problem, on one hand, the accuracy of threat assessment should be considered, and on the other hand, the real-time performance of threat assessment should be considered. Therefore, it is crucial to choose a reasonable and accurate threat assessment indicator system and an efficient threat assessment method. Since air target threat assessment is not only related to the threat capability of target, but also related to the spatial location of target and the value of our asset. This paper considers threat level, threat capability, and threat extent to establish a target threat assessment framework, and it scientifically and reasonably quantifies each threat attribute. Assessment framework is shown in Figure1.Figure 1
Target threat assessment framework.


### 2.1. Threat Level

Since the value of asset on the battlefield is different, the threat value of air target is related to the value of our asset. The higher the value of our asset, the higher the target threat level and the higher the threat value. The value of our asset is generally related to the political, economic, and military value of the asset. Therefore, refer to the importance of asset evaluation standard [17], and a threat level evaluation function is established, as shown in the following equation:(1)Ci=ln∏j=1sCij+1ln10s+1=∑j=1slnCij+1ln10s+1.In the formula,Cij represents the value of the jth evaluation factor of our asset i and its value range Cij∈1,10, and s is the number of evaluation factors. The value of the asset is usually given by the superior command.


### 2.2. Threat Capability

The target threat capability includes the penetration ability of the target to our asset and the degree of damage caused by the successful penetration. Among them, the penetration ability is related to the target type and interference ability and the degree of damage is determined by the target damage ability.


#### 2.2.1. Target Type

There are many types of air targets. Different types of targets are different in size, radar cross-section (RCS), and threat values. Air targets are mainly classified into aircraft and guided weapons according to their size. Aircraft mainly includes bombers (BA), fighter bombers (FB), and armed helicopters (AH). Guided weapons mainly include precision guided bombs (PGB), air-to-ground missiles (AGM), cruise missiles (CM), tactical ballistic missiles (TBM), and antiradiation missiles (ARM). Therefore, the domain expert group quantifies the RCS size of different types of targets, and the results are shown in Table1.Table 1
Target type quantification.
Type
BA
FB
AH
PGB
AGM
CM
TBM
ARM
Threat
0.1
0.3
0.5
0.9
0.7
0.7
0.7
0.7


#### 2.2.2. Damage Ability

The target damage ability is determined by target damage probability and our survival probability. The damage probability is obtained through the intelligence data, and the survival probability is obtained by analyzing the factors such as the invulnerability and protection of our asset.(a)
When the air target is an aircraft type, the probability of damage is related to the type, quantity, and single-damage probability of the weapon carried by the aircraft. Therefore, its damage probability is the joint damage probability of different types of weapons, and its damage ability is quantified according to the following equation:(2)TdKiTj=1−1−pPGBiNPGBBA1−1−pAGMiNAGMBA1−1−pCMiNCMBApsurvivali,Tj∈BA,1−1−pPGBiNPGBFB1−1−pAGMiNAGMFB1−1−pARMiNARMFBpsurvivali,Tj∈FB,1−1−pPGBiNPGBAH1−1−pAGMiNAGMAHpsurvivali,Tj∈AH.In the formula,dKiTj represents the damage ability of the target j to the asset i, psurvivali represents the survival probability of the asset i, pPGBi represents the probability of damage of the PGB to the asset i, and NPGBBA represents the number of PGB mounted by the BA.(b)
When the air target is a guided weapon, its damage ability is quantified by(3)TdKiTj=pPGBipsurvivali,Tj∈PGB,pAGMipsurvivali,Tj∈AGM,pTBMipsurvivali,Tj∈TBM,pARMipsurvivali,Tj∈ARM,pPGBipsurvivali,Tj∈PGB.


#### 2.2.3. Interference Ability

According to the electronic interference pod and interference means, the target interference capability is classified into five levels: super strong, strong, medium, weak, and no. According to Miller’s quantitative theory [18], it is quantized to 0.9, 0.7, 0.5, 0.3, and 0.1.


### 2.3. Threat Extent

The target threat extent mainly measures the threat of the target to our asset from the spatial situation. As shown in Figure2, T represents the incoming target, K represents our asset, D represents the target distance, and θ represents the airway angle, that is, the angle between the projection of the target velocity in the horizontal plane and the connection between the target and the asset.Figure 2
Schematic diagram of relationship between target and asset.It can be seen from Figure2 that the closer the target distance, the higher the target speed and the higher the target threat value. In addition, considering killing boundary of the target, when the distance is fixed, the smaller the airway angle, the shorter the airway short cut, the more likely our asset is to fall into the target kill zone, and the higher the target threat value. Therefore, the target threat extent is measured by the three dimensions of the airway angle, the target speed, and the target distance.


#### 2.3.1. Airway Angle Threat

The airway angle usually varies in the range of [0°, 180°], but it is generally considered that when the airway angle is within [0°, 90°], the target is approaching to the asset. At this time, the larger the airway angle, the larger the airway short cut and the lesser the threat extent. Hence, the negative exponential function is used to quantify the airway angle. The formula is as follows:(4)Tθ=e−θKiTj,θKiTj∈0°,90°,0,else.


#### 2.3.2. Speed Threat

The higher the target speed, the shorter the response time and the higher the threat value of the target. Therefore, the speed threat is quantified by(5)Tv=1−e−0.005v.


#### 2.3.3. Distance Threat

The target distance is the distance between the target and the asset. The target distance directly affects the threat value of the target. The closer to our asset, the higher the threat and the farther the distance, the lower the threat value. Therefore, distance threat is quantified by(6)TD=maxD−DimaxD−minD.


## 2.1. Threat Level

Since the value of asset on the battlefield is different, the threat value of air target is related to the value of our asset. The higher the value of our asset, the higher the target threat level and the higher the threat value. The value of our asset is generally related to the political, economic, and military value of the asset. Therefore, refer to the importance of asset evaluation standard [17], and a threat level evaluation function is established, as shown in the following equation:(1)Ci=ln∏j=1sCij+1ln10s+1=∑j=1slnCij+1ln10s+1.In the formula,Cij represents the value of the jth evaluation factor of our asset i and its value range Cij∈1,10, and s is the number of evaluation factors. The value of the asset is usually given by the superior command.


## 2.2. Threat Capability

The target threat capability includes the penetration ability of the target to our asset and the degree of damage caused by the successful penetration. Among them, the penetration ability is related to the target type and interference ability and the degree of damage is determined by the target damage ability.


### 2.2.1. Target Type

There are many types of air targets. Different types of targets are different in size, radar cross-section (RCS), and threat values. Air targets are mainly classified into aircraft and guided weapons according to their size. Aircraft mainly includes bombers (BA), fighter bombers (FB), and armed helicopters (AH). Guided weapons mainly include precision guided bombs (PGB), air-to-ground missiles (AGM), cruise missiles (CM), tactical ballistic missiles (TBM), and antiradiation missiles (ARM). Therefore, the domain expert group quantifies the RCS size of different types of targets, and the results are shown in Table1.Table 1
Target type quantification.
Type
BA
FB
AH
PGB
AGM
CM
TBM
ARM
Threat
0.1
0.3
0.5
0.9
0.7
0.7
0.7
0.7


### 2.2.2. Damage Ability

The target damage ability is determined by target damage probability and our survival probability. The damage probability is obtained through the intelligence data, and the survival probability is obtained by analyzing the factors such as the invulnerability and protection of our asset.(a)
When the air target is an aircraft type, the probability of damage is related to the type, quantity, and single-damage probability of the weapon carried by the aircraft. Therefore, its damage probability is the joint damage probability of different types of weapons, and its damage ability is quantified according to the following equation:(2)TdKiTj=1−1−pPGBiNPGBBA1−1−pAGMiNAGMBA1−1−pCMiNCMBApsurvivali,Tj∈BA,1−1−pPGBiNPGBFB1−1−pAGMiNAGMFB1−1−pARMiNARMFBpsurvivali,Tj∈FB,1−1−pPGBiNPGBAH1−1−pAGMiNAGMAHpsurvivali,Tj∈AH.In the formula,dKiTj represents the damage ability of the target j to the asset i, psurvivali represents the survival probability of the asset i, pPGBi represents the probability of damage of the PGB to the asset i, and NPGBBA represents the number of PGB mounted by the BA.(b)
When the air target is a guided weapon, its damage ability is quantified by(3)TdKiTj=pPGBipsurvivali,Tj∈PGB,pAGMipsurvivali,Tj∈AGM,pTBMipsurvivali,Tj∈TBM,pARMipsurvivali,Tj∈ARM,pPGBipsurvivali,Tj∈PGB.


### 2.2.3. Interference Ability

According to the electronic interference pod and interference means, the target interference capability is classified into five levels: super strong, strong, medium, weak, and no. According to Miller’s quantitative theory [18], it is quantized to 0.9, 0.7, 0.5, 0.3, and 0.1.


## 2.2.1. Target Type

There are many types of air targets. Different types of targets are different in size, radar cross-section (RCS), and threat values. Air targets are mainly classified into aircraft and guided weapons according to their size. Aircraft mainly includes bombers (BA), fighter bombers (FB), and armed helicopters (AH). Guided weapons mainly include precision guided bombs (PGB), air-to-ground missiles (AGM), cruise missiles (CM), tactical ballistic missiles (TBM), and antiradiation missiles (ARM). Therefore, the domain expert group quantifies the RCS size of different types of targets, and the results are shown in Table1.Table 1
Target type quantification.
Type
BA
FB
AH
PGB
AGM
CM
TBM
ARM
Threat
0.1
0.3
0.5
0.9
0.7
0.7
0.7
0.7


## 2.2.2. Damage Ability

The target damage ability is determined by target damage probability and our survival probability. The damage probability is obtained through the intelligence data, and the survival probability is obtained by analyzing the factors such as the invulnerability and protection of our asset.(a)
When the air target is an aircraft type, the probability of damage is related to the type, quantity, and single-damage probability of the weapon carried by the aircraft. Therefore, its damage probability is the joint damage probability of different types of weapons, and its damage ability is quantified according to the following equation:(2)TdKiTj=1−1−pPGBiNPGBBA1−1−pAGMiNAGMBA1−1−pCMiNCMBApsurvivali,Tj∈BA,1−1−pPGBiNPGBFB1−1−pAGMiNAGMFB1−1−pARMiNARMFBpsurvivali,Tj∈FB,1−1−pPGBiNPGBAH1−1−pAGMiNAGMAHpsurvivali,Tj∈AH.In the formula,dKiTj represents the damage ability of the target j to the asset i, psurvivali represents the survival probability of the asset i, pPGBi represents the probability of damage of the PGB to the asset i, and NPGBBA represents the number of PGB mounted by the BA.(b)
When the air target is a guided weapon, its damage ability is quantified by(3)TdKiTj=pPGBipsurvivali,Tj∈PGB,pAGMipsurvivali,Tj∈AGM,pTBMipsurvivali,Tj∈TBM,pARMipsurvivali,Tj∈ARM,pPGBipsurvivali,Tj∈PGB.


## 2.2.3. Interference Ability

According to the electronic interference pod and interference means, the target interference capability is classified into five levels: super strong, strong, medium, weak, and no. According to Miller’s quantitative theory [18], it is quantized to 0.9, 0.7, 0.5, 0.3, and 0.1.


## 2.3. Threat Extent

The target threat extent mainly measures the threat of the target to our asset from the spatial situation. As shown in Figure2, T represents the incoming target, K represents our asset, D represents the target distance, and θ represents the airway angle, that is, the angle between the projection of the target velocity in the horizontal plane and the connection between the target and the asset.Figure 2
Schematic diagram of relationship between target and asset.It can be seen from Figure2 that the closer the target distance, the higher the target speed and the higher the target threat value. In addition, considering killing boundary of the target, when the distance is fixed, the smaller the airway angle, the shorter the airway short cut, the more likely our asset is to fall into the target kill zone, and the higher the target threat value. Therefore, the target threat extent is measured by the three dimensions of the airway angle, the target speed, and the target distance.


### 2.3.1. Airway Angle Threat

The airway angle usually varies in the range of [0°, 180°], but it is generally considered that when the airway angle is within [0°, 90°], the target is approaching to the asset. At this time, the larger the airway angle, the larger the airway short cut and the lesser the threat extent. Hence, the negative exponential function is used to quantify the airway angle. The formula is as follows:(4)Tθ=e−θKiTj,θKiTj∈0°,90°,0,else.


### 2.3.2. Speed Threat

The higher the target speed, the shorter the response time and the higher the threat value of the target. Therefore, the speed threat is quantified by(5)Tv=1−e−0.005v.


### 2.3.3. Distance Threat

The target distance is the distance between the target and the asset. The target distance directly affects the threat value of the target. The closer to our asset, the higher the threat and the farther the distance, the lower the threat value. Therefore, distance threat is quantified by(6)TD=maxD−DimaxD−minD.


## 2.3.1. Airway Angle Threat

The airway angle usually varies in the range of [0°, 180°], but it is generally considered that when the airway angle is within [0°, 90°], the target is approaching to the asset. At this time, the larger the airway angle, the larger the airway short cut and the lesser the threat extent. Hence, the negative exponential function is used to quantify the airway angle. The formula is as follows:(4)Tθ=e−θKiTj,θKiTj∈0°,90°,0,else.


## 2.3.2. Speed Threat

The higher the target speed, the shorter the response time and the higher the threat value of the target. Therefore, the speed threat is quantified by(5)Tv=1−e−0.005v.


## 2.3.3. Distance Threat

The target distance is the distance between the target and the asset. The target distance directly affects the threat value of the target. The closer to our asset, the higher the threat and the farther the distance, the lower the threat value. Therefore, distance threat is quantified by(6)TD=maxD−DimaxD−minD.


## 3. Construction of IMFO-GNNM

### 3.1. GNNM

Suppose there is a sequence(7)x0=x01,x02,…,x0n.We can get a new sequence by accumulating (AGO):(8)x1=x11,x12,…,x1n=x01,x01+x02,…,x01+x02+⋯+x0n.It is easy to knowx1 is exponential growth law, so according to the idea of derivative and difference, we can construct differential equations to fit x1.For the sake of convenience, the original sequencex0 is recorded as xt, the new sequence x1 is recorded as yt, and the prediction result is zt; so, the differential equation expression of the gray neural network with n parameters is as follows:(9)dy1dt+ay1=b1y2+b2y3+⋯+bn−1yn.In the formula, the output variable isy1, the input variable is yii=2,…,n, and a and bi are parameters.Solving the differential equation we can get(10)zt=y10−b1ay2t−b2ay3t−⋯−bn−1aynte−at+b1ay2t+b2ay3t+⋯+bn−1aynt.Ifd=b1/ay2t+b2/ay3t+⋯+bn−1/aynt, we will derive(11)zt=y10−d⋅e−at1+e−at+d⋅11+e−at⋅1+e−at=y10−d⋅1−11+e−at+d⋅11+e−at⋅1+e−at=y10−d−y10⋅11+e−at+2d⋅11+e−at⋅1+e−at.Then we map it to the neural network to get a gray neural network. Its structure is shown in Figure3.Figure 3
Network structure of GNNM.Network output is(12)y1=w31c1+w32c2+⋯+w3ncn−θy1.The weight of the input layer isw11=a, the weight of the input layer and the hidden layer is w21=−x10, w2i=2bi−1/a,i=2,…,n, the weight of the hidden layer and the output layer is w3j=a+e−at,j=1,…,n, the output of the input layer is h=1/1+e−at, the output of the hidden layer is c1=hw21, ci=yithw2i,i=2,…,n, and the deviation of the output layer is θy1=1+e−atd−y10. Note that when the input variable is n − 1, the number of hidden layer nodes is n.


### 3.2. IMFO

#### 3.2.1. MFO

The MFO was inspired by the lateral positioning of moths during night flight [11]. The moths often keep the same angle with the moon when flying. Because they are far away from the moon, the moth and the moon are connected in parallel at different times. So, the moths can fly in a straight line. However, in reality, the moth is very close to the flame, and the moth still maintains the angle with the flame, so that the moth flies along the equiangular spiral to the flame. Hence, there is a “moth to flame,” as shown in Figure 4. Prof. Mirjalili was inspired by this isometric spiral function to invent the MFO algorithm.Figure 4
Lateral positioning flight and moth to flame.In the original MFO algorithm,M is the moth matrix, that is, the agent that searches for the solution, ZM is the moth fitness value matrix, F is the flame position matrix, that is, the current optimal solution, and ZF is the flame fitness value matrix, as shown in the formula (13) and (14). At the initial moment, the flame matrix and the moth matrix have the same dimensions:(13)M=m1,1m1,2⋯m1,dm2,1m2,2⋯m2,d⋮⋮⋮⋮mn,1mn,2⋯mn,d,ZM=zm1zm2⋮zmn,(14)F=f1,1f1,2⋯f1,df2,1f2,2⋯f2,d⋮⋮⋮⋮fn,1fn,2⋯fn,d,ZF=zf1zf2⋮zfn,where n is the number of moths and d is the number of dimensions.The flame positions are sorted according to the fitness values from small to large, and the moths fly around the sorted flames, along the isometric spiral of equation (15), and update their position by changing the parameter t:(15)Mi=Diebtcos2πt+Fj,where Mi is the new position of the ith moth, Di, formulated as Di=Mi−Fj, is the distance between the ith moth and the jth flame, b is the isometric spiral parameter, which determines the shape of the isometric spiral, and t is a random number between [−1, 1], controlling the distance between the moth and the flame, as shown in Figure 5. The smaller the t is, the closer the moth is to the flame. By changing t, the moth can reach the position around the flame, thus enhancing the local optimization performance of the algorithm.Figure 5
The position of moth changes with time.In order to improve the search efficiency of the algorithm, MFO adopts the flame adaptive reduction mechanism, which makes the later convergence faster, without continuing to optimize around the inferior solution. The number of flames is adaptively reduced by(16)NF=roundNmax−I×Nmax−1Imax,where NF is the current number of flames, Nmax is the maximum number of flames, I is the current number of iterations, and Imax is the maximum number of iterations. The number of flames varies with the number of iterations as shown in Figure 6.Figure 6
The number of flames decreases with the number of iterations.When the number of flames is less than the number of moths, the moths fly around the worst-fitted flame.


#### 3.2.2. Tent Chaos

The chaotic sequence [19–22] has ergodicity and randomness. The chaotic sequence can make the initial solution distribution more uniform, which is beneficial to find the global optimization solution. At present, two chaotic sequences are mainly used, namely logistic mapping and tent mapping. Shan et al. [23] proved that tent mapping is better ergodic than logistic mapping. Therefore, the tent mapping is used to optimize the initial population in this paper. Tent mapping and logistic mapping distribution are shown in Figure 7. Figure 7(a) is a tent map result of mapping 400 times from 0.01 to 0.5 according to the tent mapping formula, and Figure 7(b) is a logistic map result of mapping 200 times from 3 to 4 according to the logistic mapping formula.Figure 7
Tent map and logistic map. (a) Tent map. (b) Logistic map.
(a)
(b)The Tent mapping formula is as follows:(17)xt+1=2xt,0≤xt≤0.5,21−xt,0.5<xt≤1.


#### 3.2.3. Lévy Flights

Since the original MFO algorithm selects the current optimality for each iteration and is prone to trap in local optima, Lévy flight is introduced. Lévy flight, different from the Brownian motion, is a random motion proposed by French mathematician Paul Pierre Lévy [24–26], which is characterized by a large number of small step movements mixed with a small amount of large step movement, as shown in Figure 8. Figure 8 is the result of flying 3000 times according to the Lévy formula. It can also be seen from this picture that Lévy flight has large steps and small steps.Figure 8
Lévy flights path.The step size of each step movement and direction obeys Lévy distribution [27], as shown in the following equation:(18)L∼μ×uv1/β,μ=Γ1+β×sinπ×β/2×randΓ1+β/2×β×2β−1/21/β.Among them,Γλ is a standard gamma function, u and v obey the normal distribution, and β is a constant.The current optimal solution updates as follows:(19)xit+1=xit+αL,where α is a scale factor that controls the learning step size for each step. In this way, after moving a small step, it will also move a large step size, which balances exploration and exploitation of the search space and helps to jump out of the local optimal solution.


#### 3.2.4. Metropolis Criterion

The Metropolis criterion is the core of the simulated annealing algorithm. It accepts the current inferior solution with a certain probability, so that the algorithm has the ability to jump out of the local optimum and converge to the global optimal solution.The Metropolis criterion uses the following process (take the minimization problem as an example):(1)
Adding the Lévy flight to the current optimal solutionS1 to generate a new solution S2.(2)
Weight the new solution and the old solution according to equation (20) to get the final new solution. The new solution obtained in this way is beneficial to retain the advantages of the old solution and reduce the influence of the disturbance error:(20)S2=alpha×S2+1−alpha×S1.(3)
Calculate the increment,df=fS2−fS1, fS1 is the cost function of S1.(4)
According to the Metropolis criterion, the formula is expressed as follows:(21)P=1,df<0,e−df/T,df≥0.Ifdf<0, we accept the new solution; otherwise, we accept the new solution with the probability of e−df/T.


#### 3.2.5. IMFO

The IMFO algorithm is given in Algorithm1.Algorithm 1:Improvement moth flame optimization algorithm.
(1)
Initialize solution population using tent chaos map
(2)
iteration = 1
(3)
while (iteration ≤ Max_iteration)
(4)
OM = FitnessFunction(M)
(5)
if iteration = = 1
(6)
F = sort(M)
(7)
OF = sort(OM)
(8)
else
(9)
F = sort(Mt−1, Mt)
(10)
OF = sort(Mt-1, Mt)
(11)
end if
(12)
fori = 1 : n
(13)
forj = 1 : d
(14)
updatet
(15)
calculateD with respect to the corresponding flame
(16)
update M(i, j) using equation (15) with respect to the corresponding flame
(17)
end for
(18)
end for
(19)
update the position of the current optimal agent using Lévy-flight
(20)
F_lévy = Lévy(F)
(21)
OF_lévy = FitnessFunction(F_lévy)
(22)
using the Metropolis criterion for OF and OF_lévy
(23)
update the position best flame obtained so far
(24)
update flame number using equation (16)
(25)
iteration = iteration + 1
(26)
end while


### 3.3. IMFO-GNNM

Because the neural network learns the functional relationship between input and output by adjusting the weight connection, the learning effect depends on its weight parameters. In order to enhance the learning effect of GNNM, IMFO is used to optimize its weight parameters, and then GNNM with optimized weights is used to predict. Therefore, the IMFO-GNNM algorithm is obtained, and the algorithm flow is shown in Figure9.Figure 9
Flow chart for the proposed algorithm.


## 3.1. GNNM

Suppose there is a sequence(7)x0=x01,x02,…,x0n.We can get a new sequence by accumulating (AGO):(8)x1=x11,x12,…,x1n=x01,x01+x02,…,x01+x02+⋯+x0n.It is easy to knowx1 is exponential growth law, so according to the idea of derivative and difference, we can construct differential equations to fit x1.For the sake of convenience, the original sequencex0 is recorded as xt, the new sequence x1 is recorded as yt, and the prediction result is zt; so, the differential equation expression of the gray neural network with n parameters is as follows:(9)dy1dt+ay1=b1y2+b2y3+⋯+bn−1yn.In the formula, the output variable isy1, the input variable is yii=2,…,n, and a and bi are parameters.Solving the differential equation we can get(10)zt=y10−b1ay2t−b2ay3t−⋯−bn−1aynte−at+b1ay2t+b2ay3t+⋯+bn−1aynt.Ifd=b1/ay2t+b2/ay3t+⋯+bn−1/aynt, we will derive(11)zt=y10−d⋅e−at1+e−at+d⋅11+e−at⋅1+e−at=y10−d⋅1−11+e−at+d⋅11+e−at⋅1+e−at=y10−d−y10⋅11+e−at+2d⋅11+e−at⋅1+e−at.Then we map it to the neural network to get a gray neural network. Its structure is shown in Figure3.Figure 3
Network structure of GNNM.Network output is(12)y1=w31c1+w32c2+⋯+w3ncn−θy1.The weight of the input layer isw11=a, the weight of the input layer and the hidden layer is w21=−x10, w2i=2bi−1/a,i=2,…,n, the weight of the hidden layer and the output layer is w3j=a+e−at,j=1,…,n, the output of the input layer is h=1/1+e−at, the output of the hidden layer is c1=hw21, ci=yithw2i,i=2,…,n, and the deviation of the output layer is θy1=1+e−atd−y10. Note that when the input variable is n − 1, the number of hidden layer nodes is n.


## 3.2. IMFO

### 3.2.1. MFO

The MFO was inspired by the lateral positioning of moths during night flight [11]. The moths often keep the same angle with the moon when flying. Because they are far away from the moon, the moth and the moon are connected in parallel at different times. So, the moths can fly in a straight line. However, in reality, the moth is very close to the flame, and the moth still maintains the angle with the flame, so that the moth flies along the equiangular spiral to the flame. Hence, there is a “moth to flame,” as shown in Figure 4. Prof. Mirjalili was inspired by this isometric spiral function to invent the MFO algorithm.Figure 4
Lateral positioning flight and moth to flame.In the original MFO algorithm,M is the moth matrix, that is, the agent that searches for the solution, ZM is the moth fitness value matrix, F is the flame position matrix, that is, the current optimal solution, and ZF is the flame fitness value matrix, as shown in the formula (13) and (14). At the initial moment, the flame matrix and the moth matrix have the same dimensions:(13)M=m1,1m1,2⋯m1,dm2,1m2,2⋯m2,d⋮⋮⋮⋮mn,1mn,2⋯mn,d,ZM=zm1zm2⋮zmn,(14)F=f1,1f1,2⋯f1,df2,1f2,2⋯f2,d⋮⋮⋮⋮fn,1fn,2⋯fn,d,ZF=zf1zf2⋮zfn,where n is the number of moths and d is the number of dimensions.The flame positions are sorted according to the fitness values from small to large, and the moths fly around the sorted flames, along the isometric spiral of equation (15), and update their position by changing the parameter t:(15)Mi=Diebtcos2πt+Fj,where Mi is the new position of the ith moth, Di, formulated as Di=Mi−Fj, is the distance between the ith moth and the jth flame, b is the isometric spiral parameter, which determines the shape of the isometric spiral, and t is a random number between [−1, 1], controlling the distance between the moth and the flame, as shown in Figure 5. The smaller the t is, the closer the moth is to the flame. By changing t, the moth can reach the position around the flame, thus enhancing the local optimization performance of the algorithm.Figure 5
The position of moth changes with time.In order to improve the search efficiency of the algorithm, MFO adopts the flame adaptive reduction mechanism, which makes the later convergence faster, without continuing to optimize around the inferior solution. The number of flames is adaptively reduced by(16)NF=roundNmax−I×Nmax−1Imax,where NF is the current number of flames, Nmax is the maximum number of flames, I is the current number of iterations, and Imax is the maximum number of iterations. The number of flames varies with the number of iterations as shown in Figure 6.Figure 6
The number of flames decreases with the number of iterations.When the number of flames is less than the number of moths, the moths fly around the worst-fitted flame.


### 3.2.2. Tent Chaos

The chaotic sequence [19–22] has ergodicity and randomness. The chaotic sequence can make the initial solution distribution more uniform, which is beneficial to find the global optimization solution. At present, two chaotic sequences are mainly used, namely logistic mapping and tent mapping. Shan et al. [23] proved that tent mapping is better ergodic than logistic mapping. Therefore, the tent mapping is used to optimize the initial population in this paper. Tent mapping and logistic mapping distribution are shown in Figure 7. Figure 7(a) is a tent map result of mapping 400 times from 0.01 to 0.5 according to the tent mapping formula, and Figure 7(b) is a logistic map result of mapping 200 times from 3 to 4 according to the logistic mapping formula.Figure 7
Tent map and logistic map. (a) Tent map. (b) Logistic map.
(a)
(b)The Tent mapping formula is as follows:(17)xt+1=2xt,0≤xt≤0.5,21−xt,0.5<xt≤1.


### 3.2.3. Lévy Flights

Since the original MFO algorithm selects the current optimality for each iteration and is prone to trap in local optima, Lévy flight is introduced. Lévy flight, different from the Brownian motion, is a random motion proposed by French mathematician Paul Pierre Lévy [24–26], which is characterized by a large number of small step movements mixed with a small amount of large step movement, as shown in Figure 8. Figure 8 is the result of flying 3000 times according to the Lévy formula. It can also be seen from this picture that Lévy flight has large steps and small steps.Figure 8
Lévy flights path.The step size of each step movement and direction obeys Lévy distribution [27], as shown in the following equation:(18)L∼μ×uv1/β,μ=Γ1+β×sinπ×β/2×randΓ1+β/2×β×2β−1/21/β.Among them,Γλ is a standard gamma function, u and v obey the normal distribution, and β is a constant.The current optimal solution updates as follows:(19)xit+1=xit+αL,where α is a scale factor that controls the learning step size for each step. In this way, after moving a small step, it will also move a large step size, which balances exploration and exploitation of the search space and helps to jump out of the local optimal solution.


### 3.2.4. Metropolis Criterion

The Metropolis criterion is the core of the simulated annealing algorithm. It accepts the current inferior solution with a certain probability, so that the algorithm has the ability to jump out of the local optimum and converge to the global optimal solution.The Metropolis criterion uses the following process (take the minimization problem as an example):(1)
Adding the Lévy flight to the current optimal solutionS1 to generate a new solution S2.(2)
Weight the new solution and the old solution according to equation (20) to get the final new solution. The new solution obtained in this way is beneficial to retain the advantages of the old solution and reduce the influence of the disturbance error:(20)S2=alpha×S2+1−alpha×S1.(3)
Calculate the increment,df=fS2−fS1, fS1 is the cost function of S1.(4)
According to the Metropolis criterion, the formula is expressed as follows:(21)P=1,df<0,e−df/T,df≥0.Ifdf<0, we accept the new solution; otherwise, we accept the new solution with the probability of e−df/T.


### 3.2.5. IMFO

The IMFO algorithm is given in Algorithm1.Algorithm 1:Improvement moth flame optimization algorithm.
(1)
Initialize solution population using tent chaos map
(2)
iteration = 1
(3)
while (iteration ≤ Max_iteration)
(4)
OM = FitnessFunction(M)
(5)
if iteration = = 1
(6)
F = sort(M)
(7)
OF = sort(OM)
(8)
else
(9)
F = sort(Mt−1, Mt)
(10)
OF = sort(Mt-1, Mt)
(11)
end if
(12)
fori = 1 : n
(13)
forj = 1 : d
(14)
updatet
(15)
calculateD with respect to the corresponding flame
(16)
update M(i, j) using equation (15) with respect to the corresponding flame
(17)
end for
(18)
end for
(19)
update the position of the current optimal agent using Lévy-flight
(20)
F_lévy = Lévy(F)
(21)
OF_lévy = FitnessFunction(F_lévy)
(22)
using the Metropolis criterion for OF and OF_lévy
(23)
update the position best flame obtained so far
(24)
update flame number using equation (16)
(25)
iteration = iteration + 1
(26)
end while


## 3.2.1. MFO

The MFO was inspired by the lateral positioning of moths during night flight [11]. The moths often keep the same angle with the moon when flying. Because they are far away from the moon, the moth and the moon are connected in parallel at different times. So, the moths can fly in a straight line. However, in reality, the moth is very close to the flame, and the moth still maintains the angle with the flame, so that the moth flies along the equiangular spiral to the flame. Hence, there is a “moth to flame,” as shown in Figure 4. Prof. Mirjalili was inspired by this isometric spiral function to invent the MFO algorithm.Figure 4
Lateral positioning flight and moth to flame.In the original MFO algorithm,M is the moth matrix, that is, the agent that searches for the solution, ZM is the moth fitness value matrix, F is the flame position matrix, that is, the current optimal solution, and ZF is the flame fitness value matrix, as shown in the formula (13) and (14). At the initial moment, the flame matrix and the moth matrix have the same dimensions:(13)M=m1,1m1,2⋯m1,dm2,1m2,2⋯m2,d⋮⋮⋮⋮mn,1mn,2⋯mn,d,ZM=zm1zm2⋮zmn,(14)F=f1,1f1,2⋯f1,df2,1f2,2⋯f2,d⋮⋮⋮⋮fn,1fn,2⋯fn,d,ZF=zf1zf2⋮zfn,where n is the number of moths and d is the number of dimensions.The flame positions are sorted according to the fitness values from small to large, and the moths fly around the sorted flames, along the isometric spiral of equation (15), and update their position by changing the parameter t:(15)Mi=Diebtcos2πt+Fj,where Mi is the new position of the ith moth, Di, formulated as Di=Mi−Fj, is the distance between the ith moth and the jth flame, b is the isometric spiral parameter, which determines the shape of the isometric spiral, and t is a random number between [−1, 1], controlling the distance between the moth and the flame, as shown in Figure 5. The smaller the t is, the closer the moth is to the flame. By changing t, the moth can reach the position around the flame, thus enhancing the local optimization performance of the algorithm.Figure 5
The position of moth changes with time.In order to improve the search efficiency of the algorithm, MFO adopts the flame adaptive reduction mechanism, which makes the later convergence faster, without continuing to optimize around the inferior solution. The number of flames is adaptively reduced by(16)NF=roundNmax−I×Nmax−1Imax,where NF is the current number of flames, Nmax is the maximum number of flames, I is the current number of iterations, and Imax is the maximum number of iterations. The number of flames varies with the number of iterations as shown in Figure 6.Figure 6
The number of flames decreases with the number of iterations.When the number of flames is less than the number of moths, the moths fly around the worst-fitted flame.


## 3.2.2. Tent Chaos

The chaotic sequence [19–22] has ergodicity and randomness. The chaotic sequence can make the initial solution distribution more uniform, which is beneficial to find the global optimization solution. At present, two chaotic sequences are mainly used, namely logistic mapping and tent mapping. Shan et al. [23] proved that tent mapping is better ergodic than logistic mapping. Therefore, the tent mapping is used to optimize the initial population in this paper. Tent mapping and logistic mapping distribution are shown in Figure 7. Figure 7(a) is a tent map result of mapping 400 times from 0.01 to 0.5 according to the tent mapping formula, and Figure 7(b) is a logistic map result of mapping 200 times from 3 to 4 according to the logistic mapping formula.Figure 7
Tent map and logistic map. (a) Tent map. (b) Logistic map.
(a)
(b)The Tent mapping formula is as follows:(17)xt+1=2xt,0≤xt≤0.5,21−xt,0.5<xt≤1.


## 3.2.3. Lévy Flights

Since the original MFO algorithm selects the current optimality for each iteration and is prone to trap in local optima, Lévy flight is introduced. Lévy flight, different from the Brownian motion, is a random motion proposed by French mathematician Paul Pierre Lévy [24–26], which is characterized by a large number of small step movements mixed with a small amount of large step movement, as shown in Figure 8. Figure 8 is the result of flying 3000 times according to the Lévy formula. It can also be seen from this picture that Lévy flight has large steps and small steps.Figure 8
Lévy flights path.The step size of each step movement and direction obeys Lévy distribution [27], as shown in the following equation:(18)L∼μ×uv1/β,μ=Γ1+β×sinπ×β/2×randΓ1+β/2×β×2β−1/21/β.Among them,Γλ is a standard gamma function, u and v obey the normal distribution, and β is a constant.The current optimal solution updates as follows:(19)xit+1=xit+αL,where α is a scale factor that controls the learning step size for each step. In this way, after moving a small step, it will also move a large step size, which balances exploration and exploitation of the search space and helps to jump out of the local optimal solution.


## 3.2.4. Metropolis Criterion

The Metropolis criterion is the core of the simulated annealing algorithm. It accepts the current inferior solution with a certain probability, so that the algorithm has the ability to jump out of the local optimum and converge to the global optimal solution.The Metropolis criterion uses the following process (take the minimization problem as an example):(1)
Adding the Lévy flight to the current optimal solutionS1 to generate a new solution S2.(2)
Weight the new solution and the old solution according to equation (20) to get the final new solution. The new solution obtained in this way is beneficial to retain the advantages of the old solution and reduce the influence of the disturbance error:(20)S2=alpha×S2+1−alpha×S1.(3)
Calculate the increment,df=fS2−fS1, fS1 is the cost function of S1.(4)
According to the Metropolis criterion, the formula is expressed as follows:(21)P=1,df<0,e−df/T,df≥0.Ifdf<0, we accept the new solution; otherwise, we accept the new solution with the probability of e−df/T.


## 3.2.5. IMFO

The IMFO algorithm is given in Algorithm1.Algorithm 1:Improvement moth flame optimization algorithm.
(1)
Initialize solution population using tent chaos map
(2)
iteration = 1
(3)
while (iteration ≤ Max_iteration)
(4)
OM = FitnessFunction(M)
(5)
if iteration = = 1
(6)
F = sort(M)
(7)
OF = sort(OM)
(8)
else
(9)
F = sort(Mt−1, Mt)
(10)
OF = sort(Mt-1, Mt)
(11)
end if
(12)
fori = 1 : n
(13)
forj = 1 : d
(14)
updatet
(15)
calculateD with respect to the corresponding flame
(16)
update M(i, j) using equation (15) with respect to the corresponding flame
(17)
end for
(18)
end for
(19)
update the position of the current optimal agent using Lévy-flight
(20)
F_lévy = Lévy(F)
(21)
OF_lévy = FitnessFunction(F_lévy)
(22)
using the Metropolis criterion for OF and OF_lévy
(23)
update the position best flame obtained so far
(24)
update flame number using equation (16)
(25)
iteration = iteration + 1
(26)
end while


## 3.3. IMFO-GNNM

Because the neural network learns the functional relationship between input and output by adjusting the weight connection, the learning effect depends on its weight parameters. In order to enhance the learning effect of GNNM, IMFO is used to optimize its weight parameters, and then GNNM with optimized weights is used to predict. Therefore, the IMFO-GNNM algorithm is obtained, and the algorithm flow is shown in Figure9.Figure 9
Flow chart for the proposed algorithm.


## 4. Results and Discussion

### 4.1. Verification Tests by Benchmark Functions

This paper selects six classical benchmark functions in CEC2010 to test the performance of the algorithm, as shown in Table2. Among them, the Sphere function is a smooth monotonic function with only one global minimum, which is used to test the convergence speed of the algorithm. The Rosenbrock function is a nonconvex unimodal function with multiple local extremums, which is mainly used to test the convergence speed and optimization accuracy of the algorithm; Schwefel 2.26, Rastrigin, Ackley, and Griewank functions are complex multimodal functions with multiple local minima. They are mainly used to test the global optimization performance of the algorithm and ability to jump out of local extremums. These functions can more comprehensively examine the optimization performance of the algorithm.Table 2
Benchmark functions.
Function
Equation and range
Range
Optimum
Property
Sphere
f
X
=
∑
i
=
1
D
x
i
2
[−100, 100]
0
Unimodal
Rosenbrock
f
X
=
∑
i
=
1
D
−
1
100
x
i
+
1
−
x
i
2
2
+
x
i
−
1
2
[−30, 30]
0
Unimodal
Schwefel 2.26
f
X
=
−
∑
i
=
1
D
x
i
sin
x
i
[−500, 500]
−418.98
Multimodal
Rastrigin
f
X
=
∑
i
=
1
D
100
x
i
2
−
10
cos
2
π
x
i
+
10
[−5.12, 5.12]
0
Multimodal
Ackley
f
X
=
−
20
exp
−
0.2
1
/
n
∑
i
=
1
D
x
i
2
−
exp
1
/
n
∑
i
=
1
D
cos
2
π
x
i
+
20
+
e
[−30, 30]
0
Multimodal
Griewank
f
X
=
1
/
1000
∑
i
=
1
D
x
i
2
−
∏
i
=
1
D
cos
x
i
/
i
+
1
[−600, 600]
0
MultimodalIn order to evaluate the effectiveness of the IMFO algorithm, this paper uses the above benchmark function to test and compare it with the MFO algorithm. The simulation was carried out with MATLAB R2014a. The parameters were set as follows: the number of moths is 30, the maximum number of iterations is 1000, the helix parameter is 1, the Lévy flight parameter isβ = 1.5, and the weighting coefficient of the new solution and the old solution is 0.5. In order to reduce the error, the fitness mean obtained by continuously running 30 times is taken as the test result.The fitness convergence curves of the 10-dimensional benchmark function test are shown in Figure10–15. The fitness mean and standard deviation of the 10-dimensional and 50-dimensional benchmark function tests are shown in Table 3.Figure 10
Sphere (10 dimensions).Figure 11
Rosenbrock (10 dimensions).Figure 12
Schwefel (10 dimensions).Figure 13
Rastrigin (10 dimensions).Figure 14
Ackley (10 dimensions).Figure 15
Griewank (10 dimensions).Table 3
Optimization results comparison of benchmark functions.
Function
Dimension
Algorithm
Mean
Standard deviation
Sphere
10
MFO
8.738 × 10−35
8.043 × 10−35
IMFO
6.544 × 10−55
1.936 × 10−50
50
MFO
9.965 × 10−3
4.507 × 10−2
IMFO
3.346 × 10−50
1.672 × 10−48
Rosenbrock
10
MFO
7.971
1.463
IMFO
8.954 × 10–1
3.041 × 10−1
50
MFO
4.873 × 101
1.676
IMFO
9.580 × 10−1
1.637
Schwefel 2.26
10
MFO
−0.999 × 103
1.649 × 103
IMFO
−4.188 × 103
1.403
50
MFO
−1.392 × 104
7.090 × 102
IMFO
−2.094 × 104
5.981 × 10−2
Rastrigin
10
MFO
5.597
6.085
IMFO
4.967 × 10−16
3.455 × 10−15
50
MFO
7.355 × 101
9.458
IMFO
3.678 × 10−14
5.279 × 10−11
Ackley
10
MFO
4.209 × 10−15
2.621 × 10−16
IMFO
8.419 × 10−16
3.847 × 10−17
50
MFO
7.118 × 10−2
2.937 × 10−1
IMFO
8.882 × 10−16
5.580 × 10−18
Griewank
10
MFO
1.314 × 10−4
8.359 × 10–4
IMFO
3.177 × 10−17
3.564 × 10−16
50
MFO
1.129 × 10−4
9.423 × 10–4
IMFO
1.761 × 10−17
2.965 × 10−15It can be seen from Figure10–15 that for the unimodal function, the IMFO algorithm converges faster than the MFO algorithm, and the optimization precision is higher. For the multimodal function, the MFO algorithm basically traps in local optima and the convergence speed is slow. However, the IMFO algorithm can jump out of the local optimal trap and search for the optimal solution all the time. It has good global optimization performance and faster convergence speed.It can be seen from Table3 that the mean and standard deviation of the optimization of the IMFO algorithm are smaller than the mean and standard deviation of the MFO algorithm, which proves that the IMFO algorithm has better performance and good robustness. Especially for the Schwefel function and the Rastrigin function, the MFO algorithm converges prematurely and stagnates in advance to local optima, while the IMFO algorithm almost reaches the global minimum.In summary, because the IMFO algorithm introduces Tent chaotic sequence initialization, Lévy flight, and Metropolis criterion, the initial population distribution is more uniform. At the same time, it has better ability to jump out of local optima in function optimization, and the global optimization performance is better.


### 4.2. Threat Assessment Tests

#### 4.2.1. Test Settings

According to the threat assessment framework, the input dimension of the neural network is determined to be 7, that is, threat value, target type, damage ability, interference ability, airway angle, target speed, and target distance, and the output dimension is 1, that is, threat value. Therefore, the gray neural network structure is 1-1-8-1.The IMFO algorithm parameter setting is the same as above, dimension is 7 and maximum running times of GNNM is 200. The experimental data randomly selected 100 groups from the target threat database, of which 75 groups are training sets and 25 groups are test sets. Some data are shown in Table4.Table 4
Ten sets of target threat database data.
Serial number
Threat level
Target type
Damage ability
Interference ability
Airway angle (°)
Target speed (m/s)
Target distance (km)
Threat value
1
0.63
BA
0.922
Strong
9
400
100
0.5707
2
0.34
BA
0.796
Medium
17
450
200
0.5333
3
0.35
FB
0.523
Super strong
5
600
100
0.6895
4
0.23
FB
0.674
Medium
18
550
240
0.5587
5
0.21
AH
0.646
Weak
12
110
150
0.7142
6
0.15
PGB
0.525
No
13
300
5
0.3015
7
0.25
AGM
0.884
No
11
740
40
0.5256
8
0.17
CM
0.746
No
9
690
15
0.5344
9
0.23
TBM
0.875
No
13
880
50
0.6188
10
0.13
ARM
0.627
No
16
750
8
0.6426Since the range of indicators is not uniform after quantification, the quantized data is normalized, and the data of each index is uniformly compressed into the interval [0, 1], which is beneficial to speed up the training of the model and improve the learning outcomes. The formula is as follows:(22)x′=x−xminxmax−xmin,where x represents the actual value of a variable, xmax represents the maximum value, xmin represents the minimum value, and x′ represents the normalized value.In order to reduce the error, the experiment takes the average of 30 running results and compares it with the simulation results of the GNNM model, the original MFO-optimized GNNM model, and the IMFO-optimized GNNM model.


#### 4.2.2. Accuracy Analysis

Firstly, the parameters of GNNM are optimized by the training set to obtain the optimal initialization parameters. The error evolution curve is shown in Figure16.Figure 16
Error evolution curve.As can be seen from Figure16, in the weight optimization process, IMFO has a faster optimization speed, a smaller mean square error, and better weight parameters.After the three models are trained in the training set, the test set is input to test its accuracy. The GNNM model uses random initialization weights, the MFO-GNNM model uses MFO-optimized weight parameters, and the IMFO-GNNM model uses IMFO-optimized weight parameters. The test results are shown in Figure17.Figure 17
Assessment results of test set.
(a)
(b)
(c)The expected value in the figure represents the actual threat value. It can be seen from the figure that when the unoptimized GNNM predicts the test set, the predicted value of the individual points deviates greatly from the expected value and the accuracy is not high enough. The prediction effect of the MFO-GNNM model has improved. The IMFO-GNNM model with optimized initial weights has the best weight value, and the deviation between the predicted value and the expected value is the smallest, and the prediction effect is the best.The relative error and mean square error of the three models on the test set are compared below, as shown in Figure18 and Table 5.Figure 18
Relative error of test set.Table 5
Mean square error of test set.
Models
GNNM
MFO-GNNM
IMFO-GNNM
MSE
0.013
0.004
0.0012It can be seen from Figure18 and Table 5 that the relative error of the individual points predicted by the GNNM model reaches about 0.8, the relative error is larger, and the mean square error is 0.013. The relative error of the MFO-GNNM model decreases, and the mean square error is reduced to 0.004; the IMFO-GNNM model has a small overall prediction error and a mean square error of only 0.0012. This indicates that MFO algorithm optimization weight can improve the learning effect of the model. The IMFO algorithm can obtain better initial weight parameters, so the prediction relative error is the smallest and the accuracy is the highest, which verifies the effectiveness of the proposed method.


#### 4.2.3. Real-Time Performance Analysis

The training set is used to evaluate the running time of the three algorithms. The results are shown in Figure19.Figure 19
Elapsed time of the algorithms.As can be seen from the figure, the running time of the three models: GNNM < IMFO-GNNM < MFO-GNNM. This is because the GNNM model directly inputs the training set for training after random initialization, while both the MFO-GNNM model and the IMFO-GNNM model need to be trained in the training set to obtain the optimized weight parameters, so the model runs for a long time. At the same time, the IMFO model has been improved and the search speed is faster, so the IMFO-GNNM model runs slightly faster than the MFO-GNNM model.Considering that the neural network algorithm adopts offline training and online prediction, the evaluation time of 5, 10, 15, 20, and 25 sets of data is compared with the three algorithms. The results are shown in Figure20.Figure 20
Assessment time of the algorithms.With reference to Figure20, the evaluation time of the three models is approximative. This is because the three models use the same network structure, but with different weights. Moreover, it takes less than 10 ms to evaluate 25 sets of data, which meets the real-time requirements.In summary, the IMFO-GNNM model sacrifices the training time of the algorithm while improving its accuracy, but it does not affect the real-time assessment. Therefore, the IMFO-GNNM model can improve the accuracy of air target threat assessment and has better real-time performance.


## 4.1. Verification Tests by Benchmark Functions

This paper selects six classical benchmark functions in CEC2010 to test the performance of the algorithm, as shown in Table2. Among them, the Sphere function is a smooth monotonic function with only one global minimum, which is used to test the convergence speed of the algorithm. The Rosenbrock function is a nonconvex unimodal function with multiple local extremums, which is mainly used to test the convergence speed and optimization accuracy of the algorithm; Schwefel 2.26, Rastrigin, Ackley, and Griewank functions are complex multimodal functions with multiple local minima. They are mainly used to test the global optimization performance of the algorithm and ability to jump out of local extremums. These functions can more comprehensively examine the optimization performance of the algorithm.Table 2
Benchmark functions.
Function
Equation and range
Range
Optimum
Property
Sphere
f
X
=
∑
i
=
1
D
x
i
2
[−100, 100]
0
Unimodal
Rosenbrock
f
X
=
∑
i
=
1
D
−
1
100
x
i
+
1
−
x
i
2
2
+
x
i
−
1
2
[−30, 30]
0
Unimodal
Schwefel 2.26
f
X
=
−
∑
i
=
1
D
x
i
sin
x
i
[−500, 500]
−418.98
Multimodal
Rastrigin
f
X
=
∑
i
=
1
D
100
x
i
2
−
10
cos
2
π
x
i
+
10
[−5.12, 5.12]
0
Multimodal
Ackley
f
X
=
−
20
exp
−
0.2
1
/
n
∑
i
=
1
D
x
i
2
−
exp
1
/
n
∑
i
=
1
D
cos
2
π
x
i
+
20
+
e
[−30, 30]
0
Multimodal
Griewank
f
X
=
1
/
1000
∑
i
=
1
D
x
i
2
−
∏
i
=
1
D
cos
x
i
/
i
+
1
[−600, 600]
0
MultimodalIn order to evaluate the effectiveness of the IMFO algorithm, this paper uses the above benchmark function to test and compare it with the MFO algorithm. The simulation was carried out with MATLAB R2014a. The parameters were set as follows: the number of moths is 30, the maximum number of iterations is 1000, the helix parameter is 1, the Lévy flight parameter isβ = 1.5, and the weighting coefficient of the new solution and the old solution is 0.5. In order to reduce the error, the fitness mean obtained by continuously running 30 times is taken as the test result.The fitness convergence curves of the 10-dimensional benchmark function test are shown in Figure10–15. The fitness mean and standard deviation of the 10-dimensional and 50-dimensional benchmark function tests are shown in Table 3.Figure 10
Sphere (10 dimensions).Figure 11
Rosenbrock (10 dimensions).Figure 12
Schwefel (10 dimensions).Figure 13
Rastrigin (10 dimensions).Figure 14
Ackley (10 dimensions).Figure 15
Griewank (10 dimensions).Table 3
Optimization results comparison of benchmark functions.
Function
Dimension
Algorithm
Mean
Standard deviation
Sphere
10
MFO
8.738 × 10−35
8.043 × 10−35
IMFO
6.544 × 10−55
1.936 × 10−50
50
MFO
9.965 × 10−3
4.507 × 10−2
IMFO
3.346 × 10−50
1.672 × 10−48
Rosenbrock
10
MFO
7.971
1.463
IMFO
8.954 × 10–1
3.041 × 10−1
50
MFO
4.873 × 101
1.676
IMFO
9.580 × 10−1
1.637
Schwefel 2.26
10
MFO
−0.999 × 103
1.649 × 103
IMFO
−4.188 × 103
1.403
50
MFO
−1.392 × 104
7.090 × 102
IMFO
−2.094 × 104
5.981 × 10−2
Rastrigin
10
MFO
5.597
6.085
IMFO
4.967 × 10−16
3.455 × 10−15
50
MFO
7.355 × 101
9.458
IMFO
3.678 × 10−14
5.279 × 10−11
Ackley
10
MFO
4.209 × 10−15
2.621 × 10−16
IMFO
8.419 × 10−16
3.847 × 10−17
50
MFO
7.118 × 10−2
2.937 × 10−1
IMFO
8.882 × 10−16
5.580 × 10−18
Griewank
10
MFO
1.314 × 10−4
8.359 × 10–4
IMFO
3.177 × 10−17
3.564 × 10−16
50
MFO
1.129 × 10−4
9.423 × 10–4
IMFO
1.761 × 10−17
2.965 × 10−15It can be seen from Figure10–15 that for the unimodal function, the IMFO algorithm converges faster than the MFO algorithm, and the optimization precision is higher. For the multimodal function, the MFO algorithm basically traps in local optima and the convergence speed is slow. However, the IMFO algorithm can jump out of the local optimal trap and search for the optimal solution all the time. It has good global optimization performance and faster convergence speed.It can be seen from Table3 that the mean and standard deviation of the optimization of the IMFO algorithm are smaller than the mean and standard deviation of the MFO algorithm, which proves that the IMFO algorithm has better performance and good robustness. Especially for the Schwefel function and the Rastrigin function, the MFO algorithm converges prematurely and stagnates in advance to local optima, while the IMFO algorithm almost reaches the global minimum.In summary, because the IMFO algorithm introduces Tent chaotic sequence initialization, Lévy flight, and Metropolis criterion, the initial population distribution is more uniform. At the same time, it has better ability to jump out of local optima in function optimization, and the global optimization performance is better.


## 4.2. Threat Assessment Tests

### 4.2.1. Test Settings

According to the threat assessment framework, the input dimension of the neural network is determined to be 7, that is, threat value, target type, damage ability, interference ability, airway angle, target speed, and target distance, and the output dimension is 1, that is, threat value. Therefore, the gray neural network structure is 1-1-8-1.The IMFO algorithm parameter setting is the same as above, dimension is 7 and maximum running times of GNNM is 200. The experimental data randomly selected 100 groups from the target threat database, of which 75 groups are training sets and 25 groups are test sets. Some data are shown in Table4.Table 4
Ten sets of target threat database data.
Serial number
Threat level
Target type
Damage ability
Interference ability
Airway angle (°)
Target speed (m/s)
Target distance (km)
Threat value
1
0.63
BA
0.922
Strong
9
400
100
0.5707
2
0.34
BA
0.796
Medium
17
450
200
0.5333
3
0.35
FB
0.523
Super strong
5
600
100
0.6895
4
0.23
FB
0.674
Medium
18
550
240
0.5587
5
0.21
AH
0.646
Weak
12
110
150
0.7142
6
0.15
PGB
0.525
No
13
300
5
0.3015
7
0.25
AGM
0.884
No
11
740
40
0.5256
8
0.17
CM
0.746
No
9
690
15
0.5344
9
0.23
TBM
0.875
No
13
880
50
0.6188
10
0.13
ARM
0.627
No
16
750
8
0.6426Since the range of indicators is not uniform after quantification, the quantized data is normalized, and the data of each index is uniformly compressed into the interval [0, 1], which is beneficial to speed up the training of the model and improve the learning outcomes. The formula is as follows:(22)x′=x−xminxmax−xmin,where x represents the actual value of a variable, xmax represents the maximum value, xmin represents the minimum value, and x′ represents the normalized value.In order to reduce the error, the experiment takes the average of 30 running results and compares it with the simulation results of the GNNM model, the original MFO-optimized GNNM model, and the IMFO-optimized GNNM model.


### 4.2.2. Accuracy Analysis

Firstly, the parameters of GNNM are optimized by the training set to obtain the optimal initialization parameters. The error evolution curve is shown in Figure16.Figure 16
Error evolution curve.As can be seen from Figure16, in the weight optimization process, IMFO has a faster optimization speed, a smaller mean square error, and better weight parameters.After the three models are trained in the training set, the test set is input to test its accuracy. The GNNM model uses random initialization weights, the MFO-GNNM model uses MFO-optimized weight parameters, and the IMFO-GNNM model uses IMFO-optimized weight parameters. The test results are shown in Figure17.Figure 17
Assessment results of test set.
(a)
(b)
(c)The expected value in the figure represents the actual threat value. It can be seen from the figure that when the unoptimized GNNM predicts the test set, the predicted value of the individual points deviates greatly from the expected value and the accuracy is not high enough. The prediction effect of the MFO-GNNM model has improved. The IMFO-GNNM model with optimized initial weights has the best weight value, and the deviation between the predicted value and the expected value is the smallest, and the prediction effect is the best.The relative error and mean square error of the three models on the test set are compared below, as shown in Figure18 and Table 5.Figure 18
Relative error of test set.Table 5
Mean square error of test set.
Models
GNNM
MFO-GNNM
IMFO-GNNM
MSE
0.013
0.004
0.0012It can be seen from Figure18 and Table 5 that the relative error of the individual points predicted by the GNNM model reaches about 0.8, the relative error is larger, and the mean square error is 0.013. The relative error of the MFO-GNNM model decreases, and the mean square error is reduced to 0.004; the IMFO-GNNM model has a small overall prediction error and a mean square error of only 0.0012. This indicates that MFO algorithm optimization weight can improve the learning effect of the model. The IMFO algorithm can obtain better initial weight parameters, so the prediction relative error is the smallest and the accuracy is the highest, which verifies the effectiveness of the proposed method.


### 4.2.3. Real-Time Performance Analysis

The training set is used to evaluate the running time of the three algorithms. The results are shown in Figure19.Figure 19
Elapsed time of the algorithms.As can be seen from the figure, the running time of the three models: GNNM < IMFO-GNNM < MFO-GNNM. This is because the GNNM model directly inputs the training set for training after random initialization, while both the MFO-GNNM model and the IMFO-GNNM model need to be trained in the training set to obtain the optimized weight parameters, so the model runs for a long time. At the same time, the IMFO model has been improved and the search speed is faster, so the IMFO-GNNM model runs slightly faster than the MFO-GNNM model.Considering that the neural network algorithm adopts offline training and online prediction, the evaluation time of 5, 10, 15, 20, and 25 sets of data is compared with the three algorithms. The results are shown in Figure20.Figure 20
Assessment time of the algorithms.With reference to Figure20, the evaluation time of the three models is approximative. This is because the three models use the same network structure, but with different weights. Moreover, it takes less than 10 ms to evaluate 25 sets of data, which meets the real-time requirements.In summary, the IMFO-GNNM model sacrifices the training time of the algorithm while improving its accuracy, but it does not affect the real-time assessment. Therefore, the IMFO-GNNM model can improve the accuracy of air target threat assessment and has better real-time performance.


## 4.2.1. Test Settings

According to the threat assessment framework, the input dimension of the neural network is determined to be 7, that is, threat value, target type, damage ability, interference ability, airway angle, target speed, and target distance, and the output dimension is 1, that is, threat value. Therefore, the gray neural network structure is 1-1-8-1.The IMFO algorithm parameter setting is the same as above, dimension is 7 and maximum running times of GNNM is 200. The experimental data randomly selected 100 groups from the target threat database, of which 75 groups are training sets and 25 groups are test sets. Some data are shown in Table4.Table 4
Ten sets of target threat database data.
Serial number
Threat level
Target type
Damage ability
Interference ability
Airway angle (°)
Target speed (m/s)
Target distance (km)
Threat value
1
0.63
BA
0.922
Strong
9
400
100
0.5707
2
0.34
BA
0.796
Medium
17
450
200
0.5333
3
0.35
FB
0.523
Super strong
5
600
100
0.6895
4
0.23
FB
0.674
Medium
18
550
240
0.5587
5
0.21
AH
0.646
Weak
12
110
150
0.7142
6
0.15
PGB
0.525
No
13
300
5
0.3015
7
0.25
AGM
0.884
No
11
740
40
0.5256
8
0.17
CM
0.746
No
9
690
15
0.5344
9
0.23
TBM
0.875
No
13
880
50
0.6188
10
0.13
ARM
0.627
No
16
750
8
0.6426Since the range of indicators is not uniform after quantification, the quantized data is normalized, and the data of each index is uniformly compressed into the interval [0, 1], which is beneficial to speed up the training of the model and improve the learning outcomes. The formula is as follows:(22)x′=x−xminxmax−xmin,where x represents the actual value of a variable, xmax represents the maximum value, xmin represents the minimum value, and x′ represents the normalized value.In order to reduce the error, the experiment takes the average of 30 running results and compares it with the simulation results of the GNNM model, the original MFO-optimized GNNM model, and the IMFO-optimized GNNM model.


## 4.2.2. Accuracy Analysis

Firstly, the parameters of GNNM are optimized by the training set to obtain the optimal initialization parameters. The error evolution curve is shown in Figure16.Figure 16
Error evolution curve.As can be seen from Figure16, in the weight optimization process, IMFO has a faster optimization speed, a smaller mean square error, and better weight parameters.After the three models are trained in the training set, the test set is input to test its accuracy. The GNNM model uses random initialization weights, the MFO-GNNM model uses MFO-optimized weight parameters, and the IMFO-GNNM model uses IMFO-optimized weight parameters. The test results are shown in Figure17.Figure 17
Assessment results of test set.
(a)
(b)
(c)The expected value in the figure represents the actual threat value. It can be seen from the figure that when the unoptimized GNNM predicts the test set, the predicted value of the individual points deviates greatly from the expected value and the accuracy is not high enough. The prediction effect of the MFO-GNNM model has improved. The IMFO-GNNM model with optimized initial weights has the best weight value, and the deviation between the predicted value and the expected value is the smallest, and the prediction effect is the best.The relative error and mean square error of the three models on the test set are compared below, as shown in Figure18 and Table 5.Figure 18
Relative error of test set.Table 5
Mean square error of test set.
Models
GNNM
MFO-GNNM
IMFO-GNNM
MSE
0.013
0.004
0.0012It can be seen from Figure18 and Table 5 that the relative error of the individual points predicted by the GNNM model reaches about 0.8, the relative error is larger, and the mean square error is 0.013. The relative error of the MFO-GNNM model decreases, and the mean square error is reduced to 0.004; the IMFO-GNNM model has a small overall prediction error and a mean square error of only 0.0012. This indicates that MFO algorithm optimization weight can improve the learning effect of the model. The IMFO algorithm can obtain better initial weight parameters, so the prediction relative error is the smallest and the accuracy is the highest, which verifies the effectiveness of the proposed method.


## 4.2.3. Real-Time Performance Analysis

The training set is used to evaluate the running time of the three algorithms. The results are shown in Figure19.Figure 19
Elapsed time of the algorithms.As can be seen from the figure, the running time of the three models: GNNM < IMFO-GNNM < MFO-GNNM. This is because the GNNM model directly inputs the training set for training after random initialization, while both the MFO-GNNM model and the IMFO-GNNM model need to be trained in the training set to obtain the optimized weight parameters, so the model runs for a long time. At the same time, the IMFO model has been improved and the search speed is faster, so the IMFO-GNNM model runs slightly faster than the MFO-GNNM model.Considering that the neural network algorithm adopts offline training and online prediction, the evaluation time of 5, 10, 15, 20, and 25 sets of data is compared with the three algorithms. The results are shown in Figure20.Figure 20
Assessment time of the algorithms.With reference to Figure20, the evaluation time of the three models is approximative. This is because the three models use the same network structure, but with different weights. Moreover, it takes less than 10 ms to evaluate 25 sets of data, which meets the real-time requirements.In summary, the IMFO-GNNM model sacrifices the training time of the algorithm while improving its accuracy, but it does not affect the real-time assessment. Therefore, the IMFO-GNNM model can improve the accuracy of air target threat assessment and has better real-time performance.


## 5. Conclusions

(1)
By analyzing the actual situation of air defense operations and considering threat level, threat ability, and threat extent, a reasonable air target threat assessment framework is established, and each indicator is scientifically quantified.
(2)
The IMFO algorithm is proposed. By introducing Tent chaos, Lévy flight, and Metropolis criterion, the global optimization performance of the algorithm is improved.
(3)
Combining powerful learning ability of GNNM and excellent optimization performance of IMFO, a new air target threat assessment model based on IMFO-GNNM is established, which has high accuracy and good real-time performance for threat assessment of incoming targets.We will continue to study and propose better threat assessment models in the future. At the same time, we will study how to incorporate the proposed threat assessment model into the air defense system.



---

*Source: 4203538-2019-10-10.xml*