Uncertain multi-objective dynamic weapon-target allocation problem based on uncertainty theory

1.   Introduction

The weapon target allocation (WTA) problem is one of the crucial subjects for the automation military command and control system. The earliest research on WTA problem introduced by Manne ^[1] provides the basis WTA model, which is a classical constrained combinatorial optimization problem and has been proved to be NP-hard. In Manne's model, maximizing the probability of damage to the target is the optimal goal, and time is not a dimension considered. Since then, various models of WTA problem have been proposed and can be divided into two distinct categories: the Static WTA (SWTA) ^[2,3,4,5,6] and the Dynamic WTA (DWTA) ^{[7,8,9,10,11]}.

Most of the literature on the WTA problem focuses on optimizing a single objective, such as minimizing the operation cost ^[5], maximizing the destruction of hostile targets ^[6], minimizing the destruction of protected assets ^[9,10]. However, in practical terms, the commander usually considers more than one objective to evaluate the allocation scheme comprehensively. For example, focusing only on maximizing the destruction of targets may result in a waste of weapon ammunitions, which may lead to insufficient firepower in subsequent engagement. If multiple objectives are considered, the WTA problem becomes a multi-objective programming (MOP) problem, which generally does not have an absolutely optimal solution but a set of Pareto efficient solutions. For example, Leboucher et al. ^[12] considered three different criteria as objectives, but they transform the multi-objective functions to a single-objective function by linear weight method and get the optimal allocation scheme under a preference of decision maker, which cannot comprehensively obtain the optimal allocation schemes under various preferences. Xu ^[13] proposed a WTA model with two objective functions and obtained a set of Pareto optimal solutions, but the model is simplified. For solving the WTA problem with multiple objectives, the multi-objective evolutionary algorithm (MOEA) is an effective way. Schaffer ^[14] designed the VEGA, which is regarded as the earliest MOEA. So far, many outstanding MOEAs have been proposed. Deb ^[15] proposed the NSGA-II, which is the widely used MOEAs that use fast non-dominated sorting and the crowding distance. Recently, the MOEAs based on decomposition (MOEA/D) have become increasingly popular. They are based on the principle of decomposing a MOP problem into a set of subproblems and optimizing them simultaneously. Zhang ^[16] proposed the first version of MOEA/D. After that, many scholars have improved the decomposition method ^[17,18] and adaptive mechanism ^[19,20,21]. Shi ^[22] utilized the NSGA-II to solve a WTA problem with two objectives. Li ^[23] utilized the NSGA-II and MOEA/D to solve a WTA problem with two objectives and made a comparison. The performance of MOEA is very sensitive to the parameters involved in it, so how to design efficient algorithms based on the framework of MOEAs to solve scenario-based WTA problems is a research key.

For the literature about WTA problem mentioned above, a common assumption is that the battlefield environment is deterministic, and all the parameters in these models are known. However, it is difficult to meet in practice. On one hand, the electromagnetic interference and decoys released by the opponent make it difficult to detect the target information and estimate the target threat value accurately. On the other hand, the highly maneuverable flights of targets make it hard to assess the target attack intention correctly. What's more, some indicators, such as strategic value of protected assets and the cost of human resources, both rely heavily on human experience. Thus, the WTA problem is often indeterministic with vague or imprecise statements, which can be boiled down to problems of observing the parameters themselves, deficiency in history and statistical data and the subjectivity of human judgement, etc. Some scholars consider these uncertain factors as random variables, and then the probability theory is used to deal with them. For example, Krokhmal ^[24] utilized the CVaR constraint to study the uncertainty in the WTA problem. Ahner ^[25] regarded the number of targets for subsequent engagement as a random variable, and a two-stage stochastic programming model was proposed. Li ^[26] proposed a robust optimization model for the uncertainty of interception probability.

However, the application of probability theory is based on the law of large numbers, i.e., it needs a large scale of historical sample data as the support of credibility of results, which is difficult to meet in the indeterminate battlefield environment. For this case, we can invite experts to give an estimation and offer a distribution function for each uncertain parameter, but some surveys have shown that human beings usually overestimate the occurrence degree of unlikely events, and human belief degree may have larger variance than the real frequency of events, so if the probability theory is still used to deal with them, the conclusion may be a paradox ^[27]. In order to better describe the subjective imprecise quantity, Liu ^[28] proposed the uncertainty theory, which is a branch of mathematics based on normality, duality, subadditivity, product axioms to deal with imprecision. To date, many scholars have studied the uncertainty theory and established a complete mathematical system including the uncertain set ^[29], uncertain process ^[30] and uncertain differential equation ^[31]. Meanwhile, it has been studied in many application fields, like uncertain inference ^[29], uncertain programming ^[32], uncertain statistics ^[33], uncertain portfolio selection ^[34], uncertain optimal control ^[35]. So far, there are few relevant studies on WTA problem based on uncertainty theory.

In this paper, we aim to study the DWTA problem under uncertainty in a defensive scenario. The main contributions of this paper can be summarized as follows. Firstly, we use the uncertainty theory to deal with the uncertain factors in indeterminate battlefield environment and assume that the threat value of targets, the value of protected assets and the extra cost of operation are uncertain variables, and an uncertain multi-objective dynamic WTA (UMDWTA) model is established originally, where the value of destruction of targets, the value of surviving assets and the cost of operation are the objective functions. To make the proposed model solvable, a deterministic model named $E\sigma$-UMDWTA is derived based on expected value and standard deviation principle, which considers the minimum average cost and the minimum fluctuation of the uncertain objective functions simultaneously. Second, a virtual representation (VP) and construction procedure (CP) of feasible solution for UMDWTA problem are proposed, which can facilitate the generation of the feasible solution efficiently. Third, on the basis of VP and CP, an improved MOEA/D with adaptive weights (MOEA/D-AW) is designed to obtain the Pareto efficient solutions of the proposed model, where three new evolutionary operators are proposed to evolve the VP of feasible solution, and the adaptive weights mechanism is employed to improve the distribution uniformity of the Pareto front obtained in the iteration process. The test results show that the designed algorithm can obtain the Pareto fronts with high uniformity and convergence in solving MOP problems. Finally, a case of the UMDWTA problem is solved by the improved MOEA/D-AW, then various Pareto efficient solutions and corresponding allocation schemes are obtained, and the simulation results also verify the effectiveness of the improved MOEA/D-AW and feasibility of the UMDWTA model.

The paper is organized as follows. In Section 2, some basic knowledge of uncertainty theory and multi-objective programming is reviewed. In Section 3, the UMDWTA model is established. In Section 4, the solution method of UMDWTA is analyzed. In Section 5, the improved MOEA/D-AW is designed, and its performance of solving MOP problem is tested by four test functions. In Section 6, a case of the UMDWTA problem is solved by the designed algorithm, and the simulation results are analyzed. Finally, the main results of this paper are concluded in Section 7.

2.   Preliminary

2.1.   Uncertainty theory

Let $\Gamma$ be a nonempty set, and ${\mathscr L}$ is a $\sigma$-algebra over $\Gamma$. Each element $\Lambda$ in ${\mathscr L}$ is called an event. A set function ${{\mathscr M}}$ from ${\mathscr L}$ to [0, 1] is called an uncertain measure if it satisfied the following axioms ^[28]:

Axiom 1. (Normality axiom) ${\mathscr M}\{ \Gamma \} = 1$ for the universal set $\Gamma$;

Axiom 2. (Duality axiom) ${\mathscr M}\{ \Lambda \} + {\mathscr M}\{ {\Lambda ^c}\} = 1$ for any event $\Lambda \in {\mathscr L}$;

Axiom 3. (Subadditivity axiom) For every countable sequence of events ${\Lambda _1}, {\Lambda _2}, ... \in {\mathscr L}$, we have

The triplet $(\Gamma, {\mathscr L}, {\mathscr M})$ is called an uncertainty space. Furthermore, a product uncertain measure is as follows.

Axiom 4. (Product axiom) Let $({\Gamma _i}, {{\mathscr L}_i}, {{\mathscr M}_i})$ be the uncertainty space for $i = 1, 2, ...$. The product uncertain measure ${\mathscr M}$ is an uncertain measure satisfying

Definition 1. ^[28] An uncertain variable is a measurable function $\xi$ from an uncertainty space $(\Gamma, {\mathscr L}, {\mathscr M})$ to the set of real numbers, i.e., for any Borel set $B$ of real numbers, the set

is an event in ${\mathscr L}$.

Definition 2. ^[36] The uncertain variables ${\xi _1}, {\xi _1}, ..., {\xi _n}$ are said to be independent if

for any Borel sets ${B_1}, {B_2}, ..., {B_n}$ of real numbers.

Theorem 1. ^[28] Let ${\xi _1}, {\xi _1}, ..., {\xi _n}$ be uncertain variables, and $f$ is a real-valued measurable function. Then, $f({\xi _1}, {\xi _1}, ..., {\xi _n})$ is an uncertain variable.

Definition 3. ^[28] The uncertainty distribution $\Phi$ of an uncertain variable $\xi$ is defined by

for any real number $x \in R$.

Definition 4. ^[36] Let $\xi$ be an uncertain variable with regular uncertainty distribution $\Phi$, and then the inverse function ${\Phi ^{ - 1}}$ is called the inverse uncertainty distribution of $\xi \ $.

Definition 5. ^[28] An uncertain variable $\xi$ is called linear if it has a linear uncertainty distribution

denoted by $\xi \sim z(a, b, c)$ where $a, b, c$ are real numbers within $a < b < c$.

Definition 6. ^[28] An uncertain variable $\xi$ is called zigzag if it has a zigzag uncertainty distribution

denoted by $\xi \sim z(a, b, c)$ where $a, b, c$ are real numbers within $a < b < c$.

Definition 7. ^[28] An uncertain variable $\xi$ is called exponential if it has a exponential uncertainty distribution

denoted by $x \sim E(\lambda)$ where $\lambda$ is a constant within $\lambda > 0$.

Theorem 2. ^[33] Let ${\xi _1}, {\xi _2}, ..., {\xi _n}$ be independent uncertain variables with regular uncertainty distributions ${\Phi _1}, {\Phi _1}, ..., {\Phi _n}$, respectively. If the function $f({x_1}, {x_2}, ..., {x_n})$ is a measurable function which is strictly increasing with respect to ${x_1}, {x_2}, ..., {x_m}$ and strictly decreasing with respect to ${x_{m + 1}}, {x_{m + 2}}, ..., {x_n}$, then ${\xi} = f({\xi _1}, {\xi _2}, ..., {\xi _n})$ is an uncertain variable with inverse uncertainty distribution

Definition 8. ^[28] Let $\xi$ be an uncertain variable, and then the expected value of $\xi$ is defined by

provided that at least one of the two integrals is finite.

Theorem 3. ^[28] Let $\xi \ $ be an uncertain variable with regular uncertainty distribution $\Phi$. If the expected value exists, then

Definition 9. ^[28] Let $\xi \ $ be an uncertain variable with finite expected value $e$, and then the variance of $\xi$ is defined by

Theorem 4. ^[28] Let $\xi \ $ be an uncertain variable with regular uncertainty distribution $\Phi$ and finite expected value $e$. Then,

The above contents are some basic definitions and theorems of uncertainty theory, which are used subsequently.

2.2.   Multi-objective programming

A programming problem with multiple optimization objectives is called a multi-objective programming (MOP) problem. The general MOP model with $m$ objective functions is represented as follows:

where $\mathit{\boldsymbol{x}} = {({x_1}, {x_2}, ..., {x_n})^T}$ is the decision vector; $D$ denotes the feasible region.

Then, some related definitions of MOP model are introduced here.

Definition 10. ^[37] (Absolutely optimal solution) Let $\tilde {\mathit{\boldsymbol{x}}} \in {{D}}$, and if it satisfies

for any given $\mathit{\boldsymbol{x}} \in {{D}}$, then the solution $\tilde {\mathit{\boldsymbol{x}}}$ is called the absolute optimal solution of the MOP problem (2.1).

Definition 11. ^[37] (Pareto dominance) For a given vector $\mathit{\boldsymbol{y_1}} = ({y_{11}}, {y_{12}}, ..., {y_{1m}})\in \mathit{\boldsymbol{R}}^m$, it dominates the vector $\mathit{\boldsymbol{y_2}} = ({y_{21}}, {y_{22}}, ..., {y_{2m}})\in {R}^m$ if and only if

which is denoted as $\mathit{\boldsymbol{y}}_1 \prec \mathit{\boldsymbol{y}}_2$.

Definition 12. ^[37] (Non-dominated solution) For a given set ${\Gamma} \subset \mathit{\boldsymbol{R}}^m$ and $\mathit{\boldsymbol{y}} \in {\Gamma}$, if there is no other vector $\mathit{\boldsymbol{y}} ' \in {\Gamma}$ such that $\mathit{\boldsymbol{y}} ' \prec {\mathit{\boldsymbol{ y }}}$, then the vector $\mathit{\boldsymbol{y}}$ is called a non-dominated solution in the set ${\Gamma}$.

Definition 13. ^[37] (Pareto efficient solution) For a given MOP model (2.1), the vector ${\mathit{\boldsymbol{x}}^ * } \in {D}$ is called a Pareto efficient solution if there is no ${{\mathit{\boldsymbol{ x }}}} \in { D}$ such that $f(\mathit{\boldsymbol{x}}) \prec f({\mathit{\boldsymbol{x}}^ * })$.

Definition 14. ^[37] (Pareto optimal set) For a given MOP model (2.1), the Pareto optimal set ${{D}^ * }$ is defined as follows:

Definition 15. ^[37] (Pareto front) For a given MOP model (2.1), the Pareto front is defined as follows:

3.   Establishment of the UMDWTA model

This section describes a scenario of the air defense combat operation and establishes an uncertain multi-objective dynamic WTA (UMDWTA) model.

3.1.   Problem description

The following scenario is considered: There are $K$ ground assets to be protected (protected assets). At one point, the radars detect some hostile targets, and they are going to attack the protected assets. There are $W$ weapons that can launch missiles to intercept these targets. Since there is still a period of time before the targets break through the defense, the interception process can be divided into $S$ stages, where one stage is the minimum combat time unit. The commander can make a decision to allocate weapons to intercept targets in each stage. Figure 1 shows the scenario description.

However, due to the various uncertainties in indeterminate battlefield environment discussed in Section 1, some parameters are uncertain, like the attack intention of each target, the threat value of each target, the strategic value of each asset, the extra cost of operation.

3.2.   Objective functions of UMDWTA

Different evaluation criteria result in different weapon-target allocation schemes. Here we employ three criteria as the objectives of UMDWTA model: the value of destruction of targets, the value of surviving assets and the cost of operation. For ease of reading, Table 1 summarizes some notations that are used throughout this paper.

3.2.1.   The value of destruction of targets

Destroying the incoming hostile targets is the most direct goal in the defense, so the value of destruction of targets should be employed as an objective function, which is the most widely employed combat effectiveness evaluation criteria in relevant studies ^[7,8].

Generally, the interception of each weapon at each stage is independent of each other, so the degree of each target being destroyed after the final stage is equal to the product of the degrees to which each weapon destroys it at each stage. Meanwhile, since the state of each target is different, each target has its own threat value, denoted by ${\xi _{T, j}}$, and the destroying of the target with higher threat value contributes more to the objective. The formulation of this objective for current stage $t$ is shown as follows:

where ${{\bf{X}}^t} = [{\mathit{\boldsymbol{X}}_t}, {\mathit{\boldsymbol{X}}_{t + 1}}, ..., {\mathit{\boldsymbol{X}}_S}]$ with ${\mathit{\boldsymbol{X}}_k} = {[{x_{ij}}(k)]_{W(k) \times T(k)}}$ is the decision matrix at stage $k$ (${x_{ij}}(k) = 1$ if weapon $i$ is assigned to target $j$ at stage $k$; ${x_{ij}}(k) = 0$ otherwise); ${p_{ijh}}$ is the degree that weapon $i$ successfully intercept target $j$ at stage $t$; $T(t)$ is the number of existing targets at current stage $t$; ${\xi _{T, j}}$ is the threat value of target $j$.

In previous models, the threat value of each target is determined by many factors, e.g., the target type, the attack angle, the velocity of target, the target range ^[38]. However, in the indeterminate battlefield environment discussed above, it is difficult to obtain the state information of targets accurately, so we can consider the threat value of each target as an uncertain variable and invite experts to subjectively estimate them based on the incomplete information. Here, the experts give the minimum value and maximum value of it, thus ${\xi _{T, j}}$ is considered as a linear uncertain variable, i.e., ${\xi _{T, j}} \sim L({a_{T, j}}, {b_{T, j}})$, and its uncertainty distribution is

where ${{a_{T, j}}}$ and ${{b_{T, j}}}$ are the maximum and minimum threat value of target $j$, respectively, within $0 \le {a_{T, j}} < {b_{T, j}}$.

3.2.2.   The value of surviving assets

The commander usually wants to avoid the destruction of protected assets as much as possible, so the value of surviving protected assets after the final stage can be employed as an objective function, which is employed as the combat effectiveness evaluation criterion in relevant studies ^[9,10].

If a target survives after the final stage, we say that the target breaks through the defense of weapons, and then it is able to attack the protected assets. Therefore, the degree that asset $k$ survives (is not attacked by any target) after the final stage is

where ${\xi _{A, jk}}$ is the degree that target $j$ successfully attacks asset $k$ if it breaks through the defense of weapons.

The surviving of asset that has higher strategic value contributes more to the objective, so the formulation of the objective for current stage $t$ is shown as follows:

where ${\xi _{V, k}}$ is the strategic value of asset $k$.

However, the strategic value of each protected asset is determined by subjective factors, which has been discussed in Section 1, so it is reasonable to consider it as an uncertain variable and invite experts to eastimate it. Here, the experts give the maximum and minimum value of it, thus ${\xi _{V, k}}$ is a linear uncertain variable, i.e., ${\xi _{V, k}} \sim L({a_{V, k}}, {b_{V, k}})$, and its uncertainty distribution is

where ${{a_{V, k}}}$ and ${{b_{V, k}}}$ are the maximum and minimum value of asset $k$, respectively, within $0 \le {a_{V, k}} < {b_{V, k}}$.

Meanwhile, the degree that target $j$ successfully attacks asset $k$, denoted by ${\xi _{A, jk}}$, is affected by many uncertain factors, e.g., the attack intention of target, the type of target and the amount of ammunition carried by the target, so we can consider it as an uncertain variable and invite experts to estimate it. Here, the experts give the maximum, average and minimum value of it, thus ${\xi _{A, jk}}$ is a zigzag uncertain variable, i.e., ${\xi _{A, jk}} \sim z({a_{A, jk}}, {b_{A, jk}}, {c_{A, jk}})$, and its uncertainty distribution is

where ${{a_{A, jk}}}$, ${{b_{A, jk}}}$ and ${{c_{A, jk}}}$ are the maximum, average and minimum degree that target $j$ destroys asset $k$, respectively, within $0 \le {a_{A, jk}} < {b_{A, jk}} < {c_{A, jk}} \le 1$.

3.2.3.   The cost of operation

Only considering the maximization of the value of destruction of targets and the value of surviving protected assets may result in the waste of resources, and thus the cost of operation in all stages should be taken as an objective function.

The cost of operation mainly refers to the cost of missiles that all weapon systems launch to intercept the targets and other extra costs, e.g., the wastage of equipment, the cost of human resources. Here, we assume that a weapon can only launch one missile in a stage, so the formulation of the objective for current stage $t$ is shown as follows:

where ${{C_i}}$ is the cost of one missile of weapon $i$; ${\xi _{C}}$ is the extra cost per interception.

Since the extra cost is hard to predict objectively, we consider it as an uncertain variable and invite the experts to estimate it. Here, the experts regard it as an exponential uncertain variable of parameter $\lambda$, i.e., ${\xi _{C}} \sim E(\lambda)$, and its uncertainty distribution is

where $\lambda$ is a real number within $\lambda > 0$.

3.3.   Constraints

The constraints in the above UMDWTA problem mainly include the following four categories:

Inequation (3.4) limits the number of targets that a weapon can engage at each stage. Most actual weapons can shoot only one target at a stage, so ${n_i}$ is usually set to be 1. Besides, a special weapon that is capable of engaging multiple targets simultaneously can be viewed as multiple separate weapons. Inequation (3.5) ensures that at most $m_j$ weapons can be allocated to target $j$ at a stage, and $m_j$ usually depends primarily on weapon performance, firing strategy and the threat value of target. Inequation (3.6) reflects that the amount of ammunitions, as used by a certain weapon $i$ through all stages, cannot exceed its predefined allowable amount $N_i$. Inequation (3.7) is a very important constraint to the actual dynamic WTA problem which reflects the influence of time windows on the engagement feasibility of weapons. If weapon $i$ can be allocated to target $j$ in stage $h$, then ${f_{ijh}} = 1$; otherwise, ${f_{ijh}} = 0$.

3.4.   Mathematical model

Based on the objective functions and constraints analyzed above, the mathematical formulation of the UMDWTA problem is as follows:

Note that the maximization of ${f_{d}}({{\bf{X}}^t}, {\mathit{\boldsymbol{\xi}}})$ and ${f_{s}}({{\bf{X}}^t}, {\mathit{\boldsymbol{\xi}}})$ are equivalent to the minimization of $-{f_{d}}({{\bf{X}}^t}, {\mathit{\boldsymbol{\xi}}})$ and $-{f_{s}}({{\bf{X}}^t}, {\mathit{\boldsymbol{\xi}}})$, respectively. For simplicity, let $G$ denote the feasible region of the UMDWTA model above.

Remark 1. The decision variable matrix ${\bf{X}}^t$ represents the total stages allocation scheme for the current stage $t$. If the defender observes that some targets have been destroyed or some other dynamic event happened at the current stage, then the parameters of the UMDWTA model need to be reevaluate for the next stage $t+1$. In this case, the weapons need to be reallocated to targets, and a new decision variable matrix ${\bf{X}}^{t+1}$ will be generated. Otherwise, the decision variable matrix remains the same.

4.   Solution method of the UMDWTA model

The model (3.8) contains uncertain variables and multiple objective functions, and it is an unsolvable model because there is no natural ordered relation in uncertain space. In this section we introduced three methods which can convert the unsolvable UMDWTA model to a deterministic one.

Different real-life problems call for different meanings of valuation to satisfy its needing practical application, then result in different compromise decisions ^[39]. Since we want to minimize the average cost of objective functions, it is rational to use the expected value model, which can be represented as follows:

where $E(\cdot)$ denotes the expected value.

What's more, the stability and volatility are usually important, so minimizing the variance of objective functions should be considered, which can be represented as follows:

where $V (\cdot)$ denotes the variance.

Sometimes both expected value and standard deviation need to be taken into account simultaneously, so the following model can be established:

where $\sigma (\cdot)$ denotes the standard deviation, i.e., $\sigma (\cdot) = \sqrt {V(\cdot)}$.

Problems (4.1)–(4.3) are called the $E$-UMDWTA model, $V$-UMDWTA model and $E\sigma$- UMDWTA model, respectively. For simplicity, we use ${\mathit{\boldsymbol{ x }}}$ to replace ${{\bf{X}}^t}$ in subsequence. Note that the validity of the above models is based on the assumption that the objective functions are uncertain variables, and their variances are finite.

Through any one of the above three methods, the original uncertain model (3.8) can be converted to a deterministic one. In this paper, we choose the $E\sigma$-UMDWTA model for analysis since it considers the minimum average cost and the minimum fluctuation of the uncertain objective functions simultaneously. Hence, the Pareto efficient solution of $E\sigma$-UMDWTA model can be defined:

Definition 16. Let $\tilde {\mathit{\boldsymbol{x}}} \in {G}$, the solution $\tilde {\mathit{\boldsymbol{x}}}$ is called the expected-value standard deviation efficient solution of the UMDWTA model if it is a Pareto efficient solution of the $E\sigma$-UMDWTA model, that is, there is no $\mathit{\boldsymbol{x}} \in G$ such that $\left({E\left({\mathit{\boldsymbol{f}}({{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right), \sigma \left({\mathit{\boldsymbol{f}}({{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right)} \right) \prec \left({E\left({\mathit{\boldsymbol{f}}(\tilde {{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right), \sigma \left({\mathit{\boldsymbol{f}}(\tilde {{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right)} \right)$.

The expected-value efficient solution of $E$-UMDWTA and the variance efficient solution of $V$-UMDWTA can be similarly defined, respectively. Further, the relation among these three efficient solutions is given by Theorem 5.

Theorem 5. Let ${G_1}$ and ${G_2}$ be the set of expected-value efficient solutions of $E$-UMDWTA model and the set of variance efficient solutions of $V$-UMDWTA model, respectively, and let ${G_3}$ be the set of expected-value standard deviation efficient solutions of $E\sigma$-UMDWTA model, so we have ${G_1} \cap {G_2} \subset {G_3}$.

Proof. If ${G_1} \cap {G_2} = \emptyset$, the conclusion is clearly valid. If ${G_1} \cap {G_2} \ne \emptyset$, we prove it by contradiction: Let ${\mathit{\boldsymbol{ x }}}$ be both the expected-value efficient solution of $E$-UMDWTA model and the variance efficient solution of $V$-UMDWTA model, i.e., $\mathit{\boldsymbol{x}} \in {G_1} \cap {G_2}$. If ${\mathit{\boldsymbol{ x }}} \notin {{G_3}}$, then there is $\bar {\mathit{\boldsymbol{x}}} \in {G}$ such that $\left({E\left({\mathit{\boldsymbol{f}}(\bar{{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right), \sigma \left({\mathit{\boldsymbol{f}}({\bar{{\mathit{\boldsymbol{ x }}}}}, \mathit{\boldsymbol{\xi}})} \right)} \right) \prec \left({E\left({\mathit{\boldsymbol{f}}({{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right), \sigma \left({\mathit{\boldsymbol{f}}({{\mathit{\boldsymbol{ x }}}}, \mathit{\boldsymbol{\xi}})} \right)} \right)$, i.e.,

and there is ${i_0}({i_0} \in \{ 1, 2, 3\})$ satisfying

If $E[{f_{i_0}}(\bar {\mathit{\boldsymbol{x}}}, \mathit{\boldsymbol{\xi}})] < E[{f_{i_0}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})]$, ${\mathit{\boldsymbol{ x }}}$ is not an expected-value efficient solution of $E -$UMDWTA model, i.e., $\mathit{\boldsymbol{x}} \notin {G_1}$, which contradicts the prerequisite that $\mathit{\boldsymbol{x}} \in {G_1} \cap {G_2}$.

Otherwise, $\sigma[{f_{i_0}}(\bar {\mathit{\boldsymbol{x}}}, \mathit{\boldsymbol{\xi}})] < \sigma[{f_{i_0}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})]$, since $\sigma [{f_{i_0}}(\bar {\mathit{\boldsymbol{x}}}, \mathit{\boldsymbol{\xi}})] > 0$ and $\sigma [{f_{i_0}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})] > 0$, we have $V[{f_{i_0}}(\bar {\mathit{\boldsymbol{x}}}, \mathit{\boldsymbol{\xi}})] < V[{f_{i_0}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})]$, and ${\mathit{\boldsymbol{ x }}}$ is not a variance efficient solution of $V$-UMDWTA model, i.e., $\mathit{\boldsymbol{x}} \notin {G_2}$, which also contradicts the prerequisite that $\mathit{\boldsymbol{x}} \in {G_1} \cap {G_2}$. The theorem is proved.

Furthermore, we can assume that the uncertain variables in model (3.8) are all independent since there's no coupling between them, so the calculation formulas of $E\sigma$-UMDWTA model (4.3) can be derived subsequently.

For any given $\mathit{\boldsymbol{x}} \in \{ \mathit{\boldsymbol{x}}\left| {\mathit{\boldsymbol{x}} \in {G}, \mathit{\boldsymbol{x}} \ne \bf{0}} \right.\}$, the functions $- {f_{d}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})$, $- {f_{s}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})$ and ${{f_{c}}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}})$ are all real-value measurable functions of ${\mathit{\boldsymbol{\xi}}}$, so we can know that these three objective functions are all uncertain variables through Theorem 2.1. Additionally, we can get the following three properties:

$1)$ $- {f_{d}}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{\xi}}})$ is strictly decreasing with respect to ${\xi _{T, j}}$ for any $j = 1, 2, ..., T(t)$;

$2)$ $- {f_{s}}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{\xi}}})$ is strictly decreasing with respect to ${\xi _{V, k}}$ for any $k = 1, 2, ..., S$, and it is strictly increasing with respect to ${\xi _{A, jk}}$, for any $j = 1, 2, ..., T(t)$ and any $k = 1, 2, ..., S$;

$3)$ ${f_{c}}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{\xi}}})$ is strictly increasing with respect to ${\xi _{cos }}$.

Let $\Psi _{d}^{ - 1}(\alpha)$, $\Psi _{s}^{ - 1}(\alpha)$ and $\Psi _{c}^{ - 1}(\alpha)$ denote the inverse uncertainty distributions of $- {f_{d}}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{\xi}}})$, $- {f_{s}}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{\xi}}})$ and ${f_{c}}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{\xi}}})$, respectively, and then through Theorem 2.2, their calculation formulas are as follows:

Hence, the expected value and standard deviation of objective functions can be calculated through Theorem 3 and Theorem 4, respectively. The formulation is as follows:

To sum up, the original UMDWTA model (3.8) is converted to a determinate MOP model (4.3) called $E\sigma$-UMDWTA model and the calculation formulas (4.7)–(4.12) are derived. However, it is difficult to obtain the Pareto front and Pareto solutions of $E\sigma$-UMDWTA model analytically, so it is needed to design an efficient multi-objective evolutionary algorithm (MOEA) to solve it.

5.   Improved MOEA/D-AW for solving UMDWTA problem

5.1.   Virtual representation and construction procedure of feasible solution

The decision variable ${\bf{X}}^t$ of UMDWTA is a binary matrix with high dimensions and is restricted by inequations (3.4)–(3.7), so it is difficult to directly generate a feasible solution in the iteration process. To overcome this case, a virtual representation (VP) of UMDWTA solutions and its feasible solution construction procedure (CP) is proposed, which can generate a feasible solution quickly and efficiently.

First, the virtual representation is the permutation of the available allocation pairs which are defined as follows:

Definition 17. An available allocation pair denoted by $i$-$j$-$h$ is called an AAP, which indicates that weapon $i$ is assigned to target $j$ in stage $h$.

Example 1. Given an UMDWTA problem with two weapons, two targets and two stages, all the AAPs of it consist of following pairs: (1) $1$-$1$-$1$; (2) $1$-$1$-$2$; (3) $1$-$2$-$1$; (4) $1$-$2$-$2$; (5) $2$-$1$-$1$; (6) $2$-$1$-$2$; (7) $2$-$2$-$1$; (8) $2$-$2$-$2$.

Then, any permutation of them, like (2)-(1)-(5)-(7)-(3)-(8)-(4)-(6), is called a virtual representation (VP) of a certain UMDWTA solution. The foundation on which we employ the VP is that the generation of UMDWTA solution depends on the allocation order of AAPs. In fact, the permutation-based representation is common for the quadratic allocation problem (QAP) which is a classical NP-hard combinatorial optimization problem ^[40]. However, a certain permutation of AAPs for UMDWTA does not produce a feasible solution directly since some AAPs may cause violation of constraints. For this matter, several auxiliary variables are firstly defined which can record the state information currently, and they are used to judge which AAPs in the permutation cause violation to constraints. The auxiliary variables are defined as follows:

$1)$ ${\bf{N}_{WS}} = {[{n_{WS}}(i, h)]_{W \times S}}$: ${n_{WS}}(i, h)$ denotes the number of targets that weapon $i$ is assigned to in stage $h$;

$2)$ ${\bf{N}_{TS}} = {[{n_{TS}}(j, h)]_{T \times S}}$: ${n_{TS}}(j, h)$ denotes the number of weapons that are assigned to target $j$ in stage $h$;

$3)$ ${\bf{N}_{W}} = {[{n_{W}}(i)]_{1 \times W}}$: ${n_{W}}(i)$ denotes the number of times that weapon $i$ is assigned to the targets in all stages.

With the variables defined above, the rules for handling the constraints of inequations (3.4)–(3.6) are presented as follows:

$1)$ If ${n_{WS}}(i, h) = n_i$, weapon $i$ will not be assigned to other targets any more in stage $h$;

$2)$ If ${n_{TS}}(j, h) = m_j$, no more weapons will be assigned to target $j$ in stage $h$;

$3)$ If ${n_{W}}(i) = N_i$, weapon $i$ will not be assigned to any target in any stage.

For inequation (3.7), we can handle it directly by the engagement feasibility matrix ${\bf{F}}$, i.e., if $f_{ijh} = 0$, the AAP $i-j-h$ is not allowed. Then, the pseudocode of feasible solution CP is shown in Algorithm 1.

In the process of adding AAPs to the allocation scheme following the order of permutation, the auxiliary variables will be updated dynamically. If the updated auxiliary variables after adding an AAP violate any constraints (contradict the rules), this AAP will be unallowed and marked as an unassigned AAP (UAAP); otherwise, it will be allowed and marked as an assigned AAP (AAAP). It is easy to prove that no violation of constraints will occur due to the mechanism. Hence, for a certain permutation of AAPs, a feasible solution of UMDWTA can be generated by the CP. To the contrary, it is easy to prove that any feasible solution of UMDWTA can be generated by using the CP to a certain but not unique permutation of AAPs.

Additionally, it's easy to get the following two properties:

Property 1. Given a certain VP, there is no conflict among the AAAPs contained in it, i.e., the generated feasible solution ${\bf {X}}^t$ will not be changed by only rearranging the order of AAAPs.

Property 2. Given a certain VP, for a set of consecutive UAAPs following behind an AAAP, the generated feasible solution ${\bf {X}}^t$ will not be changed by only rearranging the order of these UAAPs.

In the following content of the algorithm design, we use VP and CP to represent and generate the feasible solutions of UMDWTA problem in the iterative process. What's more, since the collection of all VPs for a UMDWTA problem is too large to search efficiently, the search scope is limited to the collection of VPs consisting of all AAPs. This restriction is justified because the VP consisting of all AAPs has larger values of $E(-f_d(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}}))$ and $E(-f_s(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{\xi}}))$ than some VPs that do not consist of all AAPs, which is consistent with goal that decision makers generally prefer a larger destruction of targets or more surviving assets.

5.2.   The framework of improved MOEA/D

Multi-objective evolutionary algorithm based on decomposition (MOEA/D) ^[16] is one of the most popular algorithms among the MOEAs, which search for Pareto efficient solutions by decomposing the original MOP problem into a series of single objective optimization subproblems and optimizing each subproblem. Specially, when solving MOP problem with many objectives, MOEA/D has not the problem of "selection pressure" that NSGA-II has. There are several common decomposition approaches including weighted sum approach, Tchebycheff (TCH) approach, angle-based approach, boundary intersection (BI) approach ^[17], etc.

In this section, we modify the general MOEA/D framework with Ideal-Nadir TCH approach and design an improved algorithm to solve the UMDWTA problem efficiently, in which three new evolutionary operators are proposed to generate more promising offspring in the iteration process, and the adaptive weights mechanism is employed to increase the uniformity of the Pareto front adaptively. The framework of the improved MOEA/D is given in Algorithm 3.

5.2.1.   Population initialization

First, since the Ideal-Nadir TCH approach ^[41] can overcome the defects of the weighted sum approach caused by the dimension variance of objective functions and is easier to implement, the original MOP problem is decomposed to a series of subproblems by it as follows:

where ${\mathit{\boldsymbol{z}}^ * } = (z_1^ *, z_2^ *, ..., z_m^ *)$ and ${\mathit{\boldsymbol{z}}^{nad}} = ((z_1^{nad}, z_2^{nad}, ..., z_m^{nad}))$ are the ideal point and nadir point, respectively, i.e., ${{z_j}^ * } = \min \{ {f_j}(\mathit{\boldsymbol{x}})\left| {\mathit{\boldsymbol{x}} \in {G}} \right.\}$ and ${{z_j}^ * } = \max \{ {f_j}(\mathit{\boldsymbol{x}})\left| {\mathit{\boldsymbol{x}} \in {G}} \right.\}$ for any $j = 1, ..., m$; $\mathit{\boldsymbol{\lambda}} = ({\lambda _1}, {\lambda _2, ..., {\lambda _m}})\in {R^m}$ is the weight vector satisfying $\sum\nolimits_{j = 1}^m {{\lambda _j} = 1}$, and $0 \le {\lambda} _j < 1$ for any $j = 1, ..., m$.

It can be proved that the optimal solution of model (5.1) must be a weakly Pareto efficient solution of original MOP problem (4.3), and different weight vectors can lead individuals to diverse Pareto efficient solutions. Here, we use the uniform randomly (UR) method to initialize the weight vectors for subproblems (5.1). The process can be described as the following steps:

$Step\; 1.$ Randomly generate 5, 000 weight vectors to construct the set ${\mathit{\boldsymbol{\lambda}} _1}$, and initialize a vector set $\mathit{\boldsymbol{\lambda}}'$ that contains all dimensional unit weight vectors including (1, 0, ..., 0), (0, 1, ..., 0), ..., (0, 0, ..., 1);

$Step\; 2.$ Find the weight vector in ${\mathit{\boldsymbol{\lambda}}}_1$ which has the longest Euclidean distance with the vectors in ${\mathit{\boldsymbol{\lambda}}}$, and transfer it from ${\mathit{\boldsymbol{\lambda}}}_1$ to $\mathit{\boldsymbol{\lambda}}'$;

$Step\; 3.$ If the number of vectors in $\mathit{\boldsymbol{\lambda}}'$ reaches $N$ (population size), go to $Step\; 4$; Otherwise repeat $Step\; 2$;

$Step\; 4.$ Normalize the weight vectors in $\mathit{\boldsymbol{\lambda'}}$ by Eq (5.2) to obtain a series of initialized weight vectors $\mathit{\boldsymbol{\lambda}}_ 0$,

By the above method, the original MOP problem is decomposed to $N$ subproblems with $N$ weight vectors, denoted by $\mathit{\boldsymbol{W}}$ ($\mathit{\boldsymbol{W}} = (W(1), ..., W(N))$). After that, a VP containing all AAPs is randomly initialized as the initial individual for each subproblem, denoted by $VP_i$ ($i = 1, ..., N$), and its corresponding feasible solution is constructed by the CP. Then, the objective values of it are calculated. All the $N$ individuals make up the initial population.

5.2.2.   Neighborhood selection

The neighbors of individual $i$ are defined as follows:

where $\delta \in (0, 1)$ is the parameter that determines the size of the neighborhood.

In the exploitation search of each individual, we will select its neighbors as muting pool since their weight vectors are closed, i.e., the Euclidean distance is small, and their corresponding subproblems generally have similar optimal solutions, which means that the knowledge of them can be shared to each other.

5.2.3.   Evolutionary operation

For each individual, it corresponds to a subproblem of TCH approach with a weight vector, so our goal is to lead it to the optimal solution of its subproblem. In this paper, the following three evolution operators are proposed to search for better solutions for each individual (subproblem).

1) Local search operator

The first evolution operator is the local search (LS) operator, which is to generate a new solution by exchanging the order of two AAPs in the VP. As stated in Property 1 and Property 2, a feasible solution of UMDWTA can be generated by different VPs, so randomly changing the order of AAPs in a VP may leads to the same feasible solution. From the element of CP, we can see that any AAP ranking first in a VP will be marked as an AAAP, so selecting an UAAP and exchanging it with the first AAP in a VP can greatly generate a new feasible solution. Figure 2 describes the rationale of it.

We apply the LS operator to each UAAP in the current VP in turns, and the set of new generated VPs is called the neighborhood of the LS operator. Then, we find the best one among them (which has the optimal objective value for the subproblem) and replace the current VP with it. Note that the number of neighbors of the LS operator will rapidly grow with the number of AAPs in the VP, and the computing time consumed by the LS operator will be also increases rapidly. To avoid this case, a tenure parameter $t_e$ is set, and the LS operator will be released when executing $t_e$ times.

2) Crossover operator

Two individuals may be selected as parents to generate offspring. The common AAAPs shared by the two parents will remain and be moved to the head to generate offspring. This strategy follows the idea that "good" genes should be inherited. Figure 3 describes the rationale of this operator.

3) Mutation operator

To avoid individuals falling into locally optimal solutions,

the mutation operator is employed to increase the diversity of the population. The rationale of the mutation operator is similar to that of the LS operator, but the number of selected UAAPs, denoted as $m$, is randomly selected from the set $\{ 2, ..., {\rm size} ({S_{ \rm{UAAP}}})\}$. Then, the selected UAAPs will be moved to the head of the current VP in sequence. Figure 4 describes the rationale of this operation when $m = 3$.

Combining the above three evolution operators, the offspring production mechanism can be described as the following steps:

$Step\; 1.$ For the $i$-th individual in the current population, ${VP_i}$, the LS operator is employed to search for the optimal solution $VP^*$ in its neighborhood of the LS operator. If $VP^* {\ne VP_i}$, regard $VP^*$ as the offspring of ${VP_i}$, denoted by $\bar {VP}_i = VP^*$; otherwise $VP^* = {VP_i}$, which means there is no better solution in its neighborhood of the LS operator, and then go to $Step\; 2$.

$Step\; 2.$ Randomly select mating pool from the neighbors defined in Eq (5.3), then produce the offspring as follows:

where $cr({{VP_{i}}})$ and $mu({{VP_{i}}})$ are the individuals obtained by performing the crossover operator and the mutation operator on ${{VP_{i}}}$, respectively; $p_{sl}$ is the probability of performing crossover operator within $0 < p_{sl} < 1$.

In $Step\; 1$, the LS operator is employed to search for the better solution locally, and if it exists, treat it as the offspring. In $Step\; 2$, if no better solution can be found locally, the crossover operator and mutation operator will be employed to guide the individual to jump out of the local optimal solution and seek the global optimal solution.

5.2.4.   Non-dominant solutions storage

In addition to the evolutionary population, a collection called external population ($\mathit{\boldsymbol{EP}}$) is defined to store all non-dominant solutions and corresponding Pareto efficient solutions found in the iteration process. The size of $\mathit{\boldsymbol{EP}}$ is set to twice that of the evolutionary population, $2N$. The $\mathit{\boldsymbol{EP}}$ is also the output when the algorithm stops.

5.2.5.   Weight vectors adaptation

When using MOEA/D to solve complex MOP problems whose true Pareto fronts are discontinuous or not uniformly distributed, the outputs obtained may perform poorly in diversity and distribution uniformity if use unchanged uniformly distributed weights generated by UR method. To overcome this defect, this paper employs the weight vectors adaptation mechanism which uses the sparsity level among the population to adaptively adjust the weight vectors of subproblems. This mechanism can efficiently improve the diversity of the Pareto front obtained.

First, based on the vicinity distance ^[21], the sparsity level is defined as follows:

where $L_2^{NN_i^j}$ is the objective function Euclidean distance of $j$-th individual, ${VP_j}$, along with that of its $i$-th nearest individual in the population, ${\mathit{\boldsymbol{ P }}}$. The $m$ (the number of objective functions) closest Euclidean distances are used.

The sparsity level of each individual in the population is calculated among the population itself by Eq (5.5), and the individual with the lowest sparsity level will be removed. If 5$\%$ of the population (subproblems), $nus$, have not been removed yet, repeat the above process.

After that, calculate the sparsity level of each individual in $\mathit{\boldsymbol{EP}}$ among the population, and the individual with highest sparsity level in $\mathit{\boldsymbol{EP}}$, $VP_{sp}$, will be selected as the individual to be added to the current population. The new corresponding weight vector of the selected individual (subproblem), ${\mathit{\boldsymbol{\lambda}} ^{sp}}$, is generated as follows:

where $f_k^{sp}$ is $k$-th objective value of $VP_{sp}$. The process is repeated until $nus$ new individuals (subproblems) are add to the population. The pseudocode is shown in Algorithm 2.

5.2.6.   Algorithm explanation

Straightforwardly, the pseudocode of the improved MOEA/D with adaptive weights (MOEA/D-AW) is illustrated in Algorithm 3.

At the beginning, a set of weight vectors $\mathit{\boldsymbol{W}}$ generated by the UR method decompose the original MOP problem (4.3) to a series of subproblems, and then the neighborhood structure is initialized to $\mathit{\boldsymbol{B}}$. Then, initialize an individual for each subproblem in $\mathit{\boldsymbol{W}}$, and all individuals form a population ${\mathit{\boldsymbol{ P }}}$ (Line 1–4).

In each generation, the LS operator is first employed to search for the best solution locally for each individual (subproblem). If there is a better one, treat it as an offspring. Otherwise, a mating pool $\mathit{\boldsymbol{E}}$ is selected for the individual to produce an offspring by employing the crossover operator or mutation operator. Then, the offspring will update at most $nr$ individuals in $\mathit{\boldsymbol{E}}$ if the objective value of its corresponding subproblem can be improved. Meanwhile, the $\mathit{\boldsymbol{EP}}$ will be updated by combining the offspring with the current $\mathit{\boldsymbol{EP}}$ and deleting the solutions that are dominated by any other solution in the combination (Line 9–27).

After each generation, if the size of $\mathit{\boldsymbol{EP}}$ is larger than 2$N$, the sparsest individual in $\mathit{\boldsymbol{EP}}$ will be deleted (Line 28–30).

Finally, for every 5$\%$ of the evolutionary process, up to 90$\%$ of the total number of generations, the weight vectors of subproblems and population will be updated through weight vectors adaptation mechanism. The algorithm stops when the maximum number of generations is reached. The $\mathit{\boldsymbol{EP}}$ is the Pareto front and Pareto efficient solutions obtained by the algorithm (Line 31–36).

5.3.   Performance test

To verify the effectiveness of improved MOEA/D-AW in solving MOP problem, four typical test functions, ZDT1, ZDT2, ZDT3 and ZDT4 ^[42], with convex, non-convex and discontinuous Pareto fronts respectively are employed to test the performance of it in solving MOP problem firstly. Notice that above four test functions are continuous real-value functions, thus the evolution operator is replaced by original crossover and mutation operators in genetic algorithm ^[43].

The search space is set to ${[0, 1]^n}$, $n = 30$ for ZDT1, ZDT2, ZDT3, and $[0, 1]\times [-5, 5]^{n-1}$, $n = 10$ for ZDT4. The traditional MOEA/D ^[16] is employed as the comparison algorithm, and the inverse generational distance (IGD) and spacing metric (SP) ^[17] are selected as performance indexes to measure the convergence and distribution uniformity of Pareto fronts obtained by algorithms, respectively. The smaller the IGD is, the closer the Pareto front obtained by the algorithm is to the real Pareto front. The smaller the SP is, the better the distribution uniformity of the Pareto front obtained by the algorithm is. The control parameters of two algorithms are given in Table 2.

For each test function, the initial population is set to be the same, and two algorithms will run 30 times independently. Tables 3 and 4 show the performance indexes of MOEA/D-AW and MOEA/D on ZDT1, ZDT2, ZDT3 and ZDT4, respectively. It is shown that the average IGD and SP of 30 experiments obtained by MOEA/D-AW are smaller than those obtained by MOEA/D, which means that the Pareto fronts obtained by MOEA/D-AW are generally closer to the true Pareto fronts and more evenly distributed than that obtained by MOEA/D under the same parameter settings. In fact, the UR method initializes a set of uniformly distributed weight vectors, and the weight vectors adaption mechanism can adaptively adjust the weight vectors of subproblems according to the sparisity level of individuals, which can avoid the overcrowding of population in the iterative process and effectively improve the distribution uniformity of the Pareto front obtained. Undoubtedly, the MOEA/D-AW outperforms MOEA/D in solving the MOP problem, so it is reasonable to apply the improved MOEA/D-AW to solving the UMDWTA problem.

6.   Simulation case

In order to further verify the availability of the UMDWTA model and the applicability of improved MOEA/D-AW in solving the UMDWTA problem, a case is presented and solved by the algorithm in this Section.

6.1.   Relevant parameters of UMDWTA

In this case, the basic parameters in model (3.8) and model (4.3) are set as $W = 20$, $T = 10$, $S = 5$, $K = 5$, and other parameters are shown in Table 5.

The parameter $m_j$ $(j = 1, ..., T)$ is set to 1, 2 or 3 according to Eq (6.1), which follows the rules that the higher the target threat value is, the more weapons are allowed to be allocated to it.

6.2.   Experienmental parameters

The experimental parameters are shown in Table 6.

6.3.   Comparison algorithm

In order to verify the effectiveness of the algorithm, we use the improved MOEA/D without adaptive weights mechanism. The rest of the framework is the same as the improved MOEA/D-AW. Due to the large number of objective functions of the UMDWTA model, the Pareto fronts obtained by the two algorithms are not convenient to be visually represented in figure form. Here we calculate the SP ^[17] and C-metric ^[26] to give a further comparison study between improved MOEA/D-AW and improved MOEA/D on solving UMDWTA problem. The definition of C-metric is given as follows:

where $PF_1$ and $PF_2$ are the Pareto fronts obtained by two algorithms, respectively. The C-metric can reflect the quality of the Pareto front obtained by the two algorithms, i.e., if $C({PF_1}, {PF_2}) < C({PF_2}, {PF_1})$, the proportion of solutions in $PF_2$ dominated by $PF_1$ is smaller than that of solutions in $PF_1$ dominated by $PF_2$, so the quality of $PF_2$ is better than that of $PF_1$.

6.4.   Result analysis

The two algorithms will run 10 times for this case in MATLAB\_R2018b environment on a PC with Intel Core i5 CPU 2.30GHz and 12 GB internal memory. The average running time of improved MOEA/D-AW and improved MOEA/D is 21.4326 s and 19.8357 s, respectively.

The statistical results of C-metric and SP of Pareto fronts on the UMDWTA case are shown in Tables 7 and 8. Denote by $PF_1$ and $PF_2$ the Pareto front obtained by improved MOEA/D-AW and that obtained by improved MOEA/D, respectively.

There are two observations that can be derived from Tables 7 and 8. First, both improved MOEA/D-AW and improved MOEA/D can obtain the Pareto front with high quality in solving the case of UMDWTA problem, but the Pareto front obtained by improved MOEA/D-AW contains more non-dominated solution than that obtained by improved MOEA/D. Second, as to the comparison results in Table 8, the SP of Pareto front obtained by improved MOEA/D-AW is smaller than that of Pareto front obtained by improved MOEA/D averagely, i.e., the improved MOEA/D-AW outperforms the improved MOEA/D in distribution uniformity, which is consistent with former analysis that the adaptive weights mechanism can improve the distribution uniformity performance of output.

After that, we randomly pick one simulation result, and the Pareto front of $E\sigma$-UMDWTA model (4.3) and corresponding Pareto efficient solutions which are also called Pareto efficient allocation schemes are obtained. Part of the Pareto front obtained by improved MOEA/D-AW is shown in Table 9.

Table 9 shows that different objective functions of the $E\sigma$-UMDWTA model restrict each other, which is mainly reflected in two aspects. One is the mutual restriction between different evaluation indicators: The solution with smaller $E(-f_d)$ or smaller $E(-f_s)$ usually has larger $E(f_c)$ (No.1, No.2, No.10), i.e., allocation scheme with higher expected value of destruction of targets or higher expected value of surviving assets usually consumes more operation cost than other allocation schemes, which is also consistent with the actual combat situation. The second is the contradiction between the expected value and the standard deviation of the same evaluation indicator: The allocation scheme which performs better in the expected value of destruction of targets ($E(-f_d)$) or surviving assets ($E(-f_s)$) generally has greater standard deviation, i.e., the stability is not ideal.

It can be seen from Table 9 that the Pareto efficient allocation schemes of the $E\sigma$-UMDWTA model are contradictory, and there is no single allocation scheme that minimizes all the objective functions, i.e., there is no absolute optimal allocation scheme. If the decision makers hope to maximize the expected value of destruction of targets, i.e., to minimize $E(-f_d)$, they have to sacrifice the operation cost. Similarly, if they hope to maximize the expected value of surviving assets, it will also take more operation cost. This result verifies the availability of the model to some extent.

In conclusion, the improved MOEA/D-AW designed in this paper can effectively solve UMDWTA problem and obtained a series of Pareto efficient allocation schemes, and then different optimal allocation schemes can be defined under different combat missions.

There are mainly two methods for decision makers to select the optimal allocation scheme from the obtained Pareto efficient allocation schemes for different combat missions. The first method is to set different weights to each objective function according to specific requirement, which is also called the linear weighting method; the second method is to find the optimal allocation scheme of minimizing one objective when other objectives meet certain limits, which is also called the main target method.

For example, we assume that the decision makers pay more attention to the expected value of destruction of targets, so we given the weight vector $w = (0.5.0.1, 0.1, 0.1, 0.1, 0.1)$ to objective functions. Then, we can choose the optimal allocation scheme directly from the Pareto efficient allocation schemes obtained by algorithm directly and list it in Table 10. The objective value of this allocation scheme after normalization under the given weights is 58.3218. Commonly, the decision makers can adjust the weights according to the actual requirement and choose the corresponding optimal allocation scheme from the Pareto efficient allocation schemes.

7.   Conclusions

In this paper, the WTA problem was studied based on uncertainty theory originally, and an UMDWTA model was established. Then, based on the actual battlefield requirements, a solution method based on expected value and standard deviation principle was provided to convert the proposed unsolvable model to a deterministic one, which was a traditional MOP problem and had a set of Pareto efficient solutions instead of absolutely optimal solution. In order to solve the model efficiently, an improved MOEA/D with adaptive weights was designed, and the algorithm has been tested on four well-known MOP test functions, and the results showed that MOEA/D-AW outperforms MOEA/D greatly in both convergence and unform distribution. Finally, a case of UMDWTA was presented and solved by the designed algorithm, and a set of Pareto efficient solutions was obtained, which provided the decision maker a series of reference schemes. The results suggested that the objective values and the allocation schemes of UMDWTA model obtained by algorithm were consistent with the real situation, and the designed algorithm was an efficient multi-objective optimizer. The work in this paper can provide theoretical guidance and support for solving WTA problem in the indeterminate battlefield environment which is full of uncertainties.

Acknowledgements

This work is supported by the Research Fund of Fundamentals Department of Air Force Engineering University (JK2022204).

Conflict of interest

The authors disclosed no conflicts of interest in publishing this paper.