# Solving the Multisensor Resource Scheduling Problem for Missile Early Warning by a Hybrid Discrete Artificial Bee Colony Algorithm

Authors: Wei Liu; Changyun Liu; Xiangke Guo; Sheng He; Liangyou Fan
Journal: Journal of Sensors
 (2022)

---

## Abstract

Aiming at the problem of multisensor resource scheduling in missile early warning operation, a scheduling decomposition strategy for missile early warning tasks under cooperative detection is proposed. Taking the detection benefit factor, target threat factor, and handover factor as the fitness function, we establish a sensor-subtask assignment (SSA) model and propose a hybrid discrete artificial bee colony (HDABC) algorithm to solve the optimal solution of the SSA model. The HDABC algorithm has the following improvements: in the initialization stage, a sensor-subtask-based coding method is designed to reduce the solution dimension, and the heuristic rules are used to obtain excellent populations to improve the convergence speed; in the employed bee and onlooker bee stage, a food source update strategy based on discrete differential mutation (DDM) operation is proposed to improve the searchability of the algorithm, and a sorting-based adaptive probability (SAP) selection method is applied to enhance the global search and local optimization capacities. Simulation experiments were carried out in operation scenarios of different scales. Experimental results showed that the proposed HDABC algorithm can obtain the optimal scheduling schemes and had a better solving performance when solving the SSA model, especially in the medium-scale and large-scale operation scenarios.


---

## Body

## 1. Introduction

Missile early warning resource scheduling refers to dynamically determining the multisensor detection and tracking sequences of multitarget under the condition of limited sensor resources and then determining the time for tracking and resource assignment, so as to achieve continuous and stable detection and tracking of threat targets. Its essence is a kind of nonlinear combinatorial optimization decision-making problem for multisensor detection of multitarget. At present, most of the researches are aimed at single sensor resource scheduling problems, focusing on scheduling methods and algorithm optimization problems under the constraints of time, energy, and computing resources [1–4]. However, with the accelerated construction of missile early warning systems in the future, missile early warning operations will be characterized by multisensor cooperative tracking and detection. The research on multisensor detection of multitarget resource scheduling problems under the condition of resource conflict has become an urgent problem to be solved.Task priority determination and scheduling algorithm design are the main problems of multisensor resource scheduling. Therefore, the multisensor scheduling solution methods can be divided into two categories.The first category refers to the optimization of task priority determination. Task priority determines the order of resource invocation and reflects the importance of tasks. Setting task priorities ensures that important tasks will not be lost during scheduling, which requires that the determination of priority has good adaptability to environmental changes. The traditional task priority determination method adopts the past operation experience for fixed settings, but this method is not flexible enough to effectively schedule the new tasks in the scene. After that, the priority determination method is proposed, such as highest priority first (HPF) algorithm [5], earliest deadline first (EDF) algorithm [6], and its improved algorithm: modified earliest deadline first (MEDF) algorithm, highest priority, earliest deadline first (HPEDF) algorithm [7–10], etc. However, the above algorithm does not make full use of the prior information of the target in the process of priority planning, and there is a problem of too strong subjectivity caused by artificially assigning the priority of the task mode. To enhance the reliability and accuracy of resource scheduling, the multiparameter synthesis priority determination algorithm is proposed, which involves the optimal task benefit factor and target threat factor. Cheng et al. [11] established a scheduling model considering time and energy constraints from the perspective of scheduling benefits; Chen et al. [12] proposed a heuristic multibeam dwell scheduling algorithm based on maximal scheduling benefits; Zhang et al. [13] proposed an algorithm to calculate the synthesis priority of the task by combining the threat density of the target and the deadline of the task. The simulation results show the significant improvement of the comprehensive multiparameter method compared with the traditional task priority determination method. On this basis, this paper will establish a multisensor task priority determination method that is more suitable for missile early warning.The second category refers to the design of the scheduling algorithm, such as the algorithm based on the auction mechanism [14, 15], game-theoretic framework [16], and approximate dynamic programming [17, 18]. These algorithms can effectively solve the problem with smaller dimensions, but difficult to solve the problem with higher dimensions. Artificial intelligence algorithm is a popular optimization algorithm to solve the problem of target allocation and resource scheduling, such as ant colony optimization algorithm (ACO) [19], genetic algorithm (GA) [20], particle swarm optimization algorithm (PSO) [21, 22], and artificial bee colony algorithm (ABC) [23]. Through a lot of experiments [24, 25], it is proved that the ABC algorithm has better optimization ability than other intelligent optimization algorithms and is not easy to fall into local optima. At present, the ABC algorithm has been applied to many resource scheduling problems. For example, aiming at the problems of slow convergence speed and low search efficiency of weapon resource scheduling algorithm, Chang et al. [26] proposed an improved ABC algorithm, which adopted rule-based heuristic factors for initialization and improved the convergence speed and accuracy of the algorithm. Pang et al. [27] proposed an improved ABC algorithm based on double probability to obtain the sensor management scheme based on the fitness function of target detecting risk and target tracking risk. Xia et al. [28] proposed an ABC algorithm based on a jamming resource scheduling problem with few parameter adjustments, and the proposed ABC algorithm has better performance in convergence speed and accuracy. To sum up, the ABC algorithm is selected in this paper to optimize the multisensor resource scheduling problem, further supporting the better application of the ABC algorithm in the field of resource scheduling.Although the above methods are diverse, they still have the following shortcomings: (1) most task priority methods are aimed at aerodynamic targets, which do not conform to the characteristics of time-sensitive target early warning, such as ballistic missiles. (2) The convergence and accuracy of the ABC algorithm are not good when dealing with multisensor detection of multitarget problems. (3) The sensor-target assignment scheme cannot match the capabilities of the sensors well to obtain the optimal task benefit.This study is aimed at building a multisensor resource scheduling decision model to optimize the task priority strategy and improving the scheduling algorithm to improve convergence and accuracy. The main contributions of this article can be summarized as follows:(1)
Under the premise of predictable trajectory, a missile early warning task decomposition strategy based on periodic scheduling and task decomposition is adopted to transform the multisensor resource scheduling problem into a sensor-subtask assignment (SSA) optimization problem. Taking the cooperative detection of P-band ground-based early warning radar (PBR) and X-band ground-based early warning radar (XBR) as an example, a multisensor resource scheduling decision model based on target threat and detection benefit is constructed(2)
A hybrid discrete artificial bee colony (HDABC) algorithm is proposed to solve the resource assignment problem of multisensor cooperative detection. The algorithm is improved from the aspects of coding rules, heuristic initialization strategy, food source update strategy, and food source selection probability, so that the improved algorithm could be more excellent in dealing with such problems


## 2. Periodic Scheduling and Task Decomposition Strategy

The periodic scheduling and task decomposition of missile early warning tasks solve the problems of early warning resource time planning and target grouping, that is, when to generate the scheduling scheme and how to assign sensors and targets [29]. Figure 1 illustrates the missile early warning resource scheduling framework based on periodic scheduling and task decomposition.Figure 1
Missile early warning resource scheduling process framework based on periodic scheduling and task decomposition.


### 2.1. Periodic Scheduling

The duration of the scheduling period greatly influences the scheduling effect [30]: if the generated scheduling period is too long, with the increase in the tracking error of target detection, the scheduling scheme will not meet the detection reality; if the scheduling period is too short, the workload of the solution will be significantly increased. Therefore, the scheduling period should be dynamically adjusted according to the measurement results of the target and the changing trend of the task, which mainly depends on the accuracy of the prediction information of the early warning system and the complexity of the battlefield space. The higher the accuracy of the forecast and the smaller the scale of operations, the longer the required scheduling period and the higher the reliability, and vice versa; the accuracy of the scheduling plan will be affected. In addition, in emergencies such as the emergence of new targets, the addition and withdrawal of early warning resources, and the target deviation from the predicted ballistics, the period must be dynamically calculated to ensure self-adaptation to complex battlefield tasks.The method for calculating the duration of the scheduling period is as follows: assuming that the prediction accuracy of the early warning system to the target at the timet0 is Pt, the number of targets to be scheduled is Nbm, the number of sensors is Ns, the average time of scheduling scheme generation is Cal, the frequency of generating new targets in a scheduling period is Ftg, and then, the duration of the next scheduling period is as follows:
(1)T=fstPt,Cal,Ftg,Nst=t0⇒ST0=t0,t0+T,where fst· is the calculation function of the periodic scheduling length.


### 2.2. Task Decomposition Strategy

Task decomposition refers to refining the complex visible relationship between sensors and targets in a scheduling period into a sequence of subtasks that can be directly executed and completed by sensors [21]. Different task decomposition strategies often lead to different subtask sequences, which eventually lead to different scheduling schemes. At present, the most commonly used task decomposition strategies are the “longest observation time” and “start and end time division” methods [22], but they are all for the decomposition of single-target detection tasks. When the number of targets existing at the same time and in the same space is too large, the complexity of the decomposition of the aforementioned strategy increases greatly. Therefore, we adopt the task decomposition method of “minimum scheduling interval” to avoid the problem that the amount of calculation increases significantly due to excessive decomposition. Proceed as follows:Step 1.
Calculate the visibility of each sensor to the target in a scheduling period and set the minimum scheduling timeTsub according to the predicted trajectory of all targets.Step 2.
Divide the predicted trajectory based on the visibility to generatek subtasks. The time of each subtask is Tii∈k.Step 3.
If there isTi≥Tsub, the output is a subtask STj=Ti; if there is Ti<Tsub, then when Ti+Ti+1+⋯+Ti+n≥Tsub, the output is a subtask STj=∑ni=0nTi+ni.Step 4.
Until∀STj≥Tsub, output the subtask sequence; otherwise, jump to Step 3.After the aforementioned periodic scheduling and task decomposition of the missile early warning task, each early warning resource will correspond to a series of subtask sequences. At this time, the scheduling problem of early warning resources is transformed into a “sensor-subtask assignment” (SSA) problem based on the scheduling period. The scheduling scheme is periodically generated to determine which subtasks are to be detected and which early warning resources are assigned for execution, thereby greatly reducing the complex correspondence between tasks and resources.


## 2.1. Periodic Scheduling

The duration of the scheduling period greatly influences the scheduling effect [30]: if the generated scheduling period is too long, with the increase in the tracking error of target detection, the scheduling scheme will not meet the detection reality; if the scheduling period is too short, the workload of the solution will be significantly increased. Therefore, the scheduling period should be dynamically adjusted according to the measurement results of the target and the changing trend of the task, which mainly depends on the accuracy of the prediction information of the early warning system and the complexity of the battlefield space. The higher the accuracy of the forecast and the smaller the scale of operations, the longer the required scheduling period and the higher the reliability, and vice versa; the accuracy of the scheduling plan will be affected. In addition, in emergencies such as the emergence of new targets, the addition and withdrawal of early warning resources, and the target deviation from the predicted ballistics, the period must be dynamically calculated to ensure self-adaptation to complex battlefield tasks.The method for calculating the duration of the scheduling period is as follows: assuming that the prediction accuracy of the early warning system to the target at the timet0 is Pt, the number of targets to be scheduled is Nbm, the number of sensors is Ns, the average time of scheduling scheme generation is Cal, the frequency of generating new targets in a scheduling period is Ftg, and then, the duration of the next scheduling period is as follows:
(1)T=fstPt,Cal,Ftg,Nst=t0⇒ST0=t0,t0+T,where fst· is the calculation function of the periodic scheduling length.


## 2.2. Task Decomposition Strategy

Task decomposition refers to refining the complex visible relationship between sensors and targets in a scheduling period into a sequence of subtasks that can be directly executed and completed by sensors [21]. Different task decomposition strategies often lead to different subtask sequences, which eventually lead to different scheduling schemes. At present, the most commonly used task decomposition strategies are the “longest observation time” and “start and end time division” methods [22], but they are all for the decomposition of single-target detection tasks. When the number of targets existing at the same time and in the same space is too large, the complexity of the decomposition of the aforementioned strategy increases greatly. Therefore, we adopt the task decomposition method of “minimum scheduling interval” to avoid the problem that the amount of calculation increases significantly due to excessive decomposition. Proceed as follows:Step 1.
Calculate the visibility of each sensor to the target in a scheduling period and set the minimum scheduling timeTsub according to the predicted trajectory of all targets.Step 2.
Divide the predicted trajectory based on the visibility to generatek subtasks. The time of each subtask is Tii∈k.Step 3.
If there isTi≥Tsub, the output is a subtask STj=Ti; if there is Ti<Tsub, then when Ti+Ti+1+⋯+Ti+n≥Tsub, the output is a subtask STj=∑ni=0nTi+ni.Step 4.
Until∀STj≥Tsub, output the subtask sequence; otherwise, jump to Step 3.After the aforementioned periodic scheduling and task decomposition of the missile early warning task, each early warning resource will correspond to a series of subtask sequences. At this time, the scheduling problem of early warning resources is transformed into a “sensor-subtask assignment” (SSA) problem based on the scheduling period. The scheduling scheme is periodically generated to determine which subtasks are to be detected and which early warning resources are assigned for execution, thereby greatly reducing the complex correspondence between tasks and resources.


## 3. Resource Scheduling Model Based on Sensor-Subtask Assignment

When sensors have a visible relationship to the same target, in order to analyze the rationality of target detection by sensors, a resource scheduling model based on sensor-subtask assignment (SSA) is used to describe it. The SSA model is established as follows.


### 3.1. Fitness Function

The fitness function is composed of factors such as detection benefit factor, target threat factor, and target handover factor, which are expressed as follows:(2)maxJX=JX1+JX2+⋯+JXk,where the multisensor cooperative detection benefit JX is calculated as
(3)JXk=∑i=1m∑j=1nai,jk·Beni,jk·Thri,j·Hank,where Beni,jk is the detection benefit factor, Thri,j is the target threat factor, Hank is the target handover factor, and ai,jk is a decision variable, describing the detection of the i-th target’s j-th subtask by the k-th sensor, and the calculation formula is
(4)ai,jk=1,0,thek−thsensordetectsthej−thsubtaskofthei−thtarget,otherwise.


#### 3.1.1. Detection Benefit Factor

The detection benefit factor is expressed as(5)Beni,jk=l1·Disi,jk+l2·Losi,jk+l3·Prii,jk,where l1, l2, and l3 are the weight factors and l1+l2+l3=1.(1) Spatial DistanceDisi,jk.Disi,jk is to indicate the influence of the distance between the target and the sensor on the detection effect. The closer the target distance is, the better it is to improve the precision of the target tracking information, expressed as
(6)Disi,jk=Dmaxk−Di,jkDmaxk,where Dmaxk is the maximum detection range of the k-th sensor and Di,jk is the average distance of the k-th sensor for each subtask.(2) Line-of-Sight AngleLoai,jk. Loai,jk is to measure the deviation between the line-of-sight angle (LOA) of the sensor to detect the target and the optimal angle. Loai,jk affects the sensor’s ability to identify the target. The further away from the optimal angle, the worse the detection effect, expressed as
(7)Loai,j=1−αi,jk−αibestkαibestk,where αi_bestk is the optimal LOA of the k-th sensor to the i-th target and αi,jk is the LOA for each subtask of the i-th sensor.(3) Detection CoverageArci,jk. Arci,jk is to indicate the detection coverage capability of the sensor to the subtask. The longer the ballistic arc length is covered by the sensor’s detection of the subtask, the better it can avoid the target loss, expressed as
(8)Arci,jk=li,jkli_missile,where li_missile is the total predicted arc length for the i-th target and li,jk is the detection arc length of the i-th sensor for each subtask.(4) Detection PriorityPrik. Prik is to indicate the detection priority of the sensor. The higher the identification accuracy of the sensor, the higher its priority, expressed as
(9)Prik=1,whensensorkisPBR,2,whensensorkisXBR.


#### 3.1.2. Target Threat Factor

In missile early warning operations, ballistic targets usually appear as cluster targets and perform saturated strikes at strategic positions. Therefore, under the condition of limited resources, it is necessary to conduct a threat assessment on all targets to distinguish the detection priority and realize the reasonable assignment of resources. The threat factor of thei-th target is expressed as
(10)Thri=Ideiη1·Clai+η2·Fini,where η1 and η2 are the weight factors and η1+η2=1.(1) Friend or Foe InformationIdei. Idei is to indicate the friend or foe information of the target. Targets identified as ours are not detected.
(11)Idei=0,Identifiedasourtarget,1,Identifyasenemytarget.(2) Target Category InformationClai. Clai includes information such as warhead target confidence and target type and indicates the threat level of the target by matching the missile model and type. Clai is determined by matching in the feature database after the target is comprehensively identified. The greater the value, the more likely the target is to be a warhead target, and the greater the level of threat. The value range is
(12)Clai∈0,1.(3) Flight DataFini. Fini includes information such as the predicted range of each target, the speed of the shutdown point, the remaining flight time, the ballistic inclination, and the range between the striking point and the defense point, expressed as
(13)Fini=∑μi,fi·ωi,fi,where ωi,fi is the threat value corresponding to the fi-th factor, which is obtained from the prediction information of the early warning system. The value range is (0, 1), and the larger the value, the higher the threat degree; μi,t is the weight of each factor and ∑μi,fi=1.


#### 3.1.3. Target Handover Factor

The target handover factorHank is used to express the influence of the number of target handovers on the cooperative of sensors. To avoid mistaking and losing track in the process of tracking the target, the number of target handovers should be minimized. In addition, Hank is an important factor to control the stability of the solution. The calculation formula is
(14)Han=1−∑i=1m∑j=1n−1ai,jk⊗ai,j+1kCmaxk,where aijk⊗ai,j+1k is an XOR operation, indicating whether the k-th sensor switches the tracking target and Cmaxk is the maximum trackable target capacity of the k-th sensor.


### 3.2. Constraints

(1)
From the overall perspective of missile defense and resource optimization, it is necessary to ensure that each ballistic target can be detected as much as possible and that subtasks with detection overlap only occupy one sensor, that is(15)∑k=1sai,jk≤1,if∑k=1sBeni,jk>0(2)
From the limitation of tracking capacity, the number of targets tracked by the sensor cannot be higher than its target capacity, and a certain redundancy should be reserved to avoid being unable to respond to emergency due to resource overload, so it is expressed as(16)∑i=1mai,jk≤Cmaxk,where Cmaxk is the maximum trackable target capacity of the k-th sensor


### 3.3. Generation of Scheduling Scheme

Solving the SSA model can obtain the decision matricesAi of all subtasks detected by the sensor. By combining the decision matrix according to the sensor number, a scheduling scheme S=A1,A2,⋯,As within a scheduling period can be generated.
(17)S=a1,11⋯a1,n1⋯ai,j1am,11⋯am,n1⋯a1,1s⋯ai,ns⋯⋯am,ns⋯⋯am,1s⋯am,ns.


## 3.1. Fitness Function

The fitness function is composed of factors such as detection benefit factor, target threat factor, and target handover factor, which are expressed as follows:(2)maxJX=JX1+JX2+⋯+JXk,where the multisensor cooperative detection benefit JX is calculated as
(3)JXk=∑i=1m∑j=1nai,jk·Beni,jk·Thri,j·Hank,where Beni,jk is the detection benefit factor, Thri,j is the target threat factor, Hank is the target handover factor, and ai,jk is a decision variable, describing the detection of the i-th target’s j-th subtask by the k-th sensor, and the calculation formula is
(4)ai,jk=1,0,thek−thsensordetectsthej−thsubtaskofthei−thtarget,otherwise.


### 3.1.1. Detection Benefit Factor

The detection benefit factor is expressed as(5)Beni,jk=l1·Disi,jk+l2·Losi,jk+l3·Prii,jk,where l1, l2, and l3 are the weight factors and l1+l2+l3=1.(1) Spatial DistanceDisi,jk.Disi,jk is to indicate the influence of the distance between the target and the sensor on the detection effect. The closer the target distance is, the better it is to improve the precision of the target tracking information, expressed as
(6)Disi,jk=Dmaxk−Di,jkDmaxk,where Dmaxk is the maximum detection range of the k-th sensor and Di,jk is the average distance of the k-th sensor for each subtask.(2) Line-of-Sight AngleLoai,jk. Loai,jk is to measure the deviation between the line-of-sight angle (LOA) of the sensor to detect the target and the optimal angle. Loai,jk affects the sensor’s ability to identify the target. The further away from the optimal angle, the worse the detection effect, expressed as
(7)Loai,j=1−αi,jk−αibestkαibestk,where αi_bestk is the optimal LOA of the k-th sensor to the i-th target and αi,jk is the LOA for each subtask of the i-th sensor.(3) Detection CoverageArci,jk. Arci,jk is to indicate the detection coverage capability of the sensor to the subtask. The longer the ballistic arc length is covered by the sensor’s detection of the subtask, the better it can avoid the target loss, expressed as
(8)Arci,jk=li,jkli_missile,where li_missile is the total predicted arc length for the i-th target and li,jk is the detection arc length of the i-th sensor for each subtask.(4) Detection PriorityPrik. Prik is to indicate the detection priority of the sensor. The higher the identification accuracy of the sensor, the higher its priority, expressed as
(9)Prik=1,whensensorkisPBR,2,whensensorkisXBR.


### 3.1.2. Target Threat Factor

In missile early warning operations, ballistic targets usually appear as cluster targets and perform saturated strikes at strategic positions. Therefore, under the condition of limited resources, it is necessary to conduct a threat assessment on all targets to distinguish the detection priority and realize the reasonable assignment of resources. The threat factor of thei-th target is expressed as
(10)Thri=Ideiη1·Clai+η2·Fini,where η1 and η2 are the weight factors and η1+η2=1.(1) Friend or Foe InformationIdei. Idei is to indicate the friend or foe information of the target. Targets identified as ours are not detected.
(11)Idei=0,Identifiedasourtarget,1,Identifyasenemytarget.(2) Target Category InformationClai. Clai includes information such as warhead target confidence and target type and indicates the threat level of the target by matching the missile model and type. Clai is determined by matching in the feature database after the target is comprehensively identified. The greater the value, the more likely the target is to be a warhead target, and the greater the level of threat. The value range is
(12)Clai∈0,1.(3) Flight DataFini. Fini includes information such as the predicted range of each target, the speed of the shutdown point, the remaining flight time, the ballistic inclination, and the range between the striking point and the defense point, expressed as
(13)Fini=∑μi,fi·ωi,fi,where ωi,fi is the threat value corresponding to the fi-th factor, which is obtained from the prediction information of the early warning system. The value range is (0, 1), and the larger the value, the higher the threat degree; μi,t is the weight of each factor and ∑μi,fi=1.


### 3.1.3. Target Handover Factor

The target handover factorHank is used to express the influence of the number of target handovers on the cooperative of sensors. To avoid mistaking and losing track in the process of tracking the target, the number of target handovers should be minimized. In addition, Hank is an important factor to control the stability of the solution. The calculation formula is
(14)Han=1−∑i=1m∑j=1n−1ai,jk⊗ai,j+1kCmaxk,where aijk⊗ai,j+1k is an XOR operation, indicating whether the k-th sensor switches the tracking target and Cmaxk is the maximum trackable target capacity of the k-th sensor.


## 3.1.1. Detection Benefit Factor

The detection benefit factor is expressed as(5)Beni,jk=l1·Disi,jk+l2·Losi,jk+l3·Prii,jk,where l1, l2, and l3 are the weight factors and l1+l2+l3=1.(1) Spatial DistanceDisi,jk.Disi,jk is to indicate the influence of the distance between the target and the sensor on the detection effect. The closer the target distance is, the better it is to improve the precision of the target tracking information, expressed as
(6)Disi,jk=Dmaxk−Di,jkDmaxk,where Dmaxk is the maximum detection range of the k-th sensor and Di,jk is the average distance of the k-th sensor for each subtask.(2) Line-of-Sight AngleLoai,jk. Loai,jk is to measure the deviation between the line-of-sight angle (LOA) of the sensor to detect the target and the optimal angle. Loai,jk affects the sensor’s ability to identify the target. The further away from the optimal angle, the worse the detection effect, expressed as
(7)Loai,j=1−αi,jk−αibestkαibestk,where αi_bestk is the optimal LOA of the k-th sensor to the i-th target and αi,jk is the LOA for each subtask of the i-th sensor.(3) Detection CoverageArci,jk. Arci,jk is to indicate the detection coverage capability of the sensor to the subtask. The longer the ballistic arc length is covered by the sensor’s detection of the subtask, the better it can avoid the target loss, expressed as
(8)Arci,jk=li,jkli_missile,where li_missile is the total predicted arc length for the i-th target and li,jk is the detection arc length of the i-th sensor for each subtask.(4) Detection PriorityPrik. Prik is to indicate the detection priority of the sensor. The higher the identification accuracy of the sensor, the higher its priority, expressed as
(9)Prik=1,whensensorkisPBR,2,whensensorkisXBR.


## 3.1.2. Target Threat Factor

In missile early warning operations, ballistic targets usually appear as cluster targets and perform saturated strikes at strategic positions. Therefore, under the condition of limited resources, it is necessary to conduct a threat assessment on all targets to distinguish the detection priority and realize the reasonable assignment of resources. The threat factor of thei-th target is expressed as
(10)Thri=Ideiη1·Clai+η2·Fini,where η1 and η2 are the weight factors and η1+η2=1.(1) Friend or Foe InformationIdei. Idei is to indicate the friend or foe information of the target. Targets identified as ours are not detected.
(11)Idei=0,Identifiedasourtarget,1,Identifyasenemytarget.(2) Target Category InformationClai. Clai includes information such as warhead target confidence and target type and indicates the threat level of the target by matching the missile model and type. Clai is determined by matching in the feature database after the target is comprehensively identified. The greater the value, the more likely the target is to be a warhead target, and the greater the level of threat. The value range is
(12)Clai∈0,1.(3) Flight DataFini. Fini includes information such as the predicted range of each target, the speed of the shutdown point, the remaining flight time, the ballistic inclination, and the range between the striking point and the defense point, expressed as
(13)Fini=∑μi,fi·ωi,fi,where ωi,fi is the threat value corresponding to the fi-th factor, which is obtained from the prediction information of the early warning system. The value range is (0, 1), and the larger the value, the higher the threat degree; μi,t is the weight of each factor and ∑μi,fi=1.


## 3.1.3. Target Handover Factor

The target handover factorHank is used to express the influence of the number of target handovers on the cooperative of sensors. To avoid mistaking and losing track in the process of tracking the target, the number of target handovers should be minimized. In addition, Hank is an important factor to control the stability of the solution. The calculation formula is
(14)Han=1−∑i=1m∑j=1n−1ai,jk⊗ai,j+1kCmaxk,where aijk⊗ai,j+1k is an XOR operation, indicating whether the k-th sensor switches the tracking target and Cmaxk is the maximum trackable target capacity of the k-th sensor.


## 3.2. Constraints

(1)
From the overall perspective of missile defense and resource optimization, it is necessary to ensure that each ballistic target can be detected as much as possible and that subtasks with detection overlap only occupy one sensor, that is(15)∑k=1sai,jk≤1,if∑k=1sBeni,jk>0(2)
From the limitation of tracking capacity, the number of targets tracked by the sensor cannot be higher than its target capacity, and a certain redundancy should be reserved to avoid being unable to respond to emergency due to resource overload, so it is expressed as(16)∑i=1mai,jk≤Cmaxk,where Cmaxk is the maximum trackable target capacity of the k-th sensor


## 3.3. Generation of Scheduling Scheme

Solving the SSA model can obtain the decision matricesAi of all subtasks detected by the sensor. By combining the decision matrix according to the sensor number, a scheduling scheme S=A1,A2,⋯,As within a scheduling period can be generated.
(17)S=a1,11⋯a1,n1⋯ai,j1am,11⋯am,n1⋯a1,1s⋯ai,ns⋯⋯am,ns⋯⋯am,1s⋯am,ns.


## 4. Hybrid Discrete Artificial Bee Colony Algorithm

The proposed SSA model is a typical nonlinear combinatorial optimization problem, and there are plenty of algorithms to solve such problems. The artificial intelligence algorithm is one of the effective ways to solve this NP-Hard problem, which can get a satisfactory solution within a given period after iterations and optimum search, but has some disadvantages, such as slow convergence speed, low efficiency, and instability solution [26]. The artificial bee colony (ABC) algorithm is an efficient artificial intelligence algorithm by simulates honeybees’ foraging behavior [31, 32], which has been applied in many fields, such as dynamic clustering [33], shortest path problem [34], and traveling salesman problem [35]. Compared with PSO, DE, and EA [36], the ABC algorithm has the characteristics of flexible structure, fewer control parameters, strong optimization ability, and great advantages in large-scale solution problems [37–39]. Considering the efficient optimization ability of the ABC algorithm, we adopt the ABC algorithm to solve the SSA problem in the study.However, the original ABC algorithm was first developed to solve continuous problems [40] and cannot be directly applied to solve problems with discrete variables, such as the SSA model. Like other artificial intelligence algorithms, the ABC algorithm has the disadvantages of weak local convergence ability and slow convergence speed [41]. In order to solve the above problems, many studies have proposed the discrete ABC (DABC) algorithm [42, 43] and improved on it. Up to now, this improved DABC algorithm has been successfully applied in combinatorial optimization problem. For example, Masdari et al. [44] proposed the chaotic discrete ABC to solve discrete problems such as clustering of sensor nodes in the wireless sensor networks; Li et al. [45] presented a sorting-based discrete artificial bee colony algorithm to solve the flow shop scheduling problem; He et al. [46] proposed a multitask bee colony band selection algorithm with variable-size clustering to solve the multitask optimization problem in band selection.Based on the above research, we propose the hybrid discrete artificial bee colony (HDABC) algorithm to solve the SSA model. We first redefine the integer coding strategy and then improve the initialization rules, food source update strategy, and food source selection probability to improve the ABC algorithm. The solution process of missile early warning resource scheduling based on HDABC is shown in Figure2.Figure 2
The solution process of missile early warning resource scheduling based on HDABC.


### 4.1. ABC Algorithm

The ABC algorithm divides bee colony into three categories: employed bee, onlooker bee, and scout bee. The goal of the bee colony is to find the optimal food source, and the food source represents all possible solutions in the solution space and is measured by fitness value. Employed bee focuses on food source detection. Onlooker bee receives food source information shared by other bees and is responsible for mining food sources. Scout bee searches for new food sources randomly when food sources are abandoned. The algorithm process is as follows:(1)
Initialization Stage. In a D-dimensional search space, the population number is NP, and the position of each food source after the t-th iteration is(18)Xit=xi,1t,xi,2t,⋯,xi,Dt,where i=1,2,⋯,NP. The fitness value of food source is fiti=fitXit(2)
Employed Bee Stage. Each employed bee chooses a food source randomly and then generates a new food source vid by the food source update strategy. The fitness fiti′ value of the new food source is calculated, and the food source with better fitness value is retained by the greedy strategy(3)
Onlooker Bee Stage. After receiving the information of the employed bee, the onlooker bee selected several food sources according to the food source selection strategy, searched for new food sources according to the food source update strategy for fitness comparison, and selected the food sources with better fitness for retention(4)
Scout Bee Stage. If the food source did not improve after lim iterations, the food source is abandoned and recorded in the tabu list. At the same time, the bees corresponding to the food source turned into scout bee and generated a new food source according to the initialization strategy randomly(5)
Repeat steps (2) to (4) until the termination condition is met and the optimal food source location is output


### 4.2. Improvements to the Initialization Stage

#### 4.2.1. Discrete Integer Coding Strategy Based on Sensor-Subtask Sequence

The decision variableai,jk in the SSA model corresponds to the bee colony individual of the algorithm. For the problem of s sensors, m targets, and n subtasks, if 0-1 coding is used, a D-dimensional (D=s×m×n) vector will be generated, which will cause a dimensional disaster as scene complexity increases. To reduce the computational complexity of the algorithm, a discrete integer coding method based on sensor-subtask sequence is proposed as shown in Figure 3.Figure 3
Discrete integer coding strategy based on sensor-subtask sequence.In this method, the decision matrixAi is encoded as the target number corresponding to the sensor, and the position without a visual relationship is filled with 0. At this time, the coding length is determined by the maximum tracking capability Cmax of each sensor; that is, D=∑k=1sCmaxk. Through this coding method, the dimension of the algorithm is effectively reduced, and the number of targets assigned by each sensor does not exceed the target capacity, which is convenient for directly expressing the scheduling scheme.


#### 4.2.2. Heuristic Initialization Rules

The heuristic initialization rules are aimed at generating an initial feasible solution to improve the quality of the initial solution and speed up the convergence [18]. According to the characteristics of SSA model, we propose the initialization rules based on target-threat-priority and resource-balance-priority.(1) Heuristic Initialization Process Based on Target-Threat-Priority.Step 1.
To select the target with the maximum threat degree, assign it to the sensorS1 which has the highest detection priority, and remove the target.Step 2.
Continue to assign the target with the maximum threat degree to the sensorS1 and remove the target until S1 reaches the maximum capacity.Step 3.
According to the methods of Steps1 and 2, assign targets to S2-Sk in turn until there is no target or the sensor resources are saturated.Step 4.
To take the assignment scheme as a heuristic initial individual to replace the initial solution.(2) Heuristic Initialization Process Based on Resource-Balance-Priority.Step 1.
To select the target with the largest detection benefit forS1, assign it to S1, and remove the target; continue to take the target with the largest detection benefit for S2-Sk and assign the target to each sensor, and remove the target. So far, each sensor has been assigned a target with the largest detection benefit.Step 2.
Repeat Step1 until there is no target or the sensor resources are saturated.Step 3.
To take the assignment scheme as a heuristic initial individual to replace the initial solution.


### 4.3. Improvements to the Food Source Update Strategy

Since the adoption of the integer coding strategy, the traditional strategy of updating food sources in the employed bee stage and the scout bee stage [47] is not applicable. Based on Zhang et al.’s research [48], a food source update strategy based on discrete differential mutation (DDM) operation is proposed to enhance the algorithm’s search capacity. The food source update formula is as follows:
(19)Vit=φ1Xit⊕φ2⊗Xit−Xkt,where the scale factors ϕ1 and ϕ2 are random numbers in [0, 1], Xit,Xkt,i,k∈1,NP is the i-th food source and the k-th food source, respectively, and Vit is the new food resource. The operation process is divided into the following three parts:Step 1.
Calculate the part ofXit−Xkt.
The operation is defined as(20)Δi=Xit−Xkt⇔δid=xkd,xid≠xkd,randiN,xid=xkd=0,0,xid=xkd≠0,,where d=1,2,⋯,D and Δi=δi1,δi2,⋯,δiD. Through this operation, the codes of the i-th food source and the k-th food source are compared bit by bit. If the code of d-th bit is nonzero and the same, it is 0; if the bit is all zero, it will take a random integer in 1,N; if the bit is different, it will take the bit encoding of the i-th food source.Step 2.
Calculate the part ofϕ2⊗Δi.
The operation is defined as(21)Ρi=φ2⊗Δi⇔ρid=roundδd·DN,φ2>rand,0,otherwise,where Ρi=ρi1,ρi2,⋯,ρiD is the mutation operator and round· is the rounding function. After this operation, the code of a certain bit is probabilistically converted into an integer in 1,D.Step 3.
Calculate the part ofφ1Xit⊕Ρi.
The operation is defined as(22)φ1Xit⊕Ρ=Vi⇔vid=swapxi,ρd,xid,φ1≤rand,insertxi,ρd,xd,φ1>rand,ρd≠0,xid,ρd=0,where swap· is the function to swap two numbers and insert· is the function that inserts the first number before the second number. The process of mutation can be illustrated as follows: if ρd=0, remain xid unchanged; if ρd≠0 and φ1≤rand, then swap the ρd-th and d-th bits of Xit; if ρd≠0 and φ1>rand, insert the ρd-th bit of Xit into the d-th bit, and move the rest of the bits backward in turn, as shown in Figure 4.Figure 4
Food source update operation based on DDM.An illustration is given for the operation of Step3. After the comparison between φ1 and rand, ρ2=8 corresponds to swap xi2 and xi8; ρ4=3 corresponds to insert xi3 into the 4-th bit of Xit, and the other bits are moved backward in turn; ρ7=5 corresponds to swap xi7 and xi5; ρ9=6 corresponds to insert xi6 into the 9-th bit of Xit, and the other bits are moved backward in turn; ρ11=5 corresponds to insert xi11 into the 5-th bit of Xit, and the other bits are moved backward in turn; ρd=0 corresponds to the d-th bit does not change.


### 4.4. Improvements to the Onlooker Bee Stage

In the traditional ABC algorithm, the onlooker bee usually uses the roulette method to select the food source [14–16]. However, the roulette method has some shortcomings. For example, in the early stage of the iteration, due to the large differences in fitness value, food sources with low fitness will be quickly eliminated, which will destroy the individual diversity; in the later stage of iteration, due to the small difference in fitness value, the selection probability of each food source tends to 1/NP, and the ability to select dominant food sources will be reduced.Aiming at the above problems, we propose a sorting-based adaptive probability (SAP) selection method. The selection probability calculated by this method has no direct relationship with the fitness value but is only related to the order of the dominant food source and the number of iterations [49]. It is an effective method to control the selection ability in the iterative process. Its calculation formula is as follows:
(23)Pct=a+1−a·ebtr−1/TNP−1−1ebt/T−1,where r is the sequence index after sorting all the food sources according to the fitness value; the greater the fitness of the food source, the higher the ranking; a is the lowest selection probability; and b is the adaptive control coefficient. The variables a and b are set as a=0.05 and b=10 [49]. When NP=50, the change curve of the adaptive food selection probability is shown in Figure 5.Figure 5
Sorting-based adaptive probability change curve.It can be seen that the selection probability of food sources calculated by this method is balanced at the beginning of the iteration, which can increase the selection probability of poor food sources, maintain the diversity of the population, and increase the global search ability of the algorithm, and in the later stage of the iteration, which can increase the selection probability of excellent food sources and enable the algorithm to converge quickly.


## 4.1. ABC Algorithm

The ABC algorithm divides bee colony into three categories: employed bee, onlooker bee, and scout bee. The goal of the bee colony is to find the optimal food source, and the food source represents all possible solutions in the solution space and is measured by fitness value. Employed bee focuses on food source detection. Onlooker bee receives food source information shared by other bees and is responsible for mining food sources. Scout bee searches for new food sources randomly when food sources are abandoned. The algorithm process is as follows:(1)
Initialization Stage. In a D-dimensional search space, the population number is NP, and the position of each food source after the t-th iteration is(18)Xit=xi,1t,xi,2t,⋯,xi,Dt,where i=1,2,⋯,NP. The fitness value of food source is fiti=fitXit(2)
Employed Bee Stage. Each employed bee chooses a food source randomly and then generates a new food source vid by the food source update strategy. The fitness fiti′ value of the new food source is calculated, and the food source with better fitness value is retained by the greedy strategy(3)
Onlooker Bee Stage. After receiving the information of the employed bee, the onlooker bee selected several food sources according to the food source selection strategy, searched for new food sources according to the food source update strategy for fitness comparison, and selected the food sources with better fitness for retention(4)
Scout Bee Stage. If the food source did not improve after lim iterations, the food source is abandoned and recorded in the tabu list. At the same time, the bees corresponding to the food source turned into scout bee and generated a new food source according to the initialization strategy randomly(5)
Repeat steps (2) to (4) until the termination condition is met and the optimal food source location is output


## 4.2. Improvements to the Initialization Stage

### 4.2.1. Discrete Integer Coding Strategy Based on Sensor-Subtask Sequence

The decision variableai,jk in the SSA model corresponds to the bee colony individual of the algorithm. For the problem of s sensors, m targets, and n subtasks, if 0-1 coding is used, a D-dimensional (D=s×m×n) vector will be generated, which will cause a dimensional disaster as scene complexity increases. To reduce the computational complexity of the algorithm, a discrete integer coding method based on sensor-subtask sequence is proposed as shown in Figure 3.Figure 3
Discrete integer coding strategy based on sensor-subtask sequence.In this method, the decision matrixAi is encoded as the target number corresponding to the sensor, and the position without a visual relationship is filled with 0. At this time, the coding length is determined by the maximum tracking capability Cmax of each sensor; that is, D=∑k=1sCmaxk. Through this coding method, the dimension of the algorithm is effectively reduced, and the number of targets assigned by each sensor does not exceed the target capacity, which is convenient for directly expressing the scheduling scheme.


### 4.2.2. Heuristic Initialization Rules

The heuristic initialization rules are aimed at generating an initial feasible solution to improve the quality of the initial solution and speed up the convergence [18]. According to the characteristics of SSA model, we propose the initialization rules based on target-threat-priority and resource-balance-priority.(1) Heuristic Initialization Process Based on Target-Threat-Priority.Step 1.
To select the target with the maximum threat degree, assign it to the sensorS1 which has the highest detection priority, and remove the target.Step 2.
Continue to assign the target with the maximum threat degree to the sensorS1 and remove the target until S1 reaches the maximum capacity.Step 3.
According to the methods of Steps1 and 2, assign targets to S2-Sk in turn until there is no target or the sensor resources are saturated.Step 4.
To take the assignment scheme as a heuristic initial individual to replace the initial solution.(2) Heuristic Initialization Process Based on Resource-Balance-Priority.Step 1.
To select the target with the largest detection benefit forS1, assign it to S1, and remove the target; continue to take the target with the largest detection benefit for S2-Sk and assign the target to each sensor, and remove the target. So far, each sensor has been assigned a target with the largest detection benefit.Step 2.
Repeat Step1 until there is no target or the sensor resources are saturated.Step 3.
To take the assignment scheme as a heuristic initial individual to replace the initial solution.


## 4.2.1. Discrete Integer Coding Strategy Based on Sensor-Subtask Sequence

The decision variableai,jk in the SSA model corresponds to the bee colony individual of the algorithm. For the problem of s sensors, m targets, and n subtasks, if 0-1 coding is used, a D-dimensional (D=s×m×n) vector will be generated, which will cause a dimensional disaster as scene complexity increases. To reduce the computational complexity of the algorithm, a discrete integer coding method based on sensor-subtask sequence is proposed as shown in Figure 3.Figure 3
Discrete integer coding strategy based on sensor-subtask sequence.In this method, the decision matrixAi is encoded as the target number corresponding to the sensor, and the position without a visual relationship is filled with 0. At this time, the coding length is determined by the maximum tracking capability Cmax of each sensor; that is, D=∑k=1sCmaxk. Through this coding method, the dimension of the algorithm is effectively reduced, and the number of targets assigned by each sensor does not exceed the target capacity, which is convenient for directly expressing the scheduling scheme.


## 4.2.2. Heuristic Initialization Rules

The heuristic initialization rules are aimed at generating an initial feasible solution to improve the quality of the initial solution and speed up the convergence [18]. According to the characteristics of SSA model, we propose the initialization rules based on target-threat-priority and resource-balance-priority.(1) Heuristic Initialization Process Based on Target-Threat-Priority.Step 1.
To select the target with the maximum threat degree, assign it to the sensorS1 which has the highest detection priority, and remove the target.Step 2.
Continue to assign the target with the maximum threat degree to the sensorS1 and remove the target until S1 reaches the maximum capacity.Step 3.
According to the methods of Steps1 and 2, assign targets to S2-Sk in turn until there is no target or the sensor resources are saturated.Step 4.
To take the assignment scheme as a heuristic initial individual to replace the initial solution.(2) Heuristic Initialization Process Based on Resource-Balance-Priority.Step 1.
To select the target with the largest detection benefit forS1, assign it to S1, and remove the target; continue to take the target with the largest detection benefit for S2-Sk and assign the target to each sensor, and remove the target. So far, each sensor has been assigned a target with the largest detection benefit.Step 2.
Repeat Step1 until there is no target or the sensor resources are saturated.Step 3.
To take the assignment scheme as a heuristic initial individual to replace the initial solution.


## 4.3. Improvements to the Food Source Update Strategy

Since the adoption of the integer coding strategy, the traditional strategy of updating food sources in the employed bee stage and the scout bee stage [47] is not applicable. Based on Zhang et al.’s research [48], a food source update strategy based on discrete differential mutation (DDM) operation is proposed to enhance the algorithm’s search capacity. The food source update formula is as follows:
(19)Vit=φ1Xit⊕φ2⊗Xit−Xkt,where the scale factors ϕ1 and ϕ2 are random numbers in [0, 1], Xit,Xkt,i,k∈1,NP is the i-th food source and the k-th food source, respectively, and Vit is the new food resource. The operation process is divided into the following three parts:Step 1.
Calculate the part ofXit−Xkt.
The operation is defined as(20)Δi=Xit−Xkt⇔δid=xkd,xid≠xkd,randiN,xid=xkd=0,0,xid=xkd≠0,,where d=1,2,⋯,D and Δi=δi1,δi2,⋯,δiD. Through this operation, the codes of the i-th food source and the k-th food source are compared bit by bit. If the code of d-th bit is nonzero and the same, it is 0; if the bit is all zero, it will take a random integer in 1,N; if the bit is different, it will take the bit encoding of the i-th food source.Step 2.
Calculate the part ofϕ2⊗Δi.
The operation is defined as(21)Ρi=φ2⊗Δi⇔ρid=roundδd·DN,φ2>rand,0,otherwise,where Ρi=ρi1,ρi2,⋯,ρiD is the mutation operator and round· is the rounding function. After this operation, the code of a certain bit is probabilistically converted into an integer in 1,D.Step 3.
Calculate the part ofφ1Xit⊕Ρi.
The operation is defined as(22)φ1Xit⊕Ρ=Vi⇔vid=swapxi,ρd,xid,φ1≤rand,insertxi,ρd,xd,φ1>rand,ρd≠0,xid,ρd=0,where swap· is the function to swap two numbers and insert· is the function that inserts the first number before the second number. The process of mutation can be illustrated as follows: if ρd=0, remain xid unchanged; if ρd≠0 and φ1≤rand, then swap the ρd-th and d-th bits of Xit; if ρd≠0 and φ1>rand, insert the ρd-th bit of Xit into the d-th bit, and move the rest of the bits backward in turn, as shown in Figure 4.Figure 4
Food source update operation based on DDM.An illustration is given for the operation of Step3. After the comparison between φ1 and rand, ρ2=8 corresponds to swap xi2 and xi8; ρ4=3 corresponds to insert xi3 into the 4-th bit of Xit, and the other bits are moved backward in turn; ρ7=5 corresponds to swap xi7 and xi5; ρ9=6 corresponds to insert xi6 into the 9-th bit of Xit, and the other bits are moved backward in turn; ρ11=5 corresponds to insert xi11 into the 5-th bit of Xit, and the other bits are moved backward in turn; ρd=0 corresponds to the d-th bit does not change.


## 4.4. Improvements to the Onlooker Bee Stage

In the traditional ABC algorithm, the onlooker bee usually uses the roulette method to select the food source [14–16]. However, the roulette method has some shortcomings. For example, in the early stage of the iteration, due to the large differences in fitness value, food sources with low fitness will be quickly eliminated, which will destroy the individual diversity; in the later stage of iteration, due to the small difference in fitness value, the selection probability of each food source tends to 1/NP, and the ability to select dominant food sources will be reduced.Aiming at the above problems, we propose a sorting-based adaptive probability (SAP) selection method. The selection probability calculated by this method has no direct relationship with the fitness value but is only related to the order of the dominant food source and the number of iterations [49]. It is an effective method to control the selection ability in the iterative process. Its calculation formula is as follows:
(23)Pct=a+1−a·ebtr−1/TNP−1−1ebt/T−1,where r is the sequence index after sorting all the food sources according to the fitness value; the greater the fitness of the food source, the higher the ranking; a is the lowest selection probability; and b is the adaptive control coefficient. The variables a and b are set as a=0.05 and b=10 [49]. When NP=50, the change curve of the adaptive food selection probability is shown in Figure 5.Figure 5
Sorting-based adaptive probability change curve.It can be seen that the selection probability of food sources calculated by this method is balanced at the beginning of the iteration, which can increase the selection probability of poor food sources, maintain the diversity of the population, and increase the global search ability of the algorithm, and in the later stage of the iteration, which can increase the selection probability of excellent food sources and enable the algorithm to converge quickly.


## 5. Simulation Results and Comparisons

In the same experimental environment, we designed a simulation analysis experiment and an algorithm comparison experiment to verify the feasibility of the model and the performance of the algorithm. We used MATLAB 2021 and STK 11.0 to build the simulation scenario and generate target motion data. All experiments were run on a Windows 11 personal computer with a Core i7-11800H, a 2.3 GHz CPU, and 16 GB RAM.


### 5.1. Evaluation Indicators

The performance of the algorithm is evaluated by the following indicators.


#### 5.1.1. Value Rate of SchedulingEVR

EVR is defined as the ratio of the cooperative detection benefit of the subtask assignment in the scheduling scheme to the sum of the cooperative detection benefit of all subtasks in the scheduling period, which is used to reflect whether the result is optimal [50].
(24)EVR=∑i=1NeBTi∑i=1NallBTi,where BTi is the cooperative benefit composed of the detection benefit factor and the threat factor of a single subtask and Ne and Nall are the number of assigned subtasks and the number of all subtasks, respectively.


#### 5.1.2. Warhead Target Assignment RateEWT

EWT is defined as the ratio of the number of the assigned warhead to the total number of targets, which is used to reflect whether the algorithm results are valid, as follows:
(25)EWT=nwNwf,where nw is the number of warhead targets assigned in the scheduling scheme and Nwf is the total number of targets.


### 5.2. Performance Analysis of HDABC Algorithm

This experiment analyzed the performance of the HDABC algorithm when dealing with the SSA problem. Build a medium-scale missile early warning simulation scenario as follows: assume that the attacker launched 5 missiles from 3 launch platforms (the launch time ist0) to strike at D1 and D2 points of the defender. Each missile carried 2 companion decoys; that is, there were at most 45 targets at the same time. The missile parameters are shown in Tables 1 and 2. The defender deployed 2 PBR and 2 XBR for missile early warning operations. Among them, PBR can track up to 10 targets at the same time, and XBR can track up to 5 targets. The basic parameters are shown in Table 3. The simulation scenario is shown in Figure 6.Table 1
Typical trajectory information in the threat airspace.
No.Launch pointStrike pointAltitudeTime of flightOthers1_1(42.3°, -156.7°)(40.2°, -85.7°)1240.6 km1490.7 s1_2 and 1_3 are decoy targets1_4(42.2°, -156.8°)(36.4°, -86.3°)1152.3 km1459.3 s1_5 and 1_6 are decoy targets1_7(42.2°, -156.7°)(36.2°, -86.2°)1209.6 km1491.9 s1_8 and 1_9 are decoy targets1_10(42.3°, -156.8°)(40.1°, -85.6°)1215.3 km1476.4 s1_11 and 1_12 are decoy targets1_13(42.3°, -156.7°)(40.1°, -85.7°)1252.4 km1495.7 s1_14 and 1_15 are decoy targetsTable 2
Missile target threat information (t0+792 s).
No.Friend or foeTarget category informationFlight dataWarhead target confidenceTarget typePredicted rangeSpeed of shutdown pointBallistic inclinationRemaining flight timeDistance between striking point and defense pointFini1_11 (foe)0.901 (ballistic)5783.67 km6.085 km/s36.88°698.7 s12.87 km0.8481_21 (foe)0.551 (ballistic)5772.82 km6.086 km/s36.50°693.8 s24.26 km0.7621_31 (foe)0.401 (ballistic)5781.60 km6.085 km/s36.81°697.6 s17.44 km0.8041_41 (foe)0.901 (ballistic)5967.08 km6.108 km/s37.23°677.3 s5.85 km0.9101_51 (foe)0.551 (ballistic)5969.32 km6.107 km/s37.29°676.9 s10.23 km0.8781_61 (foe)0.401 (ballistic)5958.63 km6.108 km/s37.12°672.4 s14.93 km0.842Description:Fini is calculated according to formula (13), where threat value wi,fi is calculated by normalizing the ratio of the five factors in Table 2 to their global maximum or optimal value, and ui,t=0.2.Table 3
Basic parameters of radars.
NamePositionDetection rangeAzimuth angleElevation angleCapacityPBR126.2°N, 97.4°W4000 km-105~15°0 ~85°10PBR254.8°N, 106.9°W4000 km140~260°0 ~85°10XBR132.5°N, 87.3°W1500 km-98~8°10~85°5XBR245.8°N, 88.6°W1500 km170~276°10~85°5Figure 6
Medium-scale missile early warning simulation scenario.Aftert0+792 s, all targets entered the detection range of all sensors; the subtask sequence was generated according to the strategy described in Section 2, and the scheduling interval t0+792,t0+940 was taken for simulation analysis.The traditional discrete ABC algorithm (ABC), the HDABC algorithm in this paper, the HDABC algorithm without heuristic initialization rules (ABC-1), the HDABC algorithm without DDM (ABC-2), and the HDABC algorithm without SAP (ABC-3) were compared in this section, and the purpose of the comparison was to verify the effectiveness of the improved algorithm. The parameter settings of all algorithms are shown in Table4. The results of different algorithms are compared in Figure 7.Table 4
Parameter settings of all algorithms.
ParametersVariablesValuesThe search dimensionD120The food source scaleNP50The maximum number of iterationsTmax500The maximum number of the food source without updateslim50Figure 7
Comparison of solution results of different ABC algorithms.It can be seen that the result of the HDABC algorithm was the best, and the convergence performance of the other four algorithms was not as good as this algorithm. Among them, the difference between HDABC and ABC-1 was small, because the heuristic initialization rule made the HDABC algorithm has a better food source in the early iteration, and the convergence speed was improved; the big difference between HDABC and ABC-2 was that the addition of the DDM food source update strategy increased the algorithm optimization ability.Monte Carlo simulation was performed 20 times for the above simulation scenario, and the average value was taken. The simulation comparison results of the five algorithms are shown in Table5. The calculation results showed that all these ABC algorithms could obtain effective scheduling schemes, and the HDABC algorithm had obvious improvements in various indicators. However, as the complexity of the algorithm increased, the running time of the HDABC algorithm increased, but it was within the acceptable range.Table 5
Comparison table of evaluation results of ABC algorithms.
AlgorithmBest fitness valueAverage fitness valueEVREWTAverage running timeHDABC24.718324.232995.69%92.33%14.67 sABC-124.718323.875393.03%89.33%14.24 sABC-224.703823.159291.97%87.33%10.86 sABC-324.635923.053790.97%86.67%14.31 sABC24.235822.598484.35%83.67%9.83 s


### 5.3. Comparative Analysis of Different Algorithms

This experiment compared the performance of different algorithms and the HDABC algorithm when dealing with SSA problem of different scales. We used the HDABC algorithm, the IGA algorithm [20], the SAPPSO algorithm [21], and the PSO-VNS algorithm [22] for comparison. The purpose of the comparison was to verify the ability of the proposed HDABC algorithm in solving SSA problem.The scenarios were set to small-scale (1 PBR, 1 XBR, 5 targets, and 10 decoys), medium-scale (2 PBR, 2 XBR, 15 targets, and 30 decoys), and large-scale operation (4 PBR, 4 XBR, 30 targets, 60 decoys) scenarios. The small-scale and large-scale simulation scenarios are shown in Figure8, and the medium-scale scenario is the same as Figure 6. The parameter settings of all algorithms for the medium-scale scenario are shown in Table 6.Figure 8
Small-scale and large-scale missile early warning simulation scenarios.Table 6
Parameter settings of all algorithms for the medium-scale scenario.
AlgorithmsParameter settingHDABCParameter settings are the same as in Table4.IGAThe populationNC is 50, the largest genetic iterations GEN is 500, the length of chromosome m is 120, the crossover probability Pc1 is 0.9, Pc2 is 0.6, the mutation probability Pm1 is 0.1, and Pm2 is 0.01.SAPPSOThe number of particles is 50, the total number of iterations is 500, the particle dimension is 180, and the weight factor of particle update is 0.5.PSO-VNSThe population sizeP is 120, the maximum number of iterations G is 200, the inertia weight w is 0.9, the learning factor c1 is 2, c2 is 2, the period of evolutionary stagnation gen is 10, the number of elite solutions N is 25, and the maximum number of iterations in VNS MaxIter is 5.The results of different algorithms are compared in Figure9. It can be seen from the results that HDABC had a better convergence ability compared with other algorithms in dealing with problems of different scales. For medium-scale and large-scale SSA problems, the HDABC algorithm had the fast convergence speed and can get the best solution.Figure 9
Comparison of solution results of different algorithms.
(a)
Small scale(b)
Medium scale(c)
Large scaleAfter performing 20 Monte Carlo simulations on the above small-, medium-, and large-scale SSA problems and taking the average value, the box plots for comparative analysis are shown in Figure10, and the analysis results are shown in Table 7. For the small-scale SSA problem, the four algorithms can obtain the optimal solution, but the IGA algorithm was relatively poor, and the convergence result was unstable. For medium- and large-scale SSA problems, the HDABC algorithm can find the optimal solution compared with other algorithms and had better convergence, and the algorithm results were more reliable. From the experimental data of large-scale scenarios, it can be concluded that the value rate of schedule obtained by the HDABC algorithm was 16.72%, 10.74%, and 11.41% higher than IGA, PSO-VNS, and SAPPSO algorithms, and the warhead target assignment rate was 16.44%, 11.70%, and 13.55% higher, respectively. In conclusion, the advantages of the HDABC algorithm were more reflected in the solution of medium-scale and large-scale SSA problems. However, the efficiency of the HDABC algorithm needs to be further improved.Figure 10
Box plot of algorithm comparison of different scales.
(a)
Small scale(b)
Medium scale(c)
Large scaleTable 7
Comparison table of evaluation results of different algorithms.
ScaleAlgorithmBest fitness valueAverage fitness valueEVREWTAverage running timeSmall scaleHDABC13.831313.8313100.00%100.00%6.05 sIGA13.831313.294097.63%98.67%11.35 sPSO-VNS13.831317.810298.88%99.33%5.37 sSAPPSO13.831313.8313100.00%100.00%4.34 sMedium scaleHDABC24.718524.509296.22%92.67%14.32 sIGA24.647523.896493.13%89.67%37.92 sPSO-VNS24.865724.333294.86%90.33%13.23 sSAPPSO24.381824.047493.09%89.33%9.28 sLarge scaleHDABC51.260151.012993.51%91.37%31.26 sIGA49.751348.628376.79%74.93%63.35 sPSO-VNS51.209149.784382.77%79.67%29.35 sSAPPSO50.106449.579982.10%77.82%18.87 s


## 5.1. Evaluation Indicators

The performance of the algorithm is evaluated by the following indicators.


### 5.1.1. Value Rate of SchedulingEVR

EVR is defined as the ratio of the cooperative detection benefit of the subtask assignment in the scheduling scheme to the sum of the cooperative detection benefit of all subtasks in the scheduling period, which is used to reflect whether the result is optimal [50].
(24)EVR=∑i=1NeBTi∑i=1NallBTi,where BTi is the cooperative benefit composed of the detection benefit factor and the threat factor of a single subtask and Ne and Nall are the number of assigned subtasks and the number of all subtasks, respectively.


### 5.1.2. Warhead Target Assignment RateEWT

EWT is defined as the ratio of the number of the assigned warhead to the total number of targets, which is used to reflect whether the algorithm results are valid, as follows:
(25)EWT=nwNwf,where nw is the number of warhead targets assigned in the scheduling scheme and Nwf is the total number of targets.


## 5.1.1. Value Rate of SchedulingEVR

EVR is defined as the ratio of the cooperative detection benefit of the subtask assignment in the scheduling scheme to the sum of the cooperative detection benefit of all subtasks in the scheduling period, which is used to reflect whether the result is optimal [50].
(24)EVR=∑i=1NeBTi∑i=1NallBTi,where BTi is the cooperative benefit composed of the detection benefit factor and the threat factor of a single subtask and Ne and Nall are the number of assigned subtasks and the number of all subtasks, respectively.


## 5.1.2. Warhead Target Assignment RateEWT

EWT is defined as the ratio of the number of the assigned warhead to the total number of targets, which is used to reflect whether the algorithm results are valid, as follows:
(25)EWT=nwNwf,where nw is the number of warhead targets assigned in the scheduling scheme and Nwf is the total number of targets.


## 5.2. Performance Analysis of HDABC Algorithm

This experiment analyzed the performance of the HDABC algorithm when dealing with the SSA problem. Build a medium-scale missile early warning simulation scenario as follows: assume that the attacker launched 5 missiles from 3 launch platforms (the launch time ist0) to strike at D1 and D2 points of the defender. Each missile carried 2 companion decoys; that is, there were at most 45 targets at the same time. The missile parameters are shown in Tables 1 and 2. The defender deployed 2 PBR and 2 XBR for missile early warning operations. Among them, PBR can track up to 10 targets at the same time, and XBR can track up to 5 targets. The basic parameters are shown in Table 3. The simulation scenario is shown in Figure 6.Table 1
Typical trajectory information in the threat airspace.
No.Launch pointStrike pointAltitudeTime of flightOthers1_1(42.3°, -156.7°)(40.2°, -85.7°)1240.6 km1490.7 s1_2 and 1_3 are decoy targets1_4(42.2°, -156.8°)(36.4°, -86.3°)1152.3 km1459.3 s1_5 and 1_6 are decoy targets1_7(42.2°, -156.7°)(36.2°, -86.2°)1209.6 km1491.9 s1_8 and 1_9 are decoy targets1_10(42.3°, -156.8°)(40.1°, -85.6°)1215.3 km1476.4 s1_11 and 1_12 are decoy targets1_13(42.3°, -156.7°)(40.1°, -85.7°)1252.4 km1495.7 s1_14 and 1_15 are decoy targetsTable 2
Missile target threat information (t0+792 s).
No.Friend or foeTarget category informationFlight dataWarhead target confidenceTarget typePredicted rangeSpeed of shutdown pointBallistic inclinationRemaining flight timeDistance between striking point and defense pointFini1_11 (foe)0.901 (ballistic)5783.67 km6.085 km/s36.88°698.7 s12.87 km0.8481_21 (foe)0.551 (ballistic)5772.82 km6.086 km/s36.50°693.8 s24.26 km0.7621_31 (foe)0.401 (ballistic)5781.60 km6.085 km/s36.81°697.6 s17.44 km0.8041_41 (foe)0.901 (ballistic)5967.08 km6.108 km/s37.23°677.3 s5.85 km0.9101_51 (foe)0.551 (ballistic)5969.32 km6.107 km/s37.29°676.9 s10.23 km0.8781_61 (foe)0.401 (ballistic)5958.63 km6.108 km/s37.12°672.4 s14.93 km0.842Description:Fini is calculated according to formula (13), where threat value wi,fi is calculated by normalizing the ratio of the five factors in Table 2 to their global maximum or optimal value, and ui,t=0.2.Table 3
Basic parameters of radars.
NamePositionDetection rangeAzimuth angleElevation angleCapacityPBR126.2°N, 97.4°W4000 km-105~15°0 ~85°10PBR254.8°N, 106.9°W4000 km140~260°0 ~85°10XBR132.5°N, 87.3°W1500 km-98~8°10~85°5XBR245.8°N, 88.6°W1500 km170~276°10~85°5Figure 6
Medium-scale missile early warning simulation scenario.Aftert0+792 s, all targets entered the detection range of all sensors; the subtask sequence was generated according to the strategy described in Section 2, and the scheduling interval t0+792,t0+940 was taken for simulation analysis.The traditional discrete ABC algorithm (ABC), the HDABC algorithm in this paper, the HDABC algorithm without heuristic initialization rules (ABC-1), the HDABC algorithm without DDM (ABC-2), and the HDABC algorithm without SAP (ABC-3) were compared in this section, and the purpose of the comparison was to verify the effectiveness of the improved algorithm. The parameter settings of all algorithms are shown in Table4. The results of different algorithms are compared in Figure 7.Table 4
Parameter settings of all algorithms.
ParametersVariablesValuesThe search dimensionD120The food source scaleNP50The maximum number of iterationsTmax500The maximum number of the food source without updateslim50Figure 7
Comparison of solution results of different ABC algorithms.It can be seen that the result of the HDABC algorithm was the best, and the convergence performance of the other four algorithms was not as good as this algorithm. Among them, the difference between HDABC and ABC-1 was small, because the heuristic initialization rule made the HDABC algorithm has a better food source in the early iteration, and the convergence speed was improved; the big difference between HDABC and ABC-2 was that the addition of the DDM food source update strategy increased the algorithm optimization ability.Monte Carlo simulation was performed 20 times for the above simulation scenario, and the average value was taken. The simulation comparison results of the five algorithms are shown in Table5. The calculation results showed that all these ABC algorithms could obtain effective scheduling schemes, and the HDABC algorithm had obvious improvements in various indicators. However, as the complexity of the algorithm increased, the running time of the HDABC algorithm increased, but it was within the acceptable range.Table 5
Comparison table of evaluation results of ABC algorithms.
AlgorithmBest fitness valueAverage fitness valueEVREWTAverage running timeHDABC24.718324.232995.69%92.33%14.67 sABC-124.718323.875393.03%89.33%14.24 sABC-224.703823.159291.97%87.33%10.86 sABC-324.635923.053790.97%86.67%14.31 sABC24.235822.598484.35%83.67%9.83 s


## 5.3. Comparative Analysis of Different Algorithms

This experiment compared the performance of different algorithms and the HDABC algorithm when dealing with SSA problem of different scales. We used the HDABC algorithm, the IGA algorithm [20], the SAPPSO algorithm [21], and the PSO-VNS algorithm [22] for comparison. The purpose of the comparison was to verify the ability of the proposed HDABC algorithm in solving SSA problem.The scenarios were set to small-scale (1 PBR, 1 XBR, 5 targets, and 10 decoys), medium-scale (2 PBR, 2 XBR, 15 targets, and 30 decoys), and large-scale operation (4 PBR, 4 XBR, 30 targets, 60 decoys) scenarios. The small-scale and large-scale simulation scenarios are shown in Figure8, and the medium-scale scenario is the same as Figure 6. The parameter settings of all algorithms for the medium-scale scenario are shown in Table 6.Figure 8
Small-scale and large-scale missile early warning simulation scenarios.Table 6
Parameter settings of all algorithms for the medium-scale scenario.
AlgorithmsParameter settingHDABCParameter settings are the same as in Table4.IGAThe populationNC is 50, the largest genetic iterations GEN is 500, the length of chromosome m is 120, the crossover probability Pc1 is 0.9, Pc2 is 0.6, the mutation probability Pm1 is 0.1, and Pm2 is 0.01.SAPPSOThe number of particles is 50, the total number of iterations is 500, the particle dimension is 180, and the weight factor of particle update is 0.5.PSO-VNSThe population sizeP is 120, the maximum number of iterations G is 200, the inertia weight w is 0.9, the learning factor c1 is 2, c2 is 2, the period of evolutionary stagnation gen is 10, the number of elite solutions N is 25, and the maximum number of iterations in VNS MaxIter is 5.The results of different algorithms are compared in Figure9. It can be seen from the results that HDABC had a better convergence ability compared with other algorithms in dealing with problems of different scales. For medium-scale and large-scale SSA problems, the HDABC algorithm had the fast convergence speed and can get the best solution.Figure 9
Comparison of solution results of different algorithms.
(a)
Small scale(b)
Medium scale(c)
Large scaleAfter performing 20 Monte Carlo simulations on the above small-, medium-, and large-scale SSA problems and taking the average value, the box plots for comparative analysis are shown in Figure10, and the analysis results are shown in Table 7. For the small-scale SSA problem, the four algorithms can obtain the optimal solution, but the IGA algorithm was relatively poor, and the convergence result was unstable. For medium- and large-scale SSA problems, the HDABC algorithm can find the optimal solution compared with other algorithms and had better convergence, and the algorithm results were more reliable. From the experimental data of large-scale scenarios, it can be concluded that the value rate of schedule obtained by the HDABC algorithm was 16.72%, 10.74%, and 11.41% higher than IGA, PSO-VNS, and SAPPSO algorithms, and the warhead target assignment rate was 16.44%, 11.70%, and 13.55% higher, respectively. In conclusion, the advantages of the HDABC algorithm were more reflected in the solution of medium-scale and large-scale SSA problems. However, the efficiency of the HDABC algorithm needs to be further improved.Figure 10
Box plot of algorithm comparison of different scales.
(a)
Small scale(b)
Medium scale(c)
Large scaleTable 7
Comparison table of evaluation results of different algorithms.
ScaleAlgorithmBest fitness valueAverage fitness valueEVREWTAverage running timeSmall scaleHDABC13.831313.8313100.00%100.00%6.05 sIGA13.831313.294097.63%98.67%11.35 sPSO-VNS13.831317.810298.88%99.33%5.37 sSAPPSO13.831313.8313100.00%100.00%4.34 sMedium scaleHDABC24.718524.509296.22%92.67%14.32 sIGA24.647523.896493.13%89.67%37.92 sPSO-VNS24.865724.333294.86%90.33%13.23 sSAPPSO24.381824.047493.09%89.33%9.28 sLarge scaleHDABC51.260151.012993.51%91.37%31.26 sIGA49.751348.628376.79%74.93%63.35 sPSO-VNS51.209149.784382.77%79.67%29.35 sSAPPSO50.106449.579982.10%77.82%18.87 s


## 6. Conclusions

In this paper, some exploratory research is carried out on the problem of missile early warning multisensor resource scheduling. The multisensor and multitarget resource scheduling problem is transformed into a sensor-subtask assignment problem through periodic scheduling-task decomposition, and the SSA model of this problem is established to solve the cooperative scheduling problem under the operational background of the missile early warning system. Through the adaptive improvement of the ABC algorithm, it has better performance in dealing with such problems. Compared with other algorithms, the HDABC algorithm has the advantages of fast convergence speed, high solution accuracy, and good search performance; in addition, the HDABC algorithm has great advantages in solving the problems of large-scale missile early warning operation.However, this paper has the following shortcomings: the analysis is mainly carried out on the cooperative scheduling problem of ground-based early warning radars, such as PBR and XBR; the established model and simulation scenarios are simple; the efficiency of the HDABC algorithm needs to be further improved. In the follow-up work, we will further study the cooperative resource scheduling method of multisource sensors such as early warning satellites and early warning radars in complex scenarios.



---

*Source: 5094415-2022-10-25.xml*