A multi-objective scheduling method for operational coordination time using improved triangular fuzzy number representation.

In modern warfare, the comprehensiveness of combat domain and the complexity of tasks pose great challenges to operational coordination.To address this challenge, we use the improved triangular fuzzy number to express the combat mission time, first present a new multi-objective operational cooperative time scheduling model that takes the fluctuation of combat coordinative time and the time flexibility between each task into account. The resulting model is essentially a large-scale multi-objective combinatorial optimization problem, intractably complicated to solve optimally. We next propose multi-objective improved Bat algorithm based on angle decomposition (MOIBA/AD) to quickly identify high-quality solutions to the model. Our proposed algorithm improves the decomposition strategy by replacing the planar space with the angle space, which helps greatly reduce the difficulty of processing evolutionary individuals and hence the time complexity of the multi-objective evolutionary algorithm based on decomposition (MOEA/D). Moreover, the population replacement strategy is enhanced utilizing the improved bat algorithm, which helps evolutionary individuals avoid getting trapped in local optima. Computational experiments on multi-objective operational cooperative time scheduling (MOOCTS) problems of different scales demonstrate the superiority of our proposed method over four state-of-the-art multi-objective evolutionary algorithms (MOEAs), including multi-objective bat Algorithm (MOBA), MOEA/D, non-dominated sorting genetic algorithm version II (NSGA-II) and multi-objective particle swarm optimization algorithm (MOPSO). Our proposed method performs better in terms of four performance criteria, producing solutions of higher quality while keeping a better distribution of the Pareto solution set.

1 Introduction

Operational coordination refers to the coordination and cooperation among various combat forces in actions according to a unified plan when they jointly perform combat tasks, so as to ensure the coordinated actions of various combat forces and exert their overall combat effectiveness [1]. With the development of miniaturization and intelligence of the war, it is of great significance to study the combat coordination in tactical operations for improving the combat capability of the troops.

In the field of multi-agent model optimization and network task optimization, the problem of combat task scheduling has been widely studied. Multi-agent-based task scheduling [2, 3] makes rational planning for task subject through optimization algorithm, and obtains the reasonable task scheduling result. Task scheduling based on time network [4, 5] connects different task time nodes into a directed network, and obtains a reasonable task timing scheduling result by optimizing the parameters of the nodes and the whole network. Although there are many solutions to the time series planning of combat coordination, under the above-mentioned operational background, we can draw the conclusion that solving the time problem of combat coordination in Observe-Orient-Decide-Act (OODA) loop mainly has the following challenges.

When solving the operational coordination time sequence, how to fully consider the complex factors affecting operations, and reasonably model the execution time and completion time series of each operational task based on the complex factors, so as to make it closer to the actual combat;

When considering the execution time of all combat tasks, how to better consider the execution time of each combat task, and whether the solved task time and combat sequence are still valid and the robustness of the model is good enough under the sudden change of battlefield situation;

When solving the established model, how to choose a suitable algorithm to solve the model quickly and effectively.

The purpose of this paper aims to solve the problem of cooperation time and cooperation time sequence in complex combat. The major innovation and contributions of this paper can be summarized as follows. Firstly, considering the complexity of battle, a MOOCTS model is established. In this model, the improved triangular fuzzy number is used to re-describe the combat mission time, and the ranking function of the triangular fuzzy number and the calculation method of the probability degree are given, and the flexible combat time among various combat missions is considered. Secondly, the approach based on angle space decomposition improves the aggregation method of the traditional MOEA/D, and introduces the solution space optimization strategy of the traditional bat algorithm (BA), proposes an improved MOBA. Finally, the improved MOEA is used to solve the MOOCTS model, and several common MOEAs are used to solve this model. The performance and robustness of the algorithms are compared and discussed.

The outline of the paper is as follows: Section 2 will briefly review previous work. In Section 3, a MOOCTS model is established. Section 4 presents the proposed solution algorithm, MOIBA/AD. Some numerical experiment is carried in Section 5. The conclusion and future work are presented in Section 6.

2 Related work

2.1 Operational task time sequence cooperative scheduling

At present, when studying the scheduling of combat coordination, there are mainly two ways of thinking. One is based on the task chain formed by different tasks, which takes the completion cost, completion efficiency and resource utilization efficiency of combat tasks as the optimization objectives, and carries out collaborative planning through the time sequence relationship and resource constraints in the task chain, thus obtaining the overall time scheduling strategy. For example, Zhou et al. [6] took the operational time and operational resource utilization as the cooperative optimization objectives, the role-based artificial bee colony (RABC) algorithm is used to calculate the cooperative time plan. Wang et al. [7] took the operational timeliness as their optimization objective, and based on considering the constraints such as task timing, task completion effect and task execution conflict, they used the invasive weed bat and twin group optimization (IWBDA) method to plan and solve the Operational task time sequence cooperative scheduling problem. Cai et al. [8] used multi-dimensional asynchrony and inertia weight to improve the particle swarm optimization algorithm (PSO), the problem of timing relationship planning in the command and control tasks is solved. Li et al. [9] considered the priorities of different combat tasks based on the traditional collaborative planning problem, and planned the task sequence according to different priorities, which was more in line with the actual combat. In addition, based on considering the priority of tasks and resources at the same time, many scholars use greedy strategy to locally search for the solution space of the multi-imensional dynamic list scheduling model (MDLS) [10-12] and multi-PRI list dynamic scheduling model (MPLDS) [13-15], which have achieved good results in the scheduling problem of cooperative task time. However, this method does not take into account the coupling between the execution subjects and time of different and the same combat tasks in the solving process. It basically believes that similar combat times will not cause conflicts between the same task subjects. But it is obviously not applicable to simply consider the time conflicts among various task chains in actual combat situations.

Another is to combine the combat processes with the task nodes and task arcs to form a simple time network (STN) [16] or a temporal constraint network (TCN) [17], which converts the task timing relationship solving into a network time conflict detection problem. When solving this kind of problems, Tang et al. [18] modeled the combat task time series using STN, and then resolved the conflicts in the time series network using constraint relaxation. Hunsberger et al. [19] improved the traditional STN network using dynamic consistency (DC) and obtained the conditional simple temporary networks (CSTNS), which further improved the time detection performance of the traditional network. Based on STN model, Ding et al. [20] adjusted the task priority and constraint in the model to solve the problem of time conflict. Jia et al. [21] and Huntzberg et al. [22] respectively modeled the time constraints of multi-level joint tasks and multi-unmanned submarine cooperative combat tasks using TCN, and successfully detected and resolved the time conflicts of the tasks. Chen et al. [23] simplified TCN to a certain extent, and let the attack time combined with the weight matrix method by using the simple time constraint network (STCN) according to the combat mode of the air force on the shore, which successfully detected the time consistency and solves the conflict of combat time. Qi et al. [24] further expanded TCN, and considering the complexity of large-scale task time planning, proposed a Hierarchical Task Network (HTN) based on hierarchical partition, which solved the problem well. In the second way of solving this kind of problem, the relationship between actual combat mission is too simple. Many arm of the services have mutual cooperation and support relations, the simple task chain of a single arm basically does not exist.

2.2 MOEAs for operational task scheduling

In the field of solving the problem of combat task scheduling, when the first idea in Section 2.1 is used to model, and the model objective is two or more, the multi-objective optimization method is needed to solve the model. Wu et al. [25] modeled the task scheduling with the preference of decision makers, and used improved non-dominated sorting genetic algorithm version III (NSGA-III) with knee points to effectively solve the problem. Ramezani et al. [26] established a multi-task sequential scheduling model in a distributed computing environment, and two multi-objective optimization algorithms, multi-objective genetic algorithm (MOGA) and MOPSO, are used to compare the model solution and algorithm performance. The results show that the algorithm performance of MOPSO is superior to that of MOGA in solving such problems. Cao et al. [27] and Liu et al. [28] made reasonable assumptions for multi-task scheduling in a dynamic environment and established a sequential scheduling model. By using MOEA/D method, the sequential scheduling scheme of tasks was obtained. In addition, NSGA-II [29, 30], strength Pareto evolutionary algorithm version II (SPEA-II) [31] and multi-objective artificial bee colony (MOABC) [32] are also used to solve multi-task scheduling. At present, the application of various multi-objective optimization algorithms in this kind of problems is still a very hot field, which is worth further study.

In addition, it is one of the widely used methods to evaluate and rank combat tasks by using various evaluation methods. For example, an evaluation method combining the Analytical Hierarchy Process (AHP) method and fuzzy mathematics [33-35], Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) method in dynamic environment [36-38], Data Envelopment Analysis (DEA) method which is less affected by the subjective factors of decision makers [39-41], and entropy method combined with the above method [42, 43], etc.

3 MOCCTS

Generally, a combat mission requires a complete OODA loop, as shown in Fig 1. Taking a joint land-air combat operation as an example, the combat weapons under the jurisdiction of our army can be divided into three clusters, namely, one land-air combat cluster, one air combat cluster and one command and control support cluster. Each cluster is internally connected through a communication data link (dotted blue line in the figure), and the clusters are connected together through sensors and data links (dotted red line in the figure). When carrying out tasks against operational targets, the three clusters conduct unified operations, such as reconnaissance, surveillance, and strikes. We continuously observe the combat targets through various investigation means to obtain detailed information of the combat targets, and feed back the information to a decision module consisting of a command decision maker and an auxiliary decision system, then the decision module comprehensively analyzes and judges the battlefield situation and the combat target information, obtains the importance level of each combat task faced by us, and further makes decisions on combat operations. The contents of combat decision-making generally includes the assignment of weapon target at each combat time, assignment of weapon target changing with dynamic time, duration of each combat task, coordinative time sequence of combat task and so on. This paper mainly studies the solution of the latter two problems, which are represented by yellow arrows in the figure.

Fig 1

An illustration of operational coordination content in OODA loop.

3.1 Model assumptions and notation descriptions

The notation description of this model is shown in Table 1.

Table 1

Notation declaration.

Notation	Description
I	A set of sequence numbers for a combat mission, I = {1, 2, ⋯, n}, i, i′, i″ ∈ I;
J	A set of sequence numbers for the combat phase, J = {1, 2, ⋯, l}, j, j′, j″ ∈ J;
K	The sequence number set of the combat cluster, K = {1, 2, ⋯, m}, k, k′, k″ ∈ K;
Sij, Si′j′, Si″j″	A set of combat tasks that need to be completed in different combat phases;
PTijk	Duration for combat cluster k to complete combat mission i in phase j;
CTijk	Combat cluster k completes the end time of combat mission i in stage j;
Ni	Number of combat clusters required at least to complete combat mission i;
Ck	Combat consumption per unit time generated by combat cluster k;
Iij	Mission i importance factor in combat phase j, Iij ∈ [0, 1];
Rj^′,j^″i^′,i^″	Sequencing variable of combat mission, Rj^′,j^″i^′,i^″=1 indicates that task i′ in combat phase j′ occurs before task i′ in combat phase j′; otherwise, Rj^′,j^″i^′,i^″=0;
qijk	Correspondence variable between combat cluster and combat task, qijk=1 indicates that task i in operation phase j is completed by operation cluster k; otherwise, qijk=0.

In order to make the model established in this paper closer to the actual combat and convenient for solving the model, the following explanations and assumptions are put forward.

In this model, combat clusters are regarded as the basic unit of operational coordination. Generally speaking, the task coordination scheme is to plan the task order for different specialties and combat clusters. When executing a specific task, the selection of the internal strength of the cluster is selected by the cluster commander according to the actual situation.

Each combat cluster adopts the “Start-End” mechanism when executing combat tasks, that is, once a combat task starts to execute, the clusters that complete the task must execute the current task to the end, and the execution of the current task can not be interrupted.

The importance ranking of all combat missions is known prior to model solving.

Each combat cluster can only carry out one combat task simultaneously in the same combat phase.

The occurrence of our equipment failure during combat is not considered.

3.2 Model formulation

3.2.1 Objectives

This paper mainly considers the following two optimization objectives, that is, the maximum completion time of all tasks and the cost-effectiveness ratio of the whole combat operation.

The maximum completion time of all tasks refers to the end time of the last task in the whole operational process. We express it by the longest duration in all combat phases when k combat cluster performs i combat mission. Expressed as follows.

The cost-effectiveness ratio of the whole combat operation refers to the ratio of combat effectiveness to combat consumption in the course of combat. Because the whole combat mission coordination time network is a time series directed graph composed of several subtasks, in the classical network analysis method, the “entropy” is generally used to measure the influence and coordination degree of the network system on the subsystems [44]. Cost-effectiveness ratio is one of the performance indexes of the whole cooperative time series network, and we also use the entropy value to express the traditional cost-effectiveness ratio. So, the combat effectiveness is expressed by the important entropy of the combat mission. where indicates the importance of combat cluster k completing task i in phase j. Obviously, the smaller the value, the higher the combat effectiveness, and the larger the value, the lower the combat effectiveness. Therefore, in practice, we take the opposite number. Then, Considering the combat consumption C together with the combat importance entropy defined above, the cost-effectiveness ratio of the whole operation is expressed as follows. The denominator in the above formula represents the sum of our combat consumption in the whole combat stage.

3.2.2 Constraints

Constraint (4) means that each combat mission cannot be executed repeatedly during the battle process. Constraint (5) means that each battle cluster can only perform one combat task in each stage. Constraint (6) means that each combat cluster can not perform a combat task in a different combat stages. Constraint (7) indicates that the number of clusters that complete combat tasks can not be lower than the threshold of task completion. Constraint 8 is the decision variable constraint of the model, so the value of the variable can be between 0 and 1.

3.2.3 Model

The model is given as follows: For the MOCCTS model, the first objective function is multiplied by -1, i. e., from minimum optimization to maximum optimization. In the case of two maximized objective functions, in the objective space composed of the same decision variables, when one objective reaches the optimal value (maximum value), the other objective function can not reach the optimal value, which leads to the paradox that both objective functions taking the maximum value at the same time. Obviously, it is not appropriate to aggregate multiple objective functions into one objective by using the weight sum approach. From the search direction of the optimal value of the function, two objective functions form a two-dimensional objective space. In this space, a compromise solution is needed to find the solution set to deal with two objective functions simultaneously. The solution set is called Pareto front (PF) [45].

3.3 Definition and solution of cooperative time

In the combat process, there is the possibility of sudden changes in the environment, opponents and our fighting strength, which leads to certain fluctuations in the actual execution time of combat missions. If the execution time of task is fixed, it will be inconsistent with the actual situation, resulting in the solution result of the model not being operable. So, in this paper, the improved triangular fuzzy numbers (ITFN) is used to represent the execution time of combat mission, that is where respectively represent the optimistic time, the most likely time, and the pessimistic time for performing combat tasks. If the membership function is used to represent , that is

In the previous studies on TFN, Sakawa et al. [46] proposed a measurement method and a linear operation method, and Seikh et al. [47] extended TFN in an intuitionistic fuzzy environment, and obtained the definition and calculation method of triangular intuitivistic fuzzy number (TIFN). On this basis, the operation mode of TFN is constantly being explored [48, 49]. In this paper, the structure of triangular fuzzy numbers is appropriately improved according to the needs of the model firstly, then a sort function of ITFN and a calculation method of sort probability are proposed.

definition 1 TFN is improved by using flexible time to define ITFN. Assuming that ITFN of task i′ and task i″ are and respectively, we define the maximum and minimum combat elastic times between the two tasks as respectively, which means that if the maximum duration after the start of combat task i′ is , then combat task i″ must start; task i″ can only be started when the duration of task i is at least . As shown in Fig 2.

Fig 2

An illustration of using flexible time to improve TFN diagram.

, if there is a minimum duration between combat task i′ and task i″; If there is a longest duration between combat task i′ and task i″, . Especially, if the combat task i′ and task i″ have multiple minimum durations , then ; if there are multiple maximum durations between combat task i′ and task i″, .

definition 2 The calculation rules of ITFN [50]. For the convenience of the following calculations, define four operations: equation, addition (subtraction), multiplication and division.

(1) Equation. if and only if condition is met.

(2) Addition (Subtraction). .

(3) Multiplication. , especially, .

(4) Division. .

definition 3 Sorting function for ITFN. For ∀i ∈ I, j ∈ J, k ∈ K, , the function calls the ITFN sorting function of .

theorem 1 Sorting probability degree of ITFN. For ∀i′, i″ ∈ I, j′, j″ ∈ J, k′, k″ ∈ K, S is the projection on the xoy surface of the three-dimensional Cartesion coordinates when the sorting function . S′ is the projection on the xoy surface of the three-dimensional Cartesion coordinates when the sorting function , the probability degree , .

Proof 3.1 (Proof of Theorem 1) As shown in Fig 3, sorting functions and are plotted in a three-dimensional cartesian coordinate system to obtain ΔABC and ΔDEF. The intersecting line of the two is GH, i.e. . As can be seen from Fig 3, SΔ = S + S = S + S′, so, , .

Fig 3

Schematic diagram for proving theorem 1.

theorem 2 Solution formula of probability degree. For ∀i′, i″ ∈ I, j′, j″ ∈ J, k′, k″ ∈ K, sorting functions and , x, y ∈ [0, 1], let , then, the formulae for solving the ITFN probability degree under different conditions is shown in Table 2.

Table 2

Possibility solving formulae for different ITFN.

The size relationship of PTi^′j^′k^′ and PTi^″j^″k^″	Solution result of probability
{pti^′j^′k^′(1)=pti^″j^″k^″(1)pti^′j^′k^′(2)≥pti^″j^″k^″(2)pti^′j^′k^′(3)≥pti^″j^″k^″(3)	P(PTi^′j^′k^′occursbeforePTi^″j^″k^″)=1
{pti^′j^′k^′(1)≤pti^″j^″k^″(1)pti^′j^′k^′(2)≤pti^″j^″k^″(2)pti^′j^′k^′(3)≤pti^″j^″k^″(3)	P(PTi^″j^″k^″occursbeforePTi^′j^′k^′)=1
{pti^′j^′k^′(1)≤pti^″j^″k^″(1)pti^′j^′k^′(2)≤pti^″j^″k^″(2)pti^′j^′k^′(3)>pti^″j^″k^″(3)	{P(PTi^′j^′k^′occursbeforePTi^″j^″k^″)=(C-F)^2(C+E-B-F)(C+D-A-F)P(PTi^″j^″k^″occursbeforePTi^′j^′k^′)=(C+E-B-F)(C+D-A-F)-(C-F)^2(C+E-B-F)(C+D-A-F)
{pti^′j^′k^′(1)<pti^″j^″k^″(1)pti^′j^′k^′(2)>pti^″j^″k^″(2)pti^′j^′k^′(3)≤pti^″j^″k^″(3)	{P(PTi^′j^′k^′occursbeforePTi^″j^″k^″)=(B-E)^2(D+B-E-A)(F+B-E-C)P(PTi^″j^″k^″occursbeforePTi^′j^′k^′)=(D+B-E-A)(F+B-E-C)-(B-E)^2(D+B-E-A)(F+B-E-C)
{pti^′j^′k^′(1)<pti^″j^″k^″(1)pti^′j^′k^′(2)>pti^″j^″k^″(2)pti^′j^′k^′(3)>pti^″j^″k^″(3)	{P(PTi^′j^′k^′occursbeforePTi^″j^″k^″)=(D-A)^2(C+D-A-F)(B+D-A-E)P(PTi^″j^″k^″occursbeforePTi^′j^′k^′)=(C+D-A-F)(B+D-A-E)-(D-A)^2(C+D-A-F)(B+D-A-E)

Proof 3.2 (Proof of Theorem 2) The following proof illustrates only the third case in Table 2. As shown in Fig 3, the end point of the projection JI of the intersecting line GH is on the line MN, and the coordinate of the end point is (x0, y0), then the coordinate should be the solution of the following equation. The endpoint coordinate is . Similarly, it can be solved when the endpoint of the projection JI is on the line segments OM and ON, the endpoint coordinates are , . The formula for solving the possibility degree can be obtained by substitution.

4 A multi-objective improved bat algorithm based on angle decomposition

As shown above, the MOOCTS model is in fact a multi-objective optimization problem (MOP). In the literature, evolutionary multi-objective optimization algorithms (EMOAs), such as NSGA-II and MOEA/D, are reportedly effective to solve problems with two or three objective functions, but poor preformance are observed when they are used to address MaOPs. Algorithmic improvements to EMOAs are therefore needed for MaOPs. To date, the algorithmic improvements are usually made from three aspects: (1) Convergence acceleration: usually new dominance theory or modified Pareto dominance relationship is developed in this type of algorithms. See, for example, [51-53], can be significantly enhanced by different convergence enhancement mechanisms. (2) Decomposition: This type of algorithms works by decomposing the multi-objective problem into several single-objective problems. Typical examples of such algorithms include weight setting [54], decomposition mechanism [55, 56], generation of offspring and environment selection [57]. (3) Evaluation criteria: This type of algorithms improves the convergence and diversity of Pareto solutions by developing appropriate performance criteria [58].

In this section, inspired by the MOEA/D mechanism, we improved the decomposition mechanism of the algorithm and the population updating strategy, and proposed the MOIBA/AD algorithm.

In the literature, MOEA/D is one of the most widely used algorithm for solving MaOPs. The main idea is to combine multiple objectives into one single objective and solve this proxy objective to obtain the Pareto front. The three commonly used aggregation methods are weight sum approach, Tchebycheff approach [59] and penalty-based boundary intersection approach [60]. One major limitation of the three methods is that the time complexity of problem solving increases linearly (for the former method) or non-linearly (for the latter two methods) with the number of objective functions, thereby limiting the size of problems that can be solved.

For MOPs, the widely used methods for individual selection include entropy-based methods [61], grid-based methods [62] and cluster-based methods [63], amont others. Although these methods can produce desirable distributions of population on some specific multi-objective problems, they generally do not produce an uniform distribution of Pareto front on a general MaOP.

In almost all existing MOEA/D algorithms, the simulated binary crossover operator (SBX) is used for population update. However, this method has two limitations: 1) SBX reduces the population diversity. In the early stage of the search, the search space is usually small. For multi-objective problems with complex Pareto sets, SBX would probably eliminate the diversity of the population. 2) SBX often leads to inferior solutions to MOEAs.

4.1 Main motivation and innovation

Our proposed algorithm differs from existing algorithms in three aspects.

Firstly, we decompose the target space based on angle decomposition. Based on uniformly distributed unit vectors, the target space is divided into multiple sub-regions, and only the non-dominated individuals having the smallest angle with the weight vector will be retained. Through this update strategy, the Pareto front is decomposed into several fronts, and the PF of each subregion is linearly approximated without any aggregation function. This is different from the aggregation strategy of multiple strategies in traditional MOEA/D and helps speed up the calculation.

Secondly, we developed a strategy to improve the distribution of the population. We combine the individual dominance relationship with our proposed population selection and replacement strategy to obtain offspring individuals with greater diversity.

Lastly, we improve the traditional MOEA/D by taking advantage of the IBA, which combines local search with global search. Moreover, we apply random perturbations to the speed and position updating in BA, so as to prevent the search from getting trapped in local optima.

To sum up, the MOIBA/D continues the framework of the traditional MOEA/D algorithm, that is, the solution space is decomposed by appropriate decomposition strategy, and the required PF solution set is finally obtained by constantly updating the initial population (selection and replacement strategies) in the solution space. On the basis of this framework, in order to reduce the time complexity of the algorithm, we adopt the strategy of angle decomposition to the solution space, and on this basis, we put forward the strategy of population selection and updating of solution space based on angle decomposition. The IBA proposed in this paper is used in the specific implementation of the population individual optimization process.

4.2 Introduction of MOIBA/AD strategy

4.2.1 Method of target space decomposition and population classification

Given a set of uniformly distributed direction vectors λ1, λ2, ⋯, λ, optimization target dimension is M, the objective space of the multi-objective optimization problem can be decomposed into according to formula (13) and (14), and the population Pop is divided into by Eq (15). where F(x) represent fitness functions, = (Z1, Z2, ⋯, Z) represent a reference point, Z = min{f(x)|x∈Pop}, f(x) is the j optimization objective function.

In this way, the target space was divided by angle decomposition and the populations were classified, which was ready for the next update and screening. Taking the two-dimensional target space as an example, the spatial decomposition diagram was shown in Fig 4.

Fig 4

An illustration of the 2-dimensional target space decomposition.

4.2.2 Population selection and replacement strategy

is updated after the population classification is completed, in order to keep only one individual in each that is well divided.

If the sub-region contains an individual in , (1) if only one individual is contained, the individual is retained; (2) if there were more than one individual, all the individuals were subjected to non-dominated sorting [64]. The non-dominated individual was preferably selected, and then the individual with the smallest angle with the non-dominated individual was selected according to Eq (16). If the sub-region did not contain the individuals in , the individual corresponding to the target vector with the smallest angle with λ was selected from the current population and retained as the representative individual in , while the individuals in the remaining region were discarded. The population update diagram is shown in Fig 4. The solution space is decomposed into N regions by N direction vectors, and the individuals in the solution space is also classified into N parts. In the subspaces , there are multiple individuals in and , but no individual in . In this case, and make individual selection through non-dominated sorting, and the individuals are reserved as s1 and s2 marked in red. In , the individual corresponding to the target vector with the smallest angle of λ2 is reserved, and the reserved individual is the s2 marked in red in .

4.2.3 Optimization method of IBA

The traditional bat algorithm is inspired by the echolocation principle of bats looking for targets [65]. It controls the bat position and flight speed by controlling the frequency, loudness and rate of the transmitted pulses of the bats, which well combines the global optimization with the local optimization. In order to further improve the optimization ability of the algorithm and maintain the population distribution, the following improvements are made on the basis of the traditional bat algorithm.

In BA, the frequency adjustment formula for bats is: where f = 1, 2, ⋯, N, represents an individual of a bat population.

Improvement of global optimization for bats.

The improved speed update formula is:

In the iterative process of the algorithm, whether the position of the current individual is better than the average value is judged according to the average value of all individual positions in the current bat population; if the position of the current individual is better than the average value, the updating speed of the first formula in formula (18) is adopted. Otherwise, the updating speed of the second formula is adopted. In the first formula, the flying speed of the bat is added with a disturbance coefficient ω, ω = ωmax − (ωmax − ωmin) · iter/itermax, where ωmax, ωmin represent the maximum and minimum values of the disturbance coefficient, and iter, itermax represent the current number of iterations and the maximum number of iterations of the bat population. The velocity variable of the previous generation is multiplied by the velocity coefficient ω to disturb the individual velocity, so as to prevent it from flying to the local optimum and avoid the algorithm from falling into the local optimum. When the individual position of the bat is worse than the population average, the original formula is still used to update the speed to accelerate the convergence of the algorithm, which can help the individual in the disadvantaged position to fly to the optimal position quickly.

The improved location update formula is:

Compare with the traditional position update formula of BA algorithm, the periodic function is introduced to replace the fixed coefficient value, which is the amplitude of the trigonometric function, to spread the fluctuations of the bat’s acoustic frequency to the changes of the bat’s position, so that the bat’s position maintains the ability of continuous update, and thus maintains the diversity and distribution of the population.

Improvement of local optimization for bats.

In the local optimization method, the random Gaussian disturbance is added, where ε is of the same order as the current population matrix , ε ∼ N(0, 1). To avoid excessive fluctuation, represents the average loudness of all bat individuals currently emitted, and is used to adjust the search range of ε.

Local optimization and global optimization are linked by pulse transmit loudness and rate: where α, γ represent that rate of change of loudness and the rate of change of pulse emission, both of which are constant. Obviously, ∀0 < α < 1, γ > 0, t → ∞, , .

4.2.4 Population selection strategy

For individual selection of the updated population of each generation, in addition to the updating of the individual position and speed according to the updating strategy in Section 4.2.3, it is necessary to select and trim individuals according to the population size in the updated population.

In order to better explore the sparse region, the individual congestion generated by the local optimization method according to Eq (19) is sorted, and then the individual selection is conducted. In the global optimization process, the initial selection is firstly performed according to Eqs (16)–(19). Next, for the preservation of population distribution, the following selection method based on unit neighborhood vector is proposed:

Let , where is calculated as: where is the g-dimensional component of λ, and (f = 1, 2, ⋯, N) is the T unit neighborhood vectors of unit vector λ. In this way, the individual with the smallest angle with λ is selected by Eq (23), which greatly improves the population distribution. A five-dimensional target was taken as an example, and the selection diagram of the population was shown in Fig 5. We first select s1, s2, s3 and s5 for the regions , , and according to (16)–(19). In region , there is no candidate individuals available for selection. In this case, we used the proposed method based on unit neighborhood vector to select s4 as the retained individual for regions .

Fig 5

A illustration of the individual selection.

4.3 MOIBA/AD framework for MOOCTS

4.3.1 Population coding scheme and spatial transformation strategy for solving MOOCTS

In the introduction, many papers use binary code and Gray code when studying this kind of problem, but the algorithm of solving model in this paper is completely combined with BA algorithm, and its population search and optimization operation belong to the category of real number operation, so this paper adopts real number code. Since the sequence number set of the combat task is I = {1, 2, ⋯, n} and the sequence number set of the combat cluster is K = {1, 2, ⋯, m}, the code vectors are X = (x, x, ⋯, x), X = (x, x, ⋯, x) respectively.

Regarding population initialization in Algorithm 1, we propose a strategy to transform discrete space into continuous space, rather than apply random initialization, which is the most commonly used method for creating the initial population. The basic idea of our strategy is to transform the binary values in discrete space into some intermediate variables, which will then be converted into real values in continuous space. At the heart of our strategy is determining an appropriate function for spatial transformation. Our strategy consists of three steps, which are shown below.

Step1: After population initialization, we choose the sigmoid function as the function for spatial transformation and calculate a value for each individual by applying this function. The sigmoid function is given by

Step2: Save all populations in the intermediate variable Pop.

Step3: Transform the intermediate variables into binary values according to .

With this procedure, the continuous space is transformed into discrete space.

4.3.2 Algorithm step

Under the above strategy, the steps of MOIBA/AD algorithm proposed in this paper are as Algorithm 1. In the framework of MOIBA/AD algorithm, the operational task scheduling model, ITFN and other factors are added into the process of problem solving. The whole MOOCTS model solving process is shown in Fig 6.

Fig 6

A illustration of the whole MOOCTS model solving process.

4.3.3 Computational complexity of MOIBA/AD

In order to further illustrate the efficiency of the algorithm, the analysis of MOIBA/AD step shows that its time complexity depends on the construction and generation method of Pareto optimal solution set, which mainly comes from the updating of individuals and reference points in Step2. In Step1.4, the euclidean distance between direction vectors λ1, λ2, ⋯, λ is calculated, the complexity is O(N2), and the complexity in the initialization process of Step1 is kO(N ⋅ M) + O(N2). In the process of population updating, the complexity of updating individuals and reference points is O(D ⋅ N ⋅ M ⋅ T) by using the improved BA, where D is the dimension in IBA, T is the number of unit neighbor vectors, as T ≤ N, M ≤ N, the time complexity of MOIBA/AD is O(D ⋅ N ⋅ M ⋅ T), which indicates that MOIBA/AD is computationally efficient.

Algorithm 1 MOIBA/AD

Input: MOP, set of direction vectors {λ1, λ2, ⋯, λ}, the reference point Z, the maximum number of iterations itermax, initial loudness of pulse A0, initial pulse emission rate r0, change of pulse loudness rate α, change of pulse transmission frequency γ, initial pulse frequency f0, maximum pulse frequency f, minimum pulse frequency f, individual initial velocity v0, maximum and minimum values of disturbance coefficient of individual velocity ωmax,ωmin, amplitude of position disturbance function A.

Output: optimal EP;

1: Step 1 initializing;

2: Step1.1 set an external noninferior solution EP = ∅;

3: Step1.2 generating N uniformly distributed direction vectors λ1, λ2, ⋯, λ, according to the population size N, number of iterations iter = 0;

4: Step1.3 generate initial population Pop(0), initial speed v0;

5: Step1.4 calculating the euclidean distance between λ and λ, i, j = 1, 2, ⋯, N, determining T unit neighbor vectors ;

6: Step1.5 determine an initial reference point , where , j = 1, 2, ⋯, M;

7: Step2 population update;

8: Step2.1 update the individuals in Pop(iter) according to the update strategy of 3.2.2,so that only one individual is reserved in each divided P, and finally N individuals are selected;

9: Step2.2 updating the progeny population according to the formula (17)–(20), storing the updated N individuals in EP,Pop(iter) = Pop(iter) ∪ EP, and updating the reference point;

10: Step2.3 select N individuals from Pop(iter) as Pop(iter+ 1) according to the selection strategy in 3.2.4, iter = iter + 1, and update EP;

11: Step3 jud whether that iter = itermax is satisfy, if so, terminating the algorithm and outputting EP; if not, return to Step2.

5 Experiment and result analysis

The following is a practical example to solve the problem of cooperative time scheduling in land offensive operations. The algorithm proposed in this paper with some compared algorithms are used for simulation verification and result analysis. The simulation environment is CPU i7–8850H, 16.0GB RAM, operating system Windows10.

5.1 Simulation example

In this paper, based on Experiment 7 [66] in the template library A2C2 of US military combat experiment scene, the MOOCTS model established in this paper is solved and analyzed. Experiment 7 in A2C2 was a joint landing operation, its combat goal is to occupy airports and ports and clearing the way for subsequent landing troops. There are a total of 18 combat tasks (i = 1, 2, …, 18), the task network diagram is shown in Fig 7.

Fig 7

Network diagram of scenario task for simulation example.

There are 20 combat clusters (k = 1, 2, …, 20), represented by the number P1 − P20, there are 8 combat capabilities in this operation, namely AAW, ASUW, ASW, GASLT, FIRE, ARM, MINE and DES. Due to the large amount of data in the support information S1 Excel, we enumerate the combat capability of each combat cluster and the cluster capability data required for completing various combat tasks in detail. In addition, according to different combat missions, three groups of required times are set for comparative experiments. The comparison time data are shown in Table 3.

Table 3

The comparison time data of three groups of required times set for comparative experiments.

Instance	T1	T2	T3	T4	T5	T6	T7	T8	T9	T10	T11	T12	T13	T14	T15	T16	T17	T18
PT1	30	30	10	10	10	10	10	10	10	10	10	10	20	20	15	15	10	20
PT2	35	35	15	15	15	15	15	15	15	15	15	15	25	25	20	20	15	25
PT3	40	40	20	20	20	20	20	20	20	20	20	20	30	30	25	25	20	30

PT1, PT2 and PT3 respectively represent Processing Time1, Processing Time2 and Processing Time3. The three different sets of comparative data will change the scale of the comparative experiment.

The combat task durations (PT1, PT2, PT3) corresponding to three groups of combat missions with different intensities are set as the comparison conditions of the algorithm simulation results. , let , , where rand(1) is a function that generates evenly distributed random numbers between 0 and 1.

The combat flexible time is set as follows: in order to ensure sufficient combat coordination time, the SOUTH BEACH mission (T8) can only start 5 hours after the start of the NORTH BEACH mission (T7) at the earliest and must start 2 hours before the end of the T7 mission at the latest, that is , similarly, let .

5.2 Comparison algorithm

In this paper, four comparison algorithms are compared with the MOIBA/D algorithms, that is NSGA-II, MOEA/d, MOPSO and MOBA algorithm. The MOCCTS model is incorporated into the solution frameworks of the five algorithms, and the population coding mode of the comparative algorithm was consistent with the MOIBA/D algorithm. The algorithm initialization parameters were set as shown in Table 4.

Table 4

Parameter settings for 5 algorithms.

Algorithm	Parameter setting
NSGA-II	Population size: N = 1000, Mating pool size: MP = 600, Termination algebra: itermax = 5000, Crossover probability: Pc = 0.8, Variation probability: Pm = 0.2, Crossover distribution index: ηc = 1, Variation distribution index: ηm = 1;
MOEA/D	Weight vector number: N = 1000, Neighbor vector number: T = 20, Termination algebra: itermax = 5000, Crossover distribution index: ηc = 1, Variation distribution index: ηm = 1;
MOPSO	Weight vector number: N = 1000, Neighbor vector number: T = 20, Termination algebra: itermax = 5000, Crossover distribution index: ηc = 1, Variation distribution index: ηm = 1;
MOBA & MOIBA/AD	population size N = 1000, maximum iteration number itermax = 5000, initial pulse loudness A0 = 0, initial pulse emission rate r0 = 0, pulse loudness change rate and pulse emission rate change rate α = 0.9, γ = 0.05, pulse initial frequency f0 = 0, maximum and minimum pulse frequency fmax = 10, fmin = 0, individual initial velocity v0 = 0, maximum and minimum individual velocity disturbance coefficient ωmax = 1.2, ωmin = 0.1, amplitude of individual position disturbance function Am = 2.

5.3 Performance metrics

Due to the complexity of the model and the unknown results, this paper selects the following four indicators as the performance evaluation indicators of the algorithm.

Hypervolume (HV) [67]. The calculation formula is as follows: where λ is the Lebesgue measure, v is the hypervolume composed of the reference point Z and the non-dominated individual Pop(itermax), and S is the non-dominated solution set. The calculation process does not need to know the Pareto optimal frontier, so it is very practical. Among many evaluation indicators, its monotonicity is good, and the larger the value, the better the convergence and distribution of the algorithm.

Inverted generational distance (IGD) [68]. The calculation formula is as follows: where , represents the minimum euclidean distance from a point on the optimal Pareto front to an individual in the final solution set. IGD can not only reflect that convergence of multi-objective evolutionary algorithm, but also reflect the distribution and universality of the solution set. The smaller the value is, the better the convergence, distribution and universality of the algorithm will be.

Comprehensive performance evaluation index (MSS). This indicator is a further expansion of the statistical distribution of IGD indicators, and the calculation formula as: where M is the number of optimization objectives, representing the IGD value of the evaluation algorithm on the m optimization objective, and μ, δ is the mean and variance of the IGD value of the evaluation algorithm on the m optimization objective. Its value can not only reflect the convergence, distribution and universality of the algorithm reflected by IGD, but also reflect the statistical law of the final solution set.

Number of dominated solutions. The index counts the number of dominated solutions in the final solution set of MOCTS models solved by different algorithms, and the numerical results obtained are used to express the performance of the algorithms in solving the model. Obviously, the more dominant solutions, the better the performance of the algorithm will be.

5.4 Result analysis

First of all, in order to avoid the influence of random error on the experimental results, in three calculation examples, five algorithms were run independently for 20 times. The average value, best value and worst value obtained by statistics are listed in the table below. The bold part was the optimal value of the probability degree of operational coordinative time under various indexes.

As shown in Table 5, the results of MOIBA/AD algorithm in different numerical examples are generally better than the four comparison algorithms. Although the optimal values are the same as those of NSGA-II algorithm and MOEA/D algorithm when solving numerical Instances 1 and 2, MOIBA/AD algorithm is obviously more reliable from the point of view of mean value and worst value. With the expansion of numerical instances, the advantages of MOIBA/AD algorithm gradually appear. In order to intuitively reflect the scheduling results of the five algorithms on the operational coordinative time, we draw the Gantt chart of the scheduling results of the five algorithms in the operational coordinative time.

Table 5

Comparison of simulation results of all algorithms.

Instance	Algorithm	Average value	Best value	Worst value
PT1	NSGA-II	(158,166,174)	(150,160,170)	(168,176,184)
	MOEA/D	(161,165,172)	(150,160,170)	(162,170,181)
	MOPSO	(162,166,177)	(158,166,174)	(166,177,186)
	MOBA	(164,170,179)	(158,166,174)	(164,172,184)
	MOIBA/AD	(158,166,174)	(150,160,170)	(162,170,181)
PT2	NSGA-II	(230,246,254)	(220,230,240)	(232,240,251)
	MOEA/D	(228,236,244)	(220,230,240)	(243,255,269)
	MOPSO	(232,249,258)	(228,239,251)	(242,256,271)
	MOBA	(234,251,266)	(229,238,247)	(240,250,261)
	MOIBA/AD	(228,236,244)	(220,230,240)	(232,240,251)
PT3	NSGA-II	(300,313,331)	(297,306,316)	(309,319,329)
	MOEA/D	(305,311,336)	(296,307,319)	(306,318,326)
	MOPSO	(306,313,334)	(297,310,318)	(305,313,326)
	MOBA	(301,311,330)	(299,308,317)	(303,311,321)
	MOIBA/AD	(298,306,314)	(290,300,310)	(302,310,321)

In Fig 8, the durations of the same combat mission is indicated with the same color. It can be seen from the collaborative scheduling results that the scheduling time of MOIBA/AD algorithms is the least, while that of MOBA algorithms is the most, and the scheduling time of the other three algorithms is similar. Therefore, under the framework of MOIBA/AD algorithm, the performance of solving MOCCTS model is optimal.

Fig 8

Gantt chart of operational coordination time planning results of five algorithms under simulation background.

The comparison results of HV, IGD, MSS and the number of dominated solutions indexes obtained by simulation are shown in Table 6. It can be concluded that, on the whole, the performance of MOIBA/AD is superior to that of the comparison algorithm, and the advantages of MOIBA/AD are more obvious when the scale of the MOOCTS model is expanded. Specifically, on the index HV and the index MSS, the performance of MOIBA/AD is better than that of the comparison algorithm, and the effect of the algorithm is still excellent after the model scale is expanded. On the indicator IGD, the algorithm still maintains a good effect on the whole, except that the performance of the algorithm is similar to that of NSGA-II after the scale is enlarged. It can also be concluded that when NSGA-II solves large-scale multi-objective problems with sufficient number of iterations, the algorithm has excellent performance. In summary, MOIBA/AD obtains five optimal values in six groups of data under four indicators, which is very effective in solving the MOOCTS model.

Table 6

The values of HV, IGD, MSS and the number of non-dominated solutions corresponding to each of the comparison algorithms in solving three different-scale cooperative time scheduling problems.

Parameter	Instance	MOBA	MOPSO	MOEA/D	NSGA-II	MOIBA/AD
HV	PT1	0.1143−(0.0916−)	0.4892 +(0.3793 +)	0.1528 −(0.1134 +)	0.1662 (0.1091)	0.5112 +(0.3823 +)
	PT2	0.4861−(0.4555−)	0.5392 +(0.4491 +)	0.4966 −(0.4135 −)	0.5026 (0.4756)	0.5613 +(0.4624 +)
	PT3	0.4987−(0.4235−)	0.5045 +(0.4589 +)	0.4817 −(0.4145 −)	0.5137 (0.4512)	0.5823 +(0.4714 +)
IGD	PT1	1.323−(0.967−)	1.256 −(1.002 −)	0.088 +(0.061 +)	0.553 (0.144)	0.037 +(0.025 +)
	PT2	0.084−(0.012−)	0.107 −(0.044 −)	0.116 −(0.192 −)	0.053 (0.008)	0.054 =(0.013 =)
	PT3	1.023−(0.614−)	1.151 −(1.101 −)	0.188 +(0.161 +)	0.603 (0.214)	0.437 +(0.326 +)
MSS	PT1	10.231+	−1.761 +	1.433 +	23.899	−6.661 +
	PT2	17.891+	−18.231 +	6.963 +	24.691	−11.781 +
	PT3	17.232+	11.766 +	6.413 +	25.781	−16.002 +
No. of dominated solution	PT1	311−	321−	320−	348	347=
	PT2	291+	311+	330+	290	462+
	PT3	280−	348−	383+	360	382+
+/−/≈		4/8/0(4/8/0)	7/5/0(7/5/0)	7/5/0(8/4/0)		10/0/2(10/0/2)

“+ /−/ =” represents that the test algorithm is superior, equivalent and inferior to the comparison algorithm, respectively. The optimal value of each index data is expressed in bold, with the data outside brackets as the optimal value of the index and the data in brackets as the average value of the index.

The average convergence curves of indicators HV and IGD on three different scale models are shown in Figs 9 and 10. It can be seen that when solving the MOOCTS model with a small scale, the performance of the four comparison algorithms is unstable. The convergence rates of MOBA, NSGA-II and MOEA/D in the HV index are too slow, and the index still fails to converge when the maximum number of iterations is reached. The convergence rates of MOBA, NSGA-II and MOPSO in the IGD index are too slow, which are more than 2000 generations later than the comparison algorithm. However, MOIBA/AD shows good index convergence when solving the MOOCTS model, regardless of the large or small scale.

Fig 9

Average convergence curves of the HV indexes on the model of three different scales, obtained by a large-scale model with 10000 iterations (a)-(c), a small-scale model with 5000 iterations (d)-(f).

Fig 10

Average convergence curves of the IGD indexes on the model of three different scales, obtained by a large-scale model with 10000 iterations (a)-(c), a small-scale model with 5000 iterations (d)-(f).

In the simulation process, when the five algorithms solve the MOOCTS models of three different sizes, PF of optimization objective and illustration of parallel coordinates of the non-dominated fronts on 3 different scale MOOCTS models are shown in Fig 11. It can be seen that the diversity and convergence of MOIBA/AD are better than those of MOBA, NSGA-II, MOEA/D and MOPSO when solving the small-scale MOOCTS model. NSGA-II has poor diversity and convergence, MOEA/D has good diversity but poor convergence, and MOPSO has poor diversity but poor convergence; when solving large-scale MOOCTS models, the diversity and convergence of MOIBA/AD are better than those of MOBA, NSGA-II, MOEA/D and MOPSO. However, the diversity and convergence of the remaining comparison algorithms are greatly improved.

Fig 11

PF of optimization objective and illustration of parallel coordinates of the non-dominated fronts on three different scale ECMJTA models, obtained by small-scale model with 5000 iterations (a)-(c), large-scale model with 10000 iterations (d)-(f).

6 Conclusion and future work

In modern war, scheduling of operational coordinative time is one of the most important tasks in operational coordination. In this study, in order to solve this problem, we established the multi-objective operational cooperative time scheduling (MOCCTS) model. This problem is a complex multi-objective combinatorial optimization problem with the characteristics of multi-decision variables and multi-objective functions.In this paper, we propose a MOIBA/AD algorithm to solve this problem. This algorithm improves the MOEA/D by changing the aggregate function strategy used by MOEA/D to the strategy using angle space decomposition, thus greatly reducing the computational complexity. Moreover, the population renewal strategy is improved through the optimization of the IBA, considering that the IBA is not easy to enter the local optimum, the population renewal effect of the algorithm is effectively improved. Finally, MOIBA/AD is compared with four multi-objective evolutionary algorithms such as MOBA, NSGA-II, MOEA/D and MOPSO to solve three MOOCTS models with different scales. The comparison of performance indexes indicates that MOIBA/AD has good population convergence and distribution.

In conclusion, MOIBA/AD can effectively and efficiently solve the MOOCTS model of different sizes, and provide decision recommendations for the commander in the combat process.The MOOCTS problem proposed in this paper is actually an extension of the multi-task time scheduling problem, which only schedules the cooperative time in the course of combat. However, we do not consider the scheduling the remaining operational factors (i.e. operational resources, electromagnetic spectrum.) in a real and dynamic environment. Besides, in the process of solving the cooperative time sequence of combat missions, the method proposed in this paper first needs to analyze the time network among combat tasks. When the scale of combat mission expands, this step will produce great time complexity, which will greatly affect the practical application.

In the further, on the one hand, we need to further expand the MOOCTS model, take into account the factors such as combat support resource constraints and electromagnetic spectrum resource constraints, and add the comprehensive analysis of node factors on the basis of the directed graph processing of combat sequence, so that the model is closer to the actual combat. On the other hand, we will explore the algorithm of solving large-scale directed graphs and the processing multiple directed graphs in parallel in the model initialization, which can be applied to large-scale operations. In addition, the role of MOIBA/AD in solving other combinatorial optimization problems is also studied, and it is applied to more fields and problems.

Simulation data sheet.

Combat capability and cluster capability simulation data sheet of each combat cluster required for completing various combat tasks.

(XLSX)

Click here for additional data file.

Simulation process of MOIBA/AD algorithm.

(PDF)

Click here for additional data file.