Combining an Extended SMAA-2 Method with Integer Linear Programming for Task Assignment of Multi-UCAV under Multiple Uncertainties

Abstract

Uncertainty should be taken into account when establishing multiobjective task assignment models for multiple unmanned combat aerial vehicles (UCAVs) due to errors in the target information acquired by sensors, implicit preferences of the commander for operational objectives, and partially known weights of sensors. In this paper, we extend the stochastic multicriteria acceptability analysis-2 (SMAA-2) method and combine it with integer linear programming to achieve multiobjective task assignment for multi-UCAV under multiple uncertainties. We first represent the uncertain target information as normal distribution interval numbers so that the values of criteria (operational objectives) concerned can be computed based on the weighted arithmetic averaging operator. Thus, we obtain multiple criteria value matrices for each UCAV. Then, we propose a novel aggregation method to generate the final criteria value matrix based on which the holistic acceptability indices are computed by the extended SMAA-2 method. On this basis, we convert the task assignment model with uncertain parameters into an integer linear programming model without uncertainty so as to implement task assignment using the integer linear programming method. Finally, we conduct a case study and demonstrate the feasibility of the proposed method in solving the multiobjective task assignment problem multi-UCAV under multiple uncertainties.

1. Introduction

Unmanned combat aerial vehicles (UCAVs), compared with airborne weapons and airborne platforms (e.g., cruise missiles, manned fighters, and bombers), have demonstrated advantages in tactical flexibility and multitasking ability due to having zero casualties, excellent stealth performance, high flight height, long endurance, multiple hard points for missiles, and a relatively low life-cycle cost [1]. Benefiting from the rapid development of control theory and intelligent technology, UCAVs are increasingly likely to become the frontline fighter aircrafts in modern or future wars [2]. As one of the cores of mission planning, cooperative task assignment for UCAVs has been an emerging issue in terms of its application in attacking ground targets. The purpose of task assignment is to match tasks to available UCAVs without conflicts under the constraints of the battlefield environment, operational capabilities of the UCAVs and weapons on board, and target types in order to obtain the maximum operational benefit, such as maximizing damage to hostile targets and minimizing the user’s own consumption.

With regard to research on task assignment, most of the existing methods assume that the parameters related to mission scenarios are deterministic. These methods are mainly divided into two categories: methods based on exact algorithms and methods based on heuristic (metaheuristic) algorithms [3]. Exact algorithms are mostly based on the mathematical properties of task assignment models and attempt to obtain the optimal solution for given objectives. They include the mixed integer linear programming (MILP) formulation [4,5,6], the branch and bound tree search algorithm [7,8], and the dynamic programming algorithm [9] among others. Most research to date on the application of heuristic (metaheuristic) algorithms have either extracted some specific search rules based on the properties of the problem to obtain optimal or suboptimal solutions rapidly or introduced some local search mechanism to the basic algorithm framework to improve the solution quality, such as tabu search algorithms [10], auction algorithms [11], genetic algorithms [12,13,14,15], ant colony algorithms [16], and particle swarm optimization [17,18]. The algorithms mentioned above—whether belonging to exact algorithms or heuristic (metaheuristic) algorithms—have demonstrated the ability to provide optimal or suboptimal solutions for task assignment problems of UAVs or UCAVs in various mission scenarios but may become difficult to apply when there are uncertain parameters related to mission scenarios.

From what we know, uncertainties in the task assignment model for multi-UCAV exist in at least three aspects: (1) the target information obtained by various sensors, (2) the weights of operational objectives, and (3) the weights of these sensors. Considering that sensors may be affected by the dynamic battlefield environment or the antireconnaissance measures taken by hostile forces in the reconnaissance process, it is reasonable to believe that there are errors in the target information obtained by various sensors, indicating that the target information is uncertain. When there are multiple operational objectives to consider, the commander needs to weigh the importance of these objectives but can usually only provide the relative importance of operational objectives instead of the exact weights of these objectives. For example, the commander may, at most, express a clear preference for operational benefit over operational cost but would not provide the exact weights of the two operational objectives. Therefore, the weights of operational objectives are also thought to be uncertain. In addition, the uncertainty contained in the weights of sensors also needs to be considered. In this paper, the weight of a sensor is obtained by integrating its subjective weight and objective weight in a certain way. The subjective weight is usually determined in advance by experts according to the performance parameters of the sensor, while the objective weight is heavily dependent on the accuracy and credibility of the target information provided by the sensor. As the target information is uncertain, the objective weight of the sensor is unknown, which indicates that the integrated weight of the sensor is partially known.

In order to reflect the combat scenarios as accurately as possible, the above three aspects of uncertainties need to be taken into account when formulating the task assignment problem. However, in this case, the problem would become extremely complicated and can neither be solved by exact algorithms nor by heuristic (metaheuristic) algorithms. Therefore, in order to address the challenge, it is necessary and worthwhile to propose new methods. In this way, the optimal or suboptimal task assignment scheme under severe uncertainties can be found.

2. Related Work

To our knowledge, Mataric et al. first studied task assignment for multirobots in uncertain environments in 2003 [19]. The uncertainty in their research mainly included sensor noise and the randomness of the moving direction of robots. To simplify, they used a random number from a normal distribution to represent the sensor noise and a finite but small probability that robots would move in a random direction instead of the intended direction at each time step to model the randomness of the movement. Since then, various methods have been proposed to solve task assignment problems under different types of uncertainty.

In general, uncertainty is mathematically described by one or a combination of the following theories: probability theory, fuzzy set theory, rough set theory, and interval set theory [20]. In the existing literature on task assignment under uncertainty, as far as we know, probability theory has been the most widely used. The simplest way to model uncertainty based on probability theory is to assume that an uncertain parameter follows a specific distribution. Based on this assumption, other methods can be applied to task assignment. Mataric et al. assumed that sensor noise followed a normal distribution and that each robot had a small probability of moving in other directions than the intended direction. The authors then adopted the auction algorithm to assign tasks to robots [19]. Kim et al. used two certain probability values to measure the likelihood that a UAV would sense targets within its sensing range and the likelihood that a UAV would fail during the mission. The authors then proposed a distributed task assignment method based on resource welfare [21]. Jia et al. assumed that the flight velocity of a UAV was a stochastic variable that corresponded to a normal distribution with alterable variance according to the actual flight distance between each pair of targets. The authors then solved the task assignment problem in the framework of a two-stage stochastic programming model using a genetic algorithm [3]. A similar way of dealing with uncertain parameters can be found in Reference [22], where the task duration time was considered to be subject to a uniform distribution.

In multistage task assignment problems, the Markov decision process (MDP) or partially observable Markov decision process (POMDP) is applied based on the assumption that the uncertain parameters that are closely related to the state transition are subject to either predefined distributions or predefined probability values, as can be seen in References [23,24,25,26,27]. For example, in Reference [24], it was assumed that each target could move to any of its neighboring cells at each time step with the same probability and that a UAV could detect a target with a predefined probability. Based on the above assumptions, the possible states of UAVs and targets, the actions that UAVs could take, and the corresponding rewards of UAVs could be obtained in turn. Then, the expected cumulative reward earned during finite time steps could be computed. Thus, the task assignment scheme could be acquired by maximizing the expected cumulative reward.

In addition to probability theory, interval set theory has also been applied to task assignment problems. Liu et al. proposed an interval Hungarian algorithm for task assignment of UAVs when the utility of a UAV performing a task is an interval number. The proposed algorithm was an extension of the classic Hungarian algorithm and was able to compute the maximum interval of deviation for each entry in the assignment matrix, which would retain the same optimal assignment [28,29]. Hooshangi et al. used the interval VIKOR method to sort tasks, the characteristics of which were specified as interval numbers. They then applied an auction algorithm to assign tasks to agents [20]. The task assignment problem in this research was, to some extent, regarded as a decision problem with multiple uncertain criteria, and the simulation results demonstrated that the proposed method was able to solve task assignment problems with more uncertainty than the other methods mentioned above. Despite this, when the criteria values and the corresponding weights are both uncertain, the interval VIKOR method would become unavailable.

The methods introduced above have been proven to be effective in solving the task assignment problem where uncertainty exists either in the action of agents (UAVs) or in task rewards. However, when multiple uncertainties are simultaneously contained in the task assignment problem, as described in Section 1, none of these methods can be applied. Therefore, we need to put forward new methods to solve the problem. As a multicriteria decision support method, the stochastic multicriteria acceptability analysis (SMAA) is able to solve problems with unavailable criteria weight information and uncertain criteria values [30]. Hu et al. investigated the application of the SMAA-2 method, a variant of the SMAA method, in task assignment of multi-UAV when both the target information and the weights of operational objectives were uncertain. The proposed method demonstrated the availability of the SMAA-2 method in solving the task assignment problem that contained uncertainties from two aspects. However, the method was limited in its application [31] because the premise of applying the SMAA-2 method was that there could only be one criteria value matrix: the evaluation of all targets in terms of all concerned operational objectives. Note that in this paper, the criteria in the SMAA-2 method are equivalent to the operational objectives. In the task assignment model described in Section 1, there are various sensors and each sensor has the ability to acquire battlefield information, so there are multiple criteria value matrices. In this situation, if the SMAA-2 method is to be applied, these matrices would be first aggregated to generate a final criteria value matrix. Considering that the weights of sensors are partially known, the aggregation process cannot be directly performed. In summary, the original SMAA-2 method is difficult to use in such a situation, so it is in urgent need of improvement.

Inspired by References [30,31], we extend the original SMAA-2 method in this research by integrating a novel aggregation method so that multiple criteria value matrices can be aggregated to generate the final criteria value matrix based on which the SMAA-2 method can be applied. Then we combine the extended SMAA-2 method with integer linear programming (ILP) to solve task assignment problems for multi-UCAV where there are multiple sensors (or information sources) to provide uncertain information, describing multiple attributes of multiple targets. Furthermore, the weights of the criteria concerned are uncertain and the weight information of information sources is partially available. In subsequent sections, we will refer to the above combination method as ESMAA-2-ILP for convenience. The flow of the method is illustrated in Figure 1.

As can be seen from Figure 1, to achieve task assignment for multi-UCAV, we first formulate the task assignment problem with a variety of uncertain parameters to obtain the criteria formulas and then, based on the target information from K information sources, we construct K criteria value matrices for each UCAV through the above criteria formulas. For each UCAV, we apply the extended SMAA-2 method to compute the holistic acceptability index of attacking each target. Thus, we can obtain a holistic acceptability index matrix, which describes the holistic acceptability that each UCAV attacks each target. On this basis, according to the meaning of holistic acceptability indices, we convert the task assignment model established above into an integer linear programming model with an objective function to maximize the holistic acceptability index of the task assignment scheme. Finally, based on integer linear programming, we obtain the task assignment scheme. The main contributions of this research are summarized as follows:

• We propose a novel method to aggregate multiple criteria value matrices from different information sources (ISs) under the condition that the criteria weight information is unavailable and the weight information about ISs is partially available.
• We integrate the aggregation method into the original SMAA-2 method. Thus, we propose an extended SMAA-2 method.
• We propose to convert the task assignment model for multi-UCAV with a variety of uncertain parameters into an integer linear programming model to maximize the holistic acceptability index of the task assignment scheme according to the meaning of holistic acceptability indices, and then apply integer linear programming to achieve task assignment.

The rest of this paper is organized as follows. We formulate the task assignment problem with a variety of uncertain parameters in Section 3. In Section 4, we briefly introduce the original SMAA-2 method. Section 5 gives details of the extended SMAA-2 method. In Section 6, we convert the task assignment model established in Section 3 into an integer linear programming model. We examine the proposed method through a case study in Section 7. Section 8 concludes the paper.

3. Problem Formulation

Here, we consider a combat scenario where M UCAVs under the command of one command center are to attack N geographically dispersed ground targets with different values. Each UCAV is loaded with a certain number of the same kind of air-to-ground missiles to destroy targets with different lethality probabilities. Prior to the strike, K sensors, including reconnaissance satellites, reconnaissance aircrafts, and reconnaissance UAVs, conduct multiple reconnaissances on targets and then transmit acquired target information to the command center. Based on the received target information, the command center computes K criteria value matrices for each UCAV and then aggregates them to generate a final criteria value matrix according to which commanders are able to achieve task assignment for multi-UCAV. The operational process described above is briefly shown in Figure 2. The command center receives the target information transmitted by three sensors (one reconnaissance satellite, one reconnaissance aircraft, and one reconnaissance UAV) and then assigns four targets to four UCAVs.

We assume that each UCAV can only attack one target in one flight and that each target can only be attacked by one UCAV (for convenience, we assume M = N). Furthermore, we consider maximizing the damage to targets (i.e., operational benefit) while minimizing the flight consumption of UCAVs (i.e., operational cost) and then formulate the task assignment problem as follows:

$$\max F\left(x\right)={\omega }_B{O}_B+{\omega }_C{O}_C$$

(1)

$$O_B={\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N{x}_{ij}o\_b_{ij}^{\left(0\right)}}}$$

(2)

$$O_C=-{\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N{x}_{ij}o\_c_{ij}^{\left(0\right)}}}$$

(3)

subject to

$${\displaystyle \sum _{i=1}^M{x}_{ij}}=1$$

(4)

$${\displaystyle \sum _{j=1}^N{x}_{ij}}=1$$

(5)

$$x_{ij}\in \left\{0,1\right\}$$

(6)

$${\omega }_B+{\omega }_C=1,{\omega }_B,{\omega }_C\in \left[0,1\right]$$

(7)

$$o\_b_{ij}^{\left(0\right)}={\displaystyle \sum _{k=1}^K{w}_{k}o\_b_{ij}^{\left(k\right)}}$$

(8)

$$o\_c_{ij}^{\left(0\right)}={\displaystyle \sum _{k=1}^K{w}_{k}o\_c_{ij}^{\left(k\right)}}$$

(9)

$$o\_b_{ij}^{\left(k\right)}=v_j^{\left(k\right)}\left[1-{\left(1-x_{ij}{p}_{ij}\right)}^{q_i}\right]$$

(10)

$$o\_c_{ij}^{\left(k\right)}=d_{ij}^{\left(k\right)}{c}_{ij}^{fly}$$

(11)

$${\displaystyle \sum _{k=1}^K{w}_k}=1,w_k\in \left[0,1\right]$$

(12)

$$i=1,2,\cdots ,M;j=1,2,\cdots ,N;k=1,2,\cdots ,K$$

(13)

In the above model, the objective function, Equation (1), aims to maximize the combat reward, a weighted sum of operational benefit $O_B$ and operational cost $O_C$. The corresponding weights ${\omega }_B$ and ${\omega }_C$ are determined by commanders and should satisfy constraint in Equation (7). The constraints expressed in Equations (4)–(6) guarantee a one-to-one relationship between the UCAVs and the targets. The constraints in Equations (8)–(12) define the method to aggregate criteria value matrices.

Solutions for the above model are relatively easy to find when all parameters are deterministic. However, due to the complexity of actual combat, many parameters in the model are uncertain. Considering errors in the target information obtained by sensors, the following parameters $v_j^{\left(k\right)}$ and $d_{ij}^{\left(k\right)}$ ($i=1,2,\cdots ,M;j=1,2,\cdots ,N;k=1,2,\cdots ,K$) are assumed to be uncertain. Moreover, in our research, criteria weight information is unavailable and sensor weight information is partially available, which means ${\omega }_B$, ${\omega }_C$ and $w_k$ ($k=1,2,\cdots ,K$) are also uncertain.

As can be seen, to solve the task assignment problem, we need to first aggregate K criteria value matrices under serve uncertainty, and then conduct the stochastic multicriteria acceptability analysis for each UCAV in terms of all targets. Therefore, we suggest extending the original SMAA-2 method by combining a novel aggregation method. To facilitate a better understanding of the extended SMAA-2 method, we give a brief introduction to the original SMAA-2 method in the next section.

4. The Stochastic Multicriteria Acceptability Analysis-2 (SMAA-2) Method

Stochastic multicriteria acceptability analysis (SMAA), first proposed by Lahdelma et al. in 1998 [30], is a multicriteria decision support method to solve discrete problems with unavailable criteria weight information and uncertain criteria values. Based on the exploration of the weight space, the method is able to describe the valuations that would make each alternative the preferred one. During this process, uncertain criteria values are always represented by probability distributions, thus the method can compute confidence factors to describe the reliability of the analysis. The SMAA-2 method, also proposed by Lahdelma et al., extends the original SMAA by analyzing the sets of weight vectors for any rank from best to worst for each decision alternative [32]. It is introduced as follows:

Here, we consider a multicriteria decision-making problem, including a set of N alternatives $Set\_alt=\left\{A_1,A_2,\cdots ,A_N\right\}$, from which one is to be chosen. All alternatives are evaluated in terms of n criteria $Set\_cri=\left\{C_1,C_2,\cdots ,C_n\right\}$. The method maps the ith alternative $A_i$ to a utility value $U_i$ by $U_i=U\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)$, where ${\mathsf{\xi }}_i=\left[{\xi }_{i1},{\xi }_{i2},\cdots ,{\xi }_{in}\right]$ denotes the criteria value vector of alternative $A_i$, $\mathsf{\omega }=\left[{\omega }_1,{\omega }_2,\cdots ,{\omega }_n\right]$ is a weight vector for each decision-maker to quantify their subjective preferences for different criteria, and $U(\cdot )$ is a utility aggregation function. In SMAA-2, any type of utility aggregation function is available [33]. In this paper, we use the utility additive form and then obtain the following:

$$U_i=U\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)={\displaystyle \sum _{j=1}^n{\omega }_j{u}_{ij}\left({\xi }_{ij}\right)}$$

(14)

where $u_{ij}(\cdot )$ and $u_{ij}\left({\xi }_{ij}\right)$ are the utility value function and the utility value of alternative $A_i$ in terms of criterion $C_j$, respectively.

In SMAA-2, the criteria value ${\xi }_{ij}$ ($i=1,2,\cdots ,N;j=1,2,\cdots ,n$) is assumed to be a stochastic variable with an assumed or estimated joint probability density function $f\left(\mathsf{\xi }\right)$ in the criteria value space $\mathsf{\Omega }$. Similarly, the preferences of decision-makers are assumed to be unknown or partially known and can be represented by a weight distribution with density function $g\left(\mathsf{\omega }\right)$ in the feasible weight space $W$, defined as follows:

$$W=\left\{\mathsf{\omega }\in R^n:{\displaystyle {\sum }_{j=1}^n{\omega }_j=1};\forall j,1\ge {\omega }_j\ge 0\right\}$$

(15)

When the weights are totally unknown, $g\left(\mathsf{\omega }\right)$ is thought to be a uniform distribution density function and is expressed as follows:

$$g\left(\mathsf{\omega }\right)=1/\mathrm{vol}\left(W\right)$$

(16)

where $\mathrm{vol}\left(W\right)=\sqrt{n}/\left(n-1\right)!$ is the volume of the feasible weight space, an (n − 1)-dimensional simplex.

Based on the above assumptions, stochastic criteria and weight distributions can be mapped into utility distributions $U\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)$ through the utility function. The utility distributions are employed to compute the rank of each alternative, which is defined as an integer from the best rank (=1) to the worst rank (=N) by means of the rank function [34]:

$$rank\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)=1+{\displaystyle \sum _{s=1}^N\rho \left(U\left({\mathsf{\xi }}_s,\mathsf{\omega }\right)>U\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)\right)}$$

(17)

where $\rho \left(ture\right)=1$ and $\rho \left(false\right)=0$. Then, the favorable rank weight space $W_i^r\left(\mathsf{\xi }\right)$ is defined by restricting $rank\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)=r$, which means any weight $\mathsf{\omega }\in W_i^r\left(\mathsf{\xi }\right)$ can make alternative $A_i$ rank r. Thus, the favorable rank weight space $W_i^r\left(\mathsf{\xi }\right)$ can be expressed as follows:

$$W_i^r\left(\mathsf{\xi }\right)=\left\{\mathsf{\omega }\in W:rank\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)=r\right\}$$

(18)

By analyzing $W_i^r\left(\mathsf{\xi }\right)$, the SMAA-2 method obtains four descriptive measures to assist decision-makers in choosing the best alternative. These measures are the rank acceptability index, the central weight vector, the confidence factor, and the holistic acceptability index.

The rank acceptability index $b_i^r$ is defined as the expected volume of the favorable rank weight space $W_i^r\left(\mathsf{\xi }\right)$ and measures the variety of different weights which guarantee a ranking r for alternative $A_i$. This index is a multidimensional integral over criteria distributions and the favorable rank weight space, which can be computed numerically as follows:

$$b_i^r={\displaystyle {\int }_{\mathsf{\Omega }}f\left(\mathsf{\xi }\right)}{\displaystyle {\int }_{W_i^r\left(\mathsf{\xi }\right)}g\left(\mathsf{\omega }\right)}d\mathsf{\omega }d\mathsf{\xi }$$

(19)

Obviously, $b_i^r$ is between 0 and 1; 0 indicates that alternative $A_i$ will never acquire a given rank, and 1 indicates that it will always acquire the given rank r with arbitrary weights. The larger the index, the more the decision-makers are willing to accept the given rank r for alternative $A_i$. The best alternatives are those with high acceptability indices for the best ranks.

The central weight vector ${\mathsf{\omega }}_i^c$ is defined as the expected center of gravity of the favorable weight space $W_i^1\left(\mathsf{\xi }\right)=\left\{\mathsf{\omega }\in W:rank\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)=1\right\}$, and it is computed numerically as an integral of the weight vector over the criteria distributions and the weight distributions:

$${\mathsf{\omega }}_i^c=\frac{1}{b_i^1}{\displaystyle {\int }_{\mathsf{\Omega }}f\left(\mathsf{\xi }\right)}{\displaystyle {\int }_{W_i^1\left(\mathsf{\xi }\right)}g\left(\mathsf{\omega }\right)\mathsf{\omega }}d\mathsf{\omega }d\mathsf{\xi }$$

(20)

The central weight vector ${\mathsf{\omega }}_i^c$ is the best single vector representing the preference of a typical decision-maker who favors alternative $A_i$, and can help decision-makers understand the correspondence between different weight vectors and different decision-maker choices.

The confidence factor $p_i^c$ is the probability for alternative $A_i$ to be the preferred one when the central weight vector ${\mathsf{\omega }}_i^c$ is adopted. It is an integral over criteria distributions:

$$p_i^c={\displaystyle {\int }_{\mathsf{\xi }:rank\left({\mathsf{\xi }}_i,{\mathsf{\omega }}_i^c\right)=1}f\left(\mathsf{\xi }\right)d\mathsf{\xi }}$$

(21)

The confidence factor $p_i^c$ is used to measure whether the criteria data are accurate enough to discern alternative $A_i$ based on the central weight vector ${\mathsf{\omega }}_i^c$. For any given weight vector and alternative, the corresponding confidence factor can be computed in a similar manner. The confidence factor can be described as the proportion of the stochastic criteria space that makes the alternative the best one with the given weight vector.

The holistic acceptability index $e_i$ is a comprehensive index to examine the overall acceptability of alternative $A_i$ by combining all rank acceptability indices as follows:

$$e_i={\displaystyle \sum _{r=1}^m{\alpha }_r{b}_i^r}$$

(22)

where ${\alpha }_r$ is a metaweight. Lahdelma et al. [32] suggested choosing the alternative with the highest holistic acceptability index as the best one. In addition, considering the requirements for metaweights, i.e., they should be nonnegative, normalized, and nonincreasing when the rank increases (${\alpha }_1\ge {\alpha }_2\ge \cdots \ge {\alpha }_N\ge 0$), Lahdelma et al. proposed three kinds of metaweights, namely linear weights ${\alpha }_r=\left(N-r\right)/\left(N-1\right)$, inverse weights ${\alpha }_r=1/r$, and centroid weights ${\alpha }_r={\displaystyle {\sum }_{i=r}^{N}1/i}/{\displaystyle {\sum }_{i=1}^{N}1/i}$ [32]. Figure 3 illustrates the differences in linear, inverse, and centroid weights for 10 ranks. It can be observed that linear weights allocate more weight to the mediocre ranks, while inverse and centroid weights emphasize the best ranks. Compared to centroid weights, inverse weights are more evenly distributed over the worst ranks, which indicates that the aggregation by inverse weights is more insensitive to the order among the worst ranks. Therefore, in this research, we use centroid weights.

The multidimensional integrals in Equations (19)–(21) are almost impossible to compute analytically because the distributions $f\left(\mathsf{\Omega }\right)$ and $g\left(W\right)$ vary according to their application and can be arbitrarily complex. For this reason, Tervonen et al. proposed to employ the Monte Carlo simulation to compute approximate values for high-dimensional integrals, where the required number of iterations is set as 10,000 to achieve accuracy with 95% confidence. More details about Monte Carlo simulation-based SMAA-2 computations can be found in Reference [34].

Since it was first proposed, the SMAA-2 method has been applied in various fields, such as the evaluation of combined heat and power units [35], task assignment for multi-UAV [31], strategic forest planning [36], reservoir flood control operations [33], project portfolio optimization [37], and more. However, this method still suffers from a notable shortcoming: The application of the method is based on one criteria value matrix. When there are multiple criteria value matrices from different ISs, the SMAA-2 method is difficult to apply directly. In view of the above shortcoming, we extend the SMAA-2 method in this research. The details are introduced in the next section.

5. The Extended SMAA-2 Method

As can be seen from Figure 1, the extended SMAA-2 method first aggregates multiple criteria value matrices from different information sources (ISs) under the condition that the criteria weight information is unavailable and the weight information for the ISs is partially available in order to acquire the final criteria value matrix. Then, the holistic acceptability indices are computed through the Monte Carlo simulation based on the matrix. To achieve matrices aggregation under serve uncertainty, we propose a novel aggregation method. To compute the holistic acceptability indices, we integrate the aggregation method into the original SMAA-2 method.

This extension is based on the assumption that when one certain criterion of an alternative has been evaluated multiple times, the criteria values obtained are stable and tend to be some of the most possible values. According to the normal distribution interval number theory proposed by Wang et al. [38,39], we assume that the target value $v_j^{\left(k\right)}$ and the distance $d_{ij}^{\left(k\right)}$ ($i=1,2,\dots ,M$; $j=1,2,\cdots ,N$; $k=1,2,\cdots ,K$) provided by the sensors (or information sources) are represented by normal distribution interval numbers.

5.1. Normal Distribution Interval Number and Weighted Arithmetic Averaging Operator

Definition 1.

Assuming$\tilde{a}=\left[a^{-},a^{+}\right]$is an interval number, if any$x\in \left[a^{-},a^{+}\right]$follows the same normal distribution$N\left(\mu ,{\sigma }^2\right)$, then$\tilde{a}$is called a normal distribution interval number [38], denoted by$\tilde{\beta }=\left(\mu ,\sigma \right)$. If we let$I$be the set of all normal distribution interval numbers, according to the ‘$3\sigma$rule’ of normal distribution (if$x~N\left(\mu ,{\sigma }^2\right)$, then$P\left(\mu -3\sigma \le x\le \mu +3\sigma \right)=0.9974$) [40] and$\mu$and$\sigma$ are determined as follows:

$$\mu -3\sigma =a^{-},\mu +3\sigma =a^{+}$$

(23)

The results are as follows:

$$\mu =\left(a^{-}+a^{+}\right)/2,\sigma =\left(a^{+}-a^{-}\right)/6$$

(24)

Based on Definition 1 and the principle that the linear combination of a finite number of independent normal random variables still follows a normal distribution, we can define the computation rules of normal distribution interval numbers as follows:

Definition 2.

Assuming${\tilde{\beta }}_1=\left({\mu }_1,{\sigma }_1\right)$and${\tilde{\beta }}_2=\left({\mu }_2,{\sigma }_2\right)$are arbitrarily two normal distribution interval numbers, then (1)${\tilde{\beta }}_1+{\tilde{\beta }}_2=\left({\mu }_1+{\mu }_2,\sqrt{{\sigma }_1^2+{\sigma }_2^2}\right)$, (2)$\lambda {\tilde{\beta }}_1=\left(\lambda {\mu }_1,\lambda {\sigma }_1\right)$,$\lambda \ge 0$.

Furthermore, according to the addition commutative law and the multiplication allocating law of rational numbers, it is easy to prove that (1) ${\tilde{\beta }}_1+{\tilde{\beta }}_2={\tilde{\beta }}_2+{\tilde{\beta }}_1$; (2) $\lambda \left({\tilde{\beta }}_1+{\tilde{\beta }}_2\right)=\lambda {\tilde{\beta }}_1+\lambda {\tilde{\beta }}_2$, $\lambda \ge 0$; (3) ${\lambda }_1{\tilde{\beta }}_1+{\lambda }_2{\tilde{\beta }}_1=\left({\lambda }_1+{\lambda }_2\right){\tilde{\beta }}_1$, ${\lambda }_1,{\lambda }_2\ge 0$. Based on the distance measure of interval numbers in References [41,42], we define the distance between two normal distribution interval numbers as follows:

Definition 3.

Assuming${\tilde{\beta }}_1=\left({\mu }_1,{\sigma }_1\right)$and${\tilde{\beta }}_2=\left({\mu }_2,{\sigma }_2\right)$are arbitrarily two normal distribution interval numbers, then the distance between them can be computed as follows:

$$dis\left({\tilde{\beta }}_1,{\tilde{\beta }}_2\right)=\sqrt{{\left({\mu }_1-{\mu }_2\right)}^2+3\left({\sigma }_1{}^2+{\sigma }_2{}^2\right)-{\epsilon }^2/6}$$

(25)

$$\epsilon =\{\begin{array}{cc}minup-maxlow& minup\ge maxlow\\ 0& otherwise\end{array}$$

(26)

where$minup=\min \left({\mu }_1+3{\sigma }_1,{\mu }_2+3{\sigma }_2\right)$,$maxlow=\max \left({\mu }_1-3{\sigma }_1,{\mu }_2-3{\sigma }_2\right)$.

Definition 4.

${\tilde{\beta }}_i=\left({\mu }_i,{\sigma }_i\right)$$\left(i=1,2,\cdots ,m\right)$is a group of normal distribution interval numbers, and assuming WAA:$I^n\to I$

$${\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2,\cdots ,{\tilde{\beta }}_m\right)={\omega }_1{\tilde{\beta }}_1+{\omega }_2{\tilde{\beta }}_2+\cdots +{\omega }_m{\tilde{\beta }}_m$$

(27)

is always satisfied for arbitrary$\mathsf{\omega }={\left[{\omega }_1,{\omega }_2,\cdots ,{\omega }_m\right]}^{\mathrm{T}}$(${\omega }_i\in [0,1]$,${\displaystyle {\sum }_{i=1}^m{\omega }_i}=1$), then WAA is called the weighted arithmetic averaging operator for normal distribution interval numbers, and$\mathsf{\omega }$is called the weighted vector of${\tilde{\beta }}_i=\left({\mu }_i,{\sigma }_i\right)$.

Theorem 1.

If${\tilde{\beta }}_i=\left({\mu }_i,{\sigma }_i\right)$$\left(i=1,2,\cdots ,m\right)$is a group of normal distribution interval numbers, and$\mathsf{\omega }={\left[{\omega }_1,{\omega }_2,\cdots ,{\omega }_m\right]}^{\mathrm{T}}$(${\omega }_i\in [0,1]$,${\displaystyle {\sum }_{i=1}^m{\omega }_i}=1$) is the corresponding weighted vector, then${\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2,\cdots ,{\tilde{\beta }}_m\right)$is a normal distribution interval number, satisfying:

$${\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2,\cdots ,{\tilde{\beta }}_m\right)=\left({\displaystyle \sum _{i=1}^m{\omega }_i{\mu }_i},\sqrt{{\displaystyle \sum _{i=1}^m{\omega }_i^2{\sigma }_i^2}}\right)$$

(28)

Proof. 

Here, we use the mathematical induction to prove Theorem 1.

(1) When $n=2$, ${\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2\right)={\omega }_1{\tilde{\beta }}_1+{\omega }_2{\tilde{\beta }}_2$. Then, according to Definition 2, ${\omega }_1{\tilde{\beta }}_1=\left({\omega }_1{\mu }_1,{\omega }_1{\sigma }_1\right)$ and ${\omega }_2{\tilde{\beta }}_2=\left({\omega }_2{\mu }_2,{\omega }_2{\sigma }_2\right)$ are two normal distribution interval numbers. Therefore, we have:

$${\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2\right)={\omega }_1{\tilde{\beta }}_1+{\omega }_2{\tilde{\beta }}_2=\left({\omega }_1{\mu }_1,{\omega }_1{\sigma }_1\right)+\left({\omega }_2{\mu }_2,{\omega }_2{\sigma }_2\right)=\left({\omega }_1{\mu }_1+{\omega }_2{\mu }_2,\sqrt{{\omega }_1^2{\sigma }_1^2+{\omega }_2^2{\sigma }_2^2}\right)$$

(29)

(2) Assuming that when $m=k$, Theorem 1 holds:

$${\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2,\cdots ,{\tilde{\beta }}_k\right)=\left({\displaystyle \sum _{i=1}^k{\omega }_i{\mu }_i},\sqrt{{\displaystyle \sum _{i=1}^k{\omega }_i^2{\sigma }_i^2}}\right)$$

(30)

(3) When $m=k+1$, according to Definition 2:

$$\begin{array}{l}{\mathrm{WAA}}_{\mathsf{\omega }}\left({\tilde{\beta }}_1,{\tilde{\beta }}_2,\cdots ,{\tilde{\beta }}_k,{\tilde{\beta }}_{k+1}\right)=\left({\displaystyle \sum _{i=1}^k{\omega }_i{\mu }_i},\sqrt{{\displaystyle \sum _{i=1}^k{\omega }_i^2{\sigma }_i^2}}\right)+{\omega }_{k+1}\left({\mu }_{k+1},{\sigma }_{k+1}\right)\\ =\left({\displaystyle \sum _{i=1}^k{\omega }_i{\mu }_i}+{\omega }_{k+1}{\mu }_{k+1},\sqrt{{\left(\sqrt{{\displaystyle \sum _{i=1}^k{\omega }_i^2{\sigma }_i^2}}\right)}^2+{\left({\omega }_{k+1}{\sigma }_{k+1}\right)}^2}\right)\\ =\left({\displaystyle \sum _{i=1}^{k+1}{\omega }_i{\mu }_i},\sqrt{{\displaystyle \sum _{i=1}^{k+1}{\omega }_i^2{\sigma }_i^2}}\right)\end{array}$$

(31)

That means Theorem 1 holds when $m=k+1$.

Therefore, it can be concluded from Equations (1) and (2) that Theorem 1 holds. □

5.2. Aggregation Method of Criterion Value Matrices from Different Information Sources

5.2.1. WAA Operator-Based Criteria Value Matrices Aggregation

Definition 5.

If${\tilde{\beta }}_{ij}=\left({\mu }_{ij},{\sigma }_{ij}\right)$$\left(i=1,2,\cdots ,N;j=1,2,\cdots ,n\right)$are a group of normal distribution interval numbers, then the matrix$S={\left[{\tilde{\beta }}_{ij}\right]}_{N\times n}$is called a normal distribution interval number matrix.

Here, we consider a multicriteria decision-making problem, including a set of N alternatives$Set\_alt=\left\{A_1,A_2,\cdots ,A_N\right\}$from which one is to be chosen. All alternatives are evaluated in terms of n criteria$Set\_cri=\left\{C_1,C_2,\cdots ,C_n\right\}$by a set of K experts (or information sources)$Set\_is=\left\{I{S}_1,I{S}_2,\cdots ,I{S}_K\right\}$. Thus, there are K criteria value matrices$S^{\left(1\right)},S^{\left(2\right)},\cdots ,S^{\left(K\right)}$, each of which has N rows and n columns, corresponding to N alternatives and n criteria, respectively. According to the above assumption, i.e., the criteria values are represented by normal distribution interval numbers, the kth criteria value matrix$S^{\left(k\right)}={\left[{\tilde{\beta }}_{ij}^{\left(k\right)}\right]}_{N\times n}$$\left(i=1,2,\cdots ,N;j=1,2,\cdots ,n;k=1,2,\cdots ,K\right)$is a normal distribution interval number matrix where${\tilde{\beta }}_{ij}^{\left(k\right)}=\left({\mu }_{ij}^{\left(k\right)},{\sigma }_{ij}^{\left(k\right)}\right)$is a normal distribution interval number. In order to facilitate criteria value matrices aggregation, all matrices need to be normalized by Reference [43]:

$${\overline{\mu }}_{ij}^{\left(k\right)}=\{\begin{array}{cc}{\mu }_{ij}^{\left(k\right)}/\underset{i}{\max }\left({\mu }_{ij}^{\left(k\right)}\right)& C_j\mathrm{is}\mathrm{a}\mathrm{benefit}-\mathrm{type}\mathrm{criterion}\\ \underset{i}{\min }\left({\mu }_{ij}^{\left(k\right)}\right)/{\mu }_{ij}^{\left(k\right)}& C_j\mathrm{is}\mathrm{a}\cos \mathrm{t}-\mathrm{type}\mathrm{criterion}\end{array},$$

(32)

$${\overline{\sigma }}_{ij}^{\left(k\right)}={\sigma }_{ij}^{\left(k\right)}/\underset{i}{\max }\left({\mu }_{ij}^{\left(k\right)}\right).$$

(33)

In the following, unless otherwise stated, all criteria value matrices are assumed to be normalized.

Definition 6.

$S^{\left(0\right)}={\left[{\tilde{\beta }}_{ij}^{\left(0\right)}\right]}_{N\times n}$$\left(i=1,2,\cdots ,N;j=1,2,\cdots ,n\right)$is called a final criteria value matrix only when${\tilde{\beta }}_{ij}^{\left(0\right)}=\left({\mu }_{ij}^{\left(0\right)},{\sigma }_{ij}^{\left(0\right)}\right)$is a normal distribution interval number, and${\mu }_{ij}^{\left(0\right)}$and${\sigma }_{ij}^{\left(0\right)}$are computed through the WAA operator as follows:

$${\mu }_{ij}^{\left(0\right)}={\displaystyle \sum _{k=1}^t{\lambda }_k{\mu }_{ij}^{\left(k\right)}},$$

(34)

$${\sigma }_{ij}^{\left(0\right)}=\sqrt{{\displaystyle \sum _{k=1}^t{\left({\lambda }_k{\sigma }_{ij}^{\left(k\right)}\right)}^2}},$$

(35)

where $\mathsf{\lambda }=\left[{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_K\right]$ is the integrated weight vector of K information sources (ISs), satisfying ${\lambda }_k\in [0,1]$ and ${\displaystyle {\sum }_{k=1}^K{\lambda }_k}=1$. The integrated weight of the kth IS ${\lambda }_k$ is computed as follows:

$${\lambda }_k=\frac{{\eta }_k{}^{\tau }{\rho }_k{}^{1-\tau }}{{\displaystyle \sum _{k=1}^K{\eta }_k{}^{\tau }{\rho }_k{}^{1-\tau }}},$$

(36)

where $\tau \in [0,1]$ is an adjustment factor and ${\eta }_k\in [0,1]$ and ${\rho }_k\in [0,1]$ are the subjective weight and the objective weight of the kth IS, respectively, satisfying ${\displaystyle {\sum }_{k=1}^K{\eta }_k}=1$ and ${\displaystyle {\sum }_{k=1}^K{\rho }_k}=1$. The subjective weights are usually predefined and would not be adjusted once determined. The objective weight of each IS is determined by the accuracy and credibility of the criteria values it provides [44]. The iterative algorithm for computing objective weights is introduced in the following subsection.

5.2.2. Iterative Algorithm for Computing Objective Weights

As objective weights are related to the accuracy and credibility of criteria values, measuring their accuracy and credibility is critical. In this research, we use the deviation between the criteria value matrix of one IS and the final criteria value matrix as the measurement [44] and then proposed an iterative algorithm to compute the objective weights of the ISs.

Definition 7.

If we let$S^{\left(k\right)}={\left[{\tilde{\beta }}_{ij}^{\left(k\right)}\right]}_{N\times n}$and$S^{\left(0\right)}={\left[{\tilde{\beta }}_{ij}^{\left(0\right)}\right]}_{N\times n}$$\left(i=1,2,\cdots ,N;j=1,2,\cdots ,n;k=1,2,\cdots ,K\right)$be the criteria value matrix of the kth IS and the final criteria value matrix, respectively, then the deviation between these two matrices can be computed:

$$Dev\left(S^{\left(k\right)},S^{\left(0\right)}\right)={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^n{\omega }_j\cdot dis\left({\tilde{\beta }}_{ij}^{\left(k\right)},{\tilde{\beta }}_{ij}^{\left(0\right)}\right)}},$$

(37)

where ${\omega }_j$ $\left(j=1,2,\cdots ,n\right)$ is the weight of the jth criterion, satisfying ${\omega }_j\in \left[0,1\right]$ and ${\displaystyle {\sum }_{j=1}^n{\omega }_j}=1$, and $dis\left({\tilde{\beta }}_{ij}^{\left(k\right)},{\tilde{\beta }}_{ij}^{\left(0\right)}\right)$ is the distance between ${\tilde{\beta }}_{ij}^{\left(k\right)}$ and ${\tilde{\beta }}_{ij}^{\left(0\right)}$, as defined by Definition 3.

As can be seen from Definitions 7 and 3, a smaller deviation $Dev\left(S^{\left(k\right)},S^{\left(0\right)}\right)$ represents a higher similarity between $S^{\left(k\right)}$ and $S^{\left(0\right)}$, which means that the kth IS should be given a higher objective weight. Based on this perception, the objective weight ${\rho }_k$ is then computed:

$${\rho }_k=\frac{{\displaystyle \sum _{k=1}^{K}Dev\left(S^{\left(k\right)},S^{\left(0\right)}\right)}-Dev\left(S^{\left(k\right)},S^{\left(0\right)}\right)}{{\displaystyle \sum _{k=1}^{K}Dev\left(S^{\left(k\right)},S^{\left(0\right)}\right)}}.$$

(38)

The iteration end condition is set as $\Delta \rho \le \zeta$, where $\zeta$ is the error margin, a very small number (i.e., 10⁻⁸), and $\Delta \rho$ is computed as follows:

$$\Delta \rho =\sqrt{{\displaystyle \sum _{k=1}^K\left({\rho }_k^{\left(l\right)}-{\rho }_k^{\left(l-1\right)}\right)}}.$$

(39)

In Equation (39), ${\rho }_k^{\left(l\right)}$ and ${\rho }_k^{\left(l-1\right)}$ are the objective weight of the kth IS at the lth and the (l − 1)th iteration, respectively. For convenience, the initial objective weight vector is set as ${\rho }^{\left(1\right)}=\left[1/K,1/K,\cdots ,1/K\right]$. Once the objective weight vector is determined, then the integrated weight vector can be computed by Equation (36). On this basis, all criteria value matrices can be aggregated to obtain the final criteria value matrix by Equations (34) and (35). The pseudo-code of the aggregation method based on the iterative algorithm is presented in Algorithm 1. On the basis of the output of Algorithm 1 (i.e., the final criteria value matrix $S^{\left(0\right)}$), the SMAA-2 method can be subsequently applied.

Algorithm 1: Pseudo-code of the aggregation method
Input:	number of alternatives, criteria, and information sources: $N$, $n$, $K$;
	criteria value matrices: $S^{\left(1\right)},S^{\left(2\right)},\cdots ,S^{\left(t\right)}$;
	subjective weight vector $\mathsf{\eta }=\left[{\eta }_1,{\eta }_2,\cdots ,{\eta }_K\right]$ and initial objective weight vector, ${\mathsf{\rho }}^{\left(1\right)}=\left[1/K,1/K,\cdots ,1/K\right]$;
	criteria weight vector $\mathsf{\omega }=\left[{\omega }_1,{\omega }_2,\cdots ,{\omega }_n\right]$;
	error margin used in computing the objective weight vector: $\zeta$.
Output:	final criteria value matrix $S^{\left(0\right)}$.
Procedures:	$l\leftarrow 1$
	$\Delta {\rho }^{\left(l\right)}\leftarrow \infty$
	while$\Delta {\rho }^{\left(l\right)}\ge \zeta$do
	compute $\mathsf{\lambda }$ by (33)
	$l\leftarrow l+1$
	for $k\leftarrow 1$ to $K$ do
	compute $Dev\left(S^{\left(k\right)},S^{\left(0\right)}\right)$ by (22), (23), (31), (32) and (34)
	end for
	update ${\mathsf{\rho }}^{\left(l\right)}$ by (35)
	update $\Delta {\rho }^{\left(l\right)}$ by (36)
	end while
	return$S^{\left(0\right)}$

5.3. Integrating the Aggregation Method into the Original SMAA-2 Method

By integrating the aggregation method into the original SMAA-2 method, the stochastic multicriteria acceptability analysis can be performed to solve the task assignment problem formulated in Section 3. The implementation of the extended SMAA-2 method is also mainly based on the Monte Carlo simulation with 10,000 iterations to achieve accuracy with 95% confidence and consists of two phases of which the first is the computation of the rank acceptability index $b_i^r$, the central weight vector ${\mathsf{\omega }}_i^c$, and the holistic acceptability index $e_i$, and the second is the computation of the confidence factor $p_i^c$. The pseudo-code is specified in Algorithm 2.

Algorithm 2: Pseudo-code of the extended SMAA-2 method
Symbols:	$h_i^j$, number of times alternative $A_i$ is evaluated into rank j in the Monte Carlo simulation;
	$N_{ite}$, number of iterations, set as 10000;
	$r=\left[r_1,r_2,\cdots ,r_N\right]$, vector of ranks of alternatives;
	$U=\left[U_1,U_2,\cdots ,U_N\right]$, utility value vector of alternatives;
	$Agg\left(\mathsf{\omega }\right)$, function returning the final criteria value matrix $S^{\left(0\right)}$ based on the criteria weight vector $\mathsf{\omega }$ (details are presented in Algorithm 1)
	$RAN{D}_{\mathsf{\Omega }}\left(S\right)$, function returning a random criteria value matrix from criteria distribution $f\left(\mathsf{\Omega }\right)$;
	$RAN{D}_W$, function returning a random weight vector from weight distribution $g\left(W\right)$;
	$RANK\left(U\right)$, function returning a vector of ranks corresponding to $U$.
Input:	number of alternatives, criteria, and decision-makers: $N$, $n$, $K$;
	criteria value matrices: $S^{\left(1\right)},S^{\left(2\right)},\cdots ,S^{\left(K\right)}$;
	subjective weight vector and initial objective weight vector $\mathsf{\eta }=\left[{\eta }_1,{\eta }_2,\cdots ,{\eta }_K\right]$, ${\mathsf{\rho }}^{\left(1\right)}=\left[1/K,1/K,\cdots ,1/K\right]$;
	error margin used in computing the objective weight vector: $\zeta$;
	number of iterations of the Monte Carlo simulation: $N_{ite}$.
Output:	rank acceptability index $b_i^r$;
	central weight vector ${\mathsf{\omega }}_i^c$;
	confidence factor $p_i^c$;
	holistic acceptability index $e_i$.
Procedures:	// Initialize ${\mathsf{\omega }}_i^c$, $h_i^j$ and $e_i$ for $i\leftarrow 1$ to $N$ do ${\mathsf{\omega }}_i^c\leftarrow 0$ $e_i\leftarrow 0$ for $j\leftarrow 1$ to $N$ do $h_i^j\leftarrow 0$ end for end for // Phase 1, compute $b_i^r$, ${\mathsf{\omega }}_i^c$ and $e_i$ // Main loop for $k\leftarrow 1$ to $N_{ite}$ do  $\mathsf{\omega }\leftarrow RAN{D}_W$ //produce $S^{\left(0\right)}$ by $Agg(\cdot )$  $S^{\left(0\right)}\leftarrow Agg\left(\mathsf{\omega }\right)$ //$S$ is the generated random criteria value matrix $S\leftarrow RAN{D}_{\mathsf{\Omega }}\left(S^{\left(0\right)}\right)$  for $i\leftarrow 1$ to $N$ do  $U_i\leftarrow U\left({\mathsf{\xi }}_i,\mathsf{\omega }\right)$ end for $r\leftarrow RANK\left(U\right)$ for $i\leftarrow 1$ to $N$ do  $h_i^{r_i}\leftarrow h_i^{r_i}+1$  if $r_i=1$ then   ${\mathsf{\omega }}_i^c\leftarrow {\mathsf{\omega }}_i^c+\mathsf{\omega }$  end if end for end for // Compute ${\mathsf{\omega }}_i^c$, $b_i^r$ and $e_i$ for $i\leftarrow 1$ to $N$ do  if $h_i^1>0$ then ${\mathsf{\omega }}_i^c\leftarrow {\mathsf{\omega }}_i^c/h_i^1$ end if for $j\leftarrow 1$ to $N$ do $b_i^j\leftarrow h_i^j/N_{ite}$ $e_i\leftarrow e_i+{\alpha }_j{b}_i^j$ end for end for	// Phase 2, compute $p_i^c$ // Initialize $p_i^c$, and generate the final criteria value matrix $S_i^{\left(0\right)}$ corresponding to ${\mathsf{\omega }}_i^c$ for $i\leftarrow 1$ to $N$ do $p_i^c\leftarrow 0$ $S_i^{\left(0\right)}\leftarrow Agg\left({\mathsf{\omega }}_i^c\right)$ end for // Main loop for $j\leftarrow 1$ to $N_{ite}$ do  for $i\leftarrow 1$ to $N$ do $S\leftarrow RAN{D}_{\mathsf{\Omega }}\left(S_i^{\left(0\right)}\right)$ $U_i=U\left({\mathsf{\xi }}_i,{\mathsf{\omega }}_i^c\right)$ for $k\leftarrow 1$ to $N$ and $k\ne i$ do  if $U\left({\mathsf{\xi }}_k,{\mathsf{\omega }}_i^c\right)>U_i$ then   go to worse  end if end for $p_i^c\leftarrow p_i^c+1$ worse: end for end for for $i\leftarrow 1$ to $N$ do  $p_i^c\leftarrow p_i^c/N_{ite}$ end for

As can be seen from Algorithm 2, for each UCAV, we need to compute the holistic acceptability indices in terms of each target (here, a target is equivalent to an alternative) by the extended SMAA-2 method. In the next section, we convert the task assignment model established in Section 2 into an integer linear programming model. Thus, the holistic acceptability indices can be applied to task assignment.

6. Converting the Task Assignment Model Based on Holistic Acceptability Indices

The task assignment model established in Section 3 involves M UCAVs and N targets. By applying the extended SMAA-2 method M times, we obtain a holistic acceptability index matrix $E={\left[e_{ij}\right]}_{M\times N}$, where $e_{ij}$ ($i=1,2,\cdots ,M;j=1,2,\cdots ,N$) is the holistic acceptability index for assigning the jth target $T_j$ to the ith UCAV $UCA{V}_i$. Then, we can convert the established task assignment model into an integer linear programming model in which the objective function aims to maximize the holistic acceptability index of the task assignment scheme. The converted model is expressed as follows:

$$\max G\left(x\right)={\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N{x}_{ij}{e}_{ij}}}$$

(40)

subject to the following:

$${\displaystyle \sum _{i=1}^M{x}_{ij}}=1$$

(41)

$${\displaystyle \sum _{j=1}^N{x}_{ij}}=1$$

(42)

$$x_{ij}\in \left\{0,1\right\},i=1,2,\cdots ,M;j=1,2,\cdots ,N$$

(43)

The reasons the previous model can be converted in this way are explained as follows: The value of the objective function $F\left(x\right)$ ($F\left(x\right)={\omega }_B{O}_B+{\omega }_C{O}_C$) in the previous model is an absolute value of combat reward. However, due to uncertain parameters, the value is not able to be computed directly. According to the SMAA-2 method, Lahdelma et al. suggested choosing the alternative with the highest holistic acceptability index [32]. In other words, $e_{ij}$ in the task assignment model, which describes the overall acceptability for assigning $T_j$ to $UCA{V}_i$, can be interpreted as a score of the assignment. For $UCA{V}_i$, the larger the score, the more $UCA{V}_i$ should attack $T_j$ to achieve a larger combat reward given the criteria distributions $f\left(\mathsf{\Omega }\right)$ and the criteria weight distributions $g\left(W\right)$. Thus, the holistic acceptability index $e_{ij}$ can be treated as a relative value of combat reward. Therefore, it is reasonable to convert the objective function $F\left(x\right)$ in the previous model to a linear combination of $e_{ij}$, expressed as Equation (40).

In terms of constraints, as the constraints expressed in Equations (7)–(12) on the criteria weights and the aggregation method of the criterion value matrices in the previous model have been satisfied in computing the holistic acceptability indices, the converted model only contains the three constraints represented by Equations (41)–(43) that specify the range of the decision variable $x_{ij}$. Thus, there are no uncertain parameters in the converted model, which indicates an easier process of searching for solutions.

7. Case Study

In this section, we examine the proposed method through a case study. Here, we consider a combat scenario where four UCAVs under the command of one command center are to attack four geographically dispersed ground targets. The initial battlefield situation is shown in Figure 4. Before the mission begins, the command center has received target information that is independently acquired by four different sensors through multiple reconnaissances, as shown in

Table 1. Target information independently acquired by sensors 1, 2, 3, and 4.

Sensors	Target Value $({\mathit{v}}_{\mathit{j}}=[{\mathit{v}}_{\mathit{j}}^{-},{\mathit{v}}_{\mathit{j}}^{+}])$	$\mathbf{Flight}\mathbf{Distance}\mathbf{between}\mathit{U}\mathit{C}\mathit{A}{\mathit{V}}_{\mathit{i}}$$\mathbf{and}{\mathit{T}}_{\mathit{j}}$$({\mathit{d}}_{\mathit{i}\mathit{j}}=[{\mathit{d}}_{\mathit{i}\mathit{j}}^{-},{\mathit{d}}_{\mathit{i}\mathit{j}}^{+}]$$,\mathit{i}=1,2,3,4;\mathit{j}=1,2,3,4)$
Sensor 1	$\begin{array}{c}\left[7.0,9.0\right]\\ \left[6.0,8.0\right]\\ \left[3.0,7.0\right]\\ \left[2.0,5.0\right]\end{array}$	$\begin{array}{cccc}\left[725,762\right]& \left[721,755\right]& \left[595,637\right]& \left[612,651\right]\\ \left[751,792\right]& \left[720,759\right]& \left[672,715\right]& \left[661,704\right]\\ \left[642,684\right]& \left[567,609\right]& \left[666,703\right]& \left[608,647\right]\\ \left[630,672\right]& \left[578,620\right]& \left[602,642\right]& \left[565,607\right]\end{array}$
Sensor 2	$\begin{array}{c}\left[6.5,8.5\right]\\ \left[4.0,7.5\right]\\ \left[5.0,8.0\right]\\ \left[4.0,7.0\right]\end{array}$	$\begin{array}{cccc}\left[708,745\right]& \left[738,771\right]& \left[614,655\right]& \left[602,641\right]\\ \left[731,772\right]& \left[732,771\right]& \left[690,732\right]& \left[650,692\right]\\ \left[619,660\right]& \left[570,613\right]& \left[679,716\right]& \left[597,636\right]\\ \left[608,650\right]& \left[587,628\right]& \left[617,658\right]& \left[554,596\right]\end{array}$
Sensor 3	$\begin{array}{c}\left[5.0,9.0\right]\\ \left[5.0,8.5\right]\\ \left[4.0,9.0\right]\\ \left[5.0,8.0\right]\end{array}$	$\begin{array}{cccc}\left[715,752\right]& \left[709,742\right]& \left[597,638\right]& \left[611,650\right]\\ \left[739,780\right]& \left[703,742\right]& \left[671,713\right]& \left[661,703\right]\\ \left[627,669\right]& \left[545,588\right]& \left[658,695\right]& \left[609,648\right]\\ \left[616,658\right]& \left[559,601\right]& \left[596,637\right]& \left[565,607\right]\end{array}$
Sensor 4	$\begin{array}{c}\left[6.0,8.0\right]\\ \left[4.0,6.5\right]\\ \left[3.0,8.5\right]\\ \left[4.5,7.5\right]\end{array}$	$\begin{array}{cccc}\left[728,765\right]& \left[733,767\right]& \left[610,651\right]& \left[621,660\right]\\ \left[749,790\right]& \left[733,772\right]& \left[683,726\right]& \left[665,707\right]\\ \left[630,671\right]& \left[581,623\right]& \left[667,705\right]& \left[602,642\right]\\ \left[623,665\right]& \left[592,634\right]& \left[607,648\right]& \left[564,606\right]\end{array}$

. Other necessary parameters are shown in

Table 2. Scenario parameters related to unmanned combat aerial vehicles (UCAVs).

UCAVs	$\mathbf{No}.\mathbf{of}\mathbf{Missiles}\mathbf{on}\mathit{U}\mathit{C}\mathit{A}{\mathit{V}}_{\mathit{i}}$$\left({\mathit{q}}_{\mathit{i}}\right)$	$\mathbf{Flight}\mathbf{Cost}\mathbf{Per}\mathbf{Unit}\mathbf{Distance}\mathbf{of}\mathit{U}\mathit{C}\mathit{A}{\mathit{V}}_{\mathit{i}}$$\left({\mathit{c}}_{\mathit{i}\mathit{j}}^{\mathit{f}\mathit{l}\mathit{y}}\right)$	$\mathbf{Lethality}\mathbf{Probability}\mathbf{of}\mathbf{Missiles}\mathbf{on}\mathit{U}\mathit{C}\mathit{A}{\mathit{V}}_{\mathit{i}}$$\mathbf{against}{\mathit{T}}_{\mathit{j}}$$\left({\mathit{p}}_{\mathit{i}\mathit{j}}\right)$
$UCA{V}_1$	6	$\left[\begin{array}{cccc}5& 10& 5& 6\\ 4& 6& 6& 8\\ 4& 7& 10& 8\\ 9& 9& 8& 9\end{array}\right]$	$\left[\begin{array}{cccc}0.82& 0.63& 0.86& 0.68\\ 0.91& 0.47& 0.49& 0.72\\ 0.46& 0.78& 0.59& 0.55\\ 0.58& 0.57& 0.71& 0.49\end{array}\right]$
$UCA{V}_2$	8
$UCA{V}_3$	6
$UCA{V}_4$	6

and

Table 3. Parameters related to the aggregation method.

Parameters	Subjective Weight Vector of Sensors $\left(\mathsf{\eta }\right)$	Initial Objective Weight Vector of Sensors $({\mathsf{\rho }}^{\left(1\right)})$	Error Margin $\left(\mathit{\zeta }\right)$	Adjustment Factor $\left(\mathit{\tau }\right)$
Value	[0.20, 0.10, 0.35, 0.35]	[0.25, 0.25, 0.25, 0.25]	10−8	0.5

.

Based on the normal distribution interval number theory, we first present the target information in Table 1 as normal distribution interval numbers by Equation (24), as shown in

Table 4. Target information independently acquired by sensors 1, 2, 3, and 4.

Sensors	Target Value $({\mathit{v}}_{\mathit{j}}=\left({\mathit{\mu }}_{\mathit{j}},{\mathit{\sigma }}_{\mathit{j}}\right))$	$\mathbf{Flight}\mathbf{Distance}\mathbf{between}\mathit{U}\mathit{C}\mathit{A}{\mathit{V}}_{\mathit{i}}$$\mathbf{and}{\mathit{T}}_{\mathit{j}}$$({\mathit{d}}_{\mathit{i}\mathit{j}}=\left({\mathit{\mu }}_{\mathit{i}\mathit{j}},{\mathit{\sigma }}_{\mathit{i}\mathit{j}}\right)$$,\mathit{i}=1,2,3,4;\mathit{j}=1,2,3,4)$
Sensor 1	$\begin{array}{c}\left(8.00,0.33\right)\\ \left(7.00,0.33\right)\\ \left(5.00,0.67\right)\\ \left(3.50,0.50\right)\end{array}$	$\begin{array}{cccc}\left(743.4,6.26\right)& \left(738.3,5.69\right)& \left(615.9,6.90\right)& \left(631.6,6.57\right)\\ \left(771.8,6.87\right)& \left(739.3,6.56\right)& \left(693.5,7.07\right)& \left(682.5,7.03\right)\\ \left(662.9,6.94\right)& \left(588.0,7.06\right)& \left(684.2,6.14\right)& \left(627.7,6.53\right)\\ \left(651.2,7.06\right)& \left(599.1,6.93\right)& \left(621.8,6.75\right)& \left(586.0,7.00\right)\end{array}$
Sensor 2	$\begin{array}{c}\left(7.50,0.33\right)\\ \left(5.75,0.53\right)\\ \left(6.50,0.50\right)\\ \left(5.50,0.50\right)\end{array}$	$\begin{array}{cccc}\left(726.2,6.18\right)& \left(754.5,5.61\right)& \left(634.2,6.70\right)& \left(621.7,6.55\right)\\ \left(751.2,6.84\right)& \left(751.4,6.49\right)& \left(711.0,7.07\right)& \left(671.4,7.03\right)\\ \left(639.6,6.94\right)& \left(591.3,7.07\right)& \left(697.2,6.20\right)& \left(616.8,6.52\right)\\ \left(628.8,7.06\right)& \left(607.8,6.88\right)& \left(637.3,6.79\right)& \left(574.6,7.00\right)\end{array}$
Sensor 3	$\begin{array}{c}\left(7.00,0.67\right)\\ \left(6.75,0.58\right)\\ \left(6.50,0.83\right)\\ \left(6.50,0.50\right)\end{array}$	$\begin{array}{cccc}\left(733.6,6.20\right)& \left(725.7,5.58\right)& \left(617.4,6.87\right)& \left(630.6,6.58\right)\\ \left(759.4,6.85\right)& \left(722.4,6.51\right)& \left(691.7,7.07\right)& \left(682.2,7.04\right)\\ \left(647.9,6.95\right)& \left(566.4,7.06\right)& \left(676.7,6.19\right)& \left(628.8,6.52\right)\\ \left(637.3,7.06\right)& \left(579.8,6.90\right)& \left(616.8,6.79\right)& \left(586.4,6.99\right)\end{array}$
Sensor 4	$\begin{array}{c}\left(7.00,0.33\right)\\ \left(5.25,0.42\right)\\ \left(5.75,0.92\right)\\ \left(6.00,0.50\right)\end{array}$	$\begin{array}{cccc}\left(746.8,6.15\right)& \left(749.8,5.74\right)& \left(630.8,6.86\right)& \left(640.5,6.49\right)\\ \left(769.3,6.81\right)& \left(752.4,6.58\right)& \left(704.4,7.07\right)& \left(686.2,7.01\right)\\ \left(650.5,6.98\right)& \left(602.1,7.06\right)& \left(685.9,6.24\right)& \left(621.7,6.61\right)\\ \left(644.1,7.05\right)& \left(613.0,6.93\right)& \left(627.9,6.82\right)& \left(584.9,7.02\right)\end{array}$

. Then we compute the corresponding criteria value matrix $S^{\left(i,k\right)}$ for each UCAV using Equations (10) and (11). Finally, we normalize $S^{\left(i,k\right)}$ using Equations (32) and (33) and obtain the normalized matrix ${\overline{S}}^{\left(i,k\right)}$, as shown in

Table 5. Normalized criteria value matrices.

${\overline{\mathbf{S}}}^{\left(\mathit{i},\mathit{k}\right)}$	Sensor 1	Sensor 2	Sensor 3	Sensor 4
UCAV1	$\begin{array}{ll}(1.000,0.042)\hfill & (0.829,0.004)\hfill \\ (0.873,0.042)\hfill & (0.417,0.008)\hfill \\ (0.625,0.083)\hfill & (1.000,0.005)\hfill \\ (0.637,0.062)\hfill & (0.813,0.005)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.044)\hfill & (0.873,0.004)\hfill \\ (0.765,0.078)\hfill & (0.420,0.007)\hfill \\ (0.867,0.067)\hfill & (1.000,0.005)\hfill \\ (0.733,0.067)\hfill & (0.850,0.005)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.095)\hfill & (0.842,0.004)\hfill \\ (0.962,0.083)\hfill & (0.425,0.008)\hfill \\ (0.929,0.119)\hfill & (1.000,0.005)\hfill \\ (0.928,0.071)\hfill & (0.816,0.005)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.048)\hfill & (0.845,0.004)\hfill \\ (0.748,0.059)\hfill & (0.421,0.008)\hfill \\ (0.822,0.131)\hfill & (1.000,0.005)\hfill \\ (0.856,0.071)\hfill & (0.821,0.005)\hfill \end{array}$
UCAV2	$\begin{array}{ll}(1.000,0.042)\hfill & (1.000,0.005)\hfill \\ (0.870,0.041)\hfill & (0.696,0.007)\hfill \\ (0.622,0.083)\hfill & (0.742,0.008)\hfill \\ (0.438,0.063)\hfill & (0.565,0.010)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.044)\hfill & (1.000,0.005)\hfill \\ (0.762,0.077)\hfill & (0.667,0.007)\hfill \\ (0.863,0.066)\hfill & (0.704,0.008)\hfill \\ (0.733,0.067)\hfill & (0.559,0.011)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.095)\hfill & (1.000,0.005)\hfill \\ (0.958,0.083)\hfill & (0.701,0.007)\hfill \\ (0.924,0.119)\hfill & (0.732,0.008)\hfill \\ (0.929,0.071)\hfill & (0.557,0.010)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.048)\hfill & (1.000,0.005)\hfill \\ (0.745,0.059)\hfill & (0.682,0.007)\hfill \\ (0.818,0.130)\hfill & (0.728,0.008)\hfill \\ (0.857,0.071)\hfill & (0.561,0.010)\hfill \end{array}$
UCAV3	$\begin{array}{ll}(1.000,0.042)\hfill & (1.000,0.004)\hfill \\ (0.897,0.043)\hfill & (0.644,0.007)\hfill \\ (0.638,0.085)\hfill & (0.388,0.009)\hfill \\ (0.445,0.094)\hfill & (0.528,0.008)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.044)\hfill & (1.000,0.004)\hfill \\ (0.786,0.080)\hfill & (0.618,0.007)\hfill \\ (0.885,0.068)\hfill & (0.367,0.009)\hfill \\ (0.746,0.068)\hfill & (1.000,0.008)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.095)\hfill & (1.000,0.004)\hfill \\ (0.989,0.085)\hfill & (0.654,0.007)\hfill \\ (0.948,0.122)\hfill & (0.383,0.009)\hfill \\ (0.944,0.073)\hfill & (1.000,0.008)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.048)\hfill & (1.000,0.004)\hfill \\ (0.769,0.061)\hfill & (0.617,0.007)\hfill \\ (0.838,0.134)\hfill & (0.379,0.009)\hfill \\ (0.872,0.073)\hfill & (0.523,0.008)\hfill \end{array}$
UCAV4	$\begin{array}{ll}(1.000,0.042)\hfill & (0.849,0.011)\hfill \\ (0.874,0.042)\hfill & (0.923,0.011)\hfill \\ (0.628,0.084)\hfill & (1.000,0.009)\hfill \\ (0.432,0.062)\hfill & (0.943,0.011)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.044)\hfill & (0.901,0.011)\hfill \\ (0.766,0.078)\hfill & (0.932,0.011)\hfill \\ (0.871,0.067)\hfill & (1.000,0.010)\hfill \\ (0.724,0.066)\hfill & (0.986,0.011)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.095)\hfill & (0.860,0.011)\hfill \\ (0.964,0.083)\hfill & (0.946,0.011)\hfill \\ (0.933,0.120)\hfill & (1.000,0.010)\hfill \\ (0.917,0.071)\hfill & (0.935,0.011)\hfill \end{array}$	$\begin{array}{ll}(1.000,0.048)\hfill & (0.867,0.011)\hfill \\ (0.745,0.060)\hfill & (0.911,0.011)\hfill \\ (0.826,0.132)\hfill & (1.000,0.009)\hfill \\ (0.847,0.071)\hfill & (0.954,0.011)\hfill \end{array}$

.

Based on the normalized criteria value matrices ${\overline{S}}^{\left(i,k\right)}$ ($i=1,2,3,4;k=1,2,3,4$) in Table 5, we can apply the extended SMAA-2 method to compute the holistic acceptability index matrix $E={\left[e_{ij}\right]}_{M\times N}$. To better demonstrate our research, we consider two different preferences of commanders: (1) commanders are more concerned with the operational benefit, i.e., ${\omega }_B\ge {\omega }_C$ and (2) commanders are more concerned with the operational cost, i.e., ${\omega }_B<{\omega }_C$.

The holistic acceptability index matrices in the above two case are illustrated in Figure 5a,b, respectively.

Table 6. The holistic acceptability indices when ${\omega }_B\ge {\omega }_C$ and ${\omega }_B<{\omega }_C$.

${\mathit{e}}_{\mathit{i}\mathit{j}}$	${\mathit{\omega }}_{\mathit{B}}\ge {\mathit{\omega }}_{\mathit{C}}$				${\mathit{e}}_{\mathit{i}\mathit{j}}$	${\mathit{\omega }}_{\mathit{B}}<{\mathit{\omega }}_{\mathit{C}}$
	${\mathit{T}}_1$	${\mathit{T}}_2$	${\mathit{T}}_3$	${\mathit{T}}_4$		${\mathit{T}}_1$	${\mathit{T}}_2$	${\mathit{T}}_3$	${\mathit{T}}_4$
$UCA{V}_1$	0.9664	0.1899	0.5141	0.2496	$UCA{V}_1$	0.5689	0.1200	0.9511	0.2800
$UCA{V}_2$	0.9999	0.4003	0.3808	0.1390	$UCA{V}_2$	1.0000	0.3094	0.4906	0.1200
$UCA{V}_3$	0.9996	0.4901	0.2272	0.2031	$UCA{V}_3$	1.0000	0.5200	0.1216	0.2784
$UCA{V}_4$	0.9724	0.3728	0.4219	0.1530	$UCA{V}_4$	0.3646	0.3100	0.9176	0.3278

gives the specific values of two matrices. It may be noted that the holistic acceptability index $e_{ij}$ in Table 6 does not indicate the actual combat reward of $UCA{V}_i$ attacking $T_j$, but the relative value of the combat reward under the given distributions of operational objectives and the given distributions of operational objective weights. For $UCA{V}_i$, the larger $e_{ij}$, the more $UCA{V}_i$ should attack $T_j$ to achieve a larger combat reward. Based on the holistic acceptability index matrices in Table 6, task assignment schemes are obtained in the two cases, as shown in

Table 7. Task assignment schemes when ${\omega }_B\ge {\omega }_C$ and ${\omega }_B<{\omega }_C$.$x_{ij}$.

	${\mathit{\omega }}_{\mathit{B}}\ge {\mathit{\omega }}_{\mathit{C}}$				${\mathit{x}}_{\mathit{i}\mathit{j}}$	${\mathit{\omega }}_{\mathit{B}}<{\mathit{\omega }}_{\mathit{C}}$
	${\mathit{T}}_1$	${\mathit{T}}_2$	${\mathit{T}}_3$	${\mathit{T}}_4$		${\mathit{T}}_1$	${\mathit{T}}_2$	${\mathit{T}}_3$	${\mathit{T}}_4$
$UCA{V}_1$	0	0	0	1	$UCA{V}_1$	0	0	1	0
$UCA{V}_2$	1	0	0	0	$UCA{V}_2$	1	0	0	0
$UCA{V}_3$	0	1	0	0	$UCA{V}_3$	0	1	0	0
$UCA{V}_4$	0	0	1	0	$UCA{V}_4$	0	0	0	1

. Figure 6 shows the task assignment diagrams corresponding to the schemes in Table 7. As can be seen from Figure 5a and Table 6, when commanders are more concerned with the operational benefit, ${\omega }_B\ge {\omega }_C$, attacking $T_1$ can bring the largest holistic acceptability index for each UCAV. The reason for this can be found in Table 4 and Table 5. It can be seen from the second column of Table 4 that the average value of $T_1$ is the largest according to the target value information obtained by each sensor. Moreover, in Table 5, for each UCAV, the utilities (or normalized operational benefits) obtained by attacking $T_1$ are the largest, all being equal to 1. However, due to the constraints expressed by Equations (41) and (42) in Section 6, $T_1$ is assigned to $UCA{V}_2$ in order to maximize the holistic acceptability index of the task assignment scheme shown in the left part of Table 7.

When commanders are more concerned with the operational cost, ${\omega }_B<{\omega }_C$, the holistic acceptability indices of $UCA{V}_2$ attacking $T_1$ and $UCA{V}_3$ attacking $T_1$ are almost unchanged, both being equal to 1, which indicates that for both $UCA{V}_2$ and $UCA{V}_3$, attacking $T_1$ is the best choice. On the other hand, for $UCA{V}_1$ and $UCA{V}_4$, the best choice is attacking $T_3$ instead of $T_1$. The reason for the change is that commanders become more inclined to consider the operational cost in this case. Due to the same reason, the task assignment scheme varies, as shown in the right part of Table 7.

It can be seen from the above analysis that although the actual operational reward cannot be computed due to multiple uncertainties in the task assignment problem, we can still achieve multiobjective task assignment for multi-UCAV. Specifically, we have transformed the actual operational reward of assigning $T_j$ to $UCA{V}_i$ into the holistic acceptability index $e_{ij}$ of the assignment, as shown in Table 6. By maximizing the sum of the holistic acceptability indices of all assignments, we have found the optimal task assignment scheme. In addition, it can also be seen that as the preferences of commanders varies, the holistic acceptability index of the same assignment changes. For example, in both cases, $UCA{V}_1$ attacks $T_1$, indicating that the actual operational benefit and the actual operational cost have remained unchanged. However, due to the variation of the preferences of commanders, the utility values of the actual operational benefit and the actual operational cost have changed, leading to the holistic acceptability index of assigning $T_1$ to $UCA{V}_1$ falling from 0.9664 to 0.5689. A similar situation could also be seen in other assignments. Ultimately, the changes in the holistic acceptability indices of all assignments result in the adjustment of the task assignment scheme, which indicates that that the preferences of commanders have an important impact on operational decision-making.

In summary, the proposed method in this research is able to solve the multiobjective task assignment problem for multi-UCAV under multiple uncertainties and reveal, to some extent, the influence of the variation in the preferences of commanders on operational decision-making. In our research, we only considered two objectives: the operational benefit and the operational cost. However, more objectives can be considered, such as maximizing the operational benefit, minimizing the operational cost, and minimizing the operational loss of UCAVs. In this case, we would need to establish more complex models and acquire more target information.

8. Conclusions

In this paper, we investigated the task assignment problem for multi-UCAV under multiple uncertainties, including uncertain target information provided by multiple sensors, uncertain weights of the criteria or objectives concerned, and partially known weights of sensors. To find solutions for the problem, we extended the SMAA-2 method and then combined it with the integer linear programming method. Through a case study, the feasibility of the proposed method was demonstrated.

The most significant contribution of our research is the extension of the original SMAA-2 method. The new method possesses the following strengths and weaknesses:

• Due to the integration of the WAA operator-based aggregation method, the extended method is able to perform stochastic multicriteria acceptability analysis under circumstances of multiple criteria value matrices. This is also the most prominent advantage of the extended method over the original SMAA-2 method.
• The extended method employs the iterative algorithm to compute the objective weight of a sensor based on the accuracy and credibility of the target information it provides so as to minimize the subjective judgment errors of experts on the sensor performance. By this means, the actual importance of a sensor can be reflected to the greatest extent so that the final criteria value matrix obtained could be as accurate as possible. Through these measures, the credibility and accuracy of task assignment schemes can be further increased.
• The foundation of the extended method is the stochastic multicriteria acceptability analysis, so it inherits the advantage of the SMAA technique in that it does not require decision-makers to express their preferences explicitly [30]. This property is particularly useful in operational planning because the commander could only express implicit preferences for multiple operational objectives most of the time.
• The extended method is based on the assumption that the uncertain target information can be represented by normal distribution interval numbers. This assumption may limit the application of the method. When the uncertain target information is represented by fuzzy sets or rough sets, the extended method cannot be applied.
• The extended method employs a linear utility value function to map criteria values to utility values. Although Lahdelma et al. have proved it may often be safe to apply a linear utility value function, they also mentioned that when the convexity or concavity of the utility value function is severe, the results would be slightly different [45]. This means that in some special cases, the linear utility value function may not be the best choice, and we would need to construct a more appropriate one.

In future research, we will focus on improving the weaknesses listed above in the extended SMAA-2 method so that the combination method can be used to solve similar multiobjective resource allocation problems under serve uncertainties in other fields.