# Weapon Target Assignment Model for Small Unit Ground Combat Using Mixed Integer Nonlinear Program and Lagrangian Relaxation

Authors: Yeonseung Kim; Jaeyoung Cho; Sangwoo Han; Namsuk Cho
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The importance of the weapon target assignment problem is increasing as weapon systems are expected to become intelligent and unmanned on the future battlefield. This study discusses the weapon-target assignment in small-scale ground force combat, which has not been addressed much in other studies. In ground combat, target assignment means the start of an engagement, so assignment should be made taking into account not only the probability of killing an enemy but also the probability of surviving friendly forces. We introduce the WTAG (WTA ground) model that reflects these characteristics and describes methodologies for solving the problem including mixed integer nonlinear problems and Lagrangian relaxation approaches. Our finding is that the algorithm presented can efficiently provide the bounds of the problem in a practical size.


---

## Body

## 1. Introduction

In the future battlefield, the use of unmanned weapon systems will be further expanded. UAV (unmanned aerial vehicle) and drones are currently being used on the battlefield, and the use of unmanned weapons on the ground, called UGV (unmanned ground vehicle), will play an important role in the near future. Unmanned weapon system has to make various decisions autonomously, in the chaotic battle field, so the high level of situational judgement is required. Accordingly, the movement formation and optimal movement path to achieve the reconnaissance goal of the UAV cluster, stealth effect from detection, relay communication, tracking, and hitting targets are being studied recently. In particular, the target assignment for UAV has been in the spotlight in academic and military communities.On the other hand, UGV is being developed to assist infantry small combat, such as reconnaissance of dangerous areas, ammunition delivery, casualty transport, and close combat with the enemy. Here, in order to efficiently perform the close combat next to the maneuvering infantry, it is important which targets are allocated to UGV in a very short time. In addition, most recently, studies on counter-UAV, protecting high-value assets from attacks by unmanned weapon system, have also been conducted naturally with the development of UAV research. That is, if targets are assigned to attack UAV, a defense system has to be established in consideration of the assignments, from the defense point of view.Considering these points, the most basic but important decision is target assignment. Target assignment is simply a process of matching weapons with targets. This problem may be solved with a few simple rules such as “nearest target first” or “follow predefined priority of targets.” However, in an effort to find an exact answer to this problem, an optimization problem, called the WTA (weapon target assignment), has been studied since the 1950s in the military operations research field. WTA is finding the best weapon-target assignment solution for minimizing the total expected survival value of the targets.In this research, we introduce the WTA model for small-scale ground combat. The important characteristics that distinguish ground combat from other battles are as follows.First, a common but not necessary assumption to the general WTA problem is “fire and forget” which means that we are concentrating on the assignment before shooting but do not pay much attention to what happens after that. For this reason, weapons were mainly thought of as missiles and artillery such that they had no direct contact with the target. However, if the infantry is assumed, the assignment of weapons to targets will inevitably entail engagement. Most ground combat takes place at close range, and the fact that friendly forces can see and attack the enemy means that the enemy can also do the same. Second, in ground combat, target assignment can be made differently depending on the operational stage. For example, in security area operations, it is more important to protect one’s combat power, so the commander will minimize unnecessary engagement. On the other hand, in the target area operations, the commander will pour our available combat power for a decisive victory. Therefore, in different operational stages, it is necessary to construct the model for a different purpose. Lastly, in small unit ground combat, target assignment decision must be made faster. In the case of missile or artillery fire assignment, there is time to carefully plan for allocation. However, in ground battle, commanders have to judge the situation quickly and respond when they encounter unexpected environment.In this study, we propose the WTAG (weapon target assignment optimization model for small unit ground battle) that reflects the above considerations. The remainder of the paper is organized as follows. In Section2, we introduce previous WTA studies and explain the uniqueness of our research. In Section 3, we describe the WTAG problem and suggest solution methods. In Section 4, we provide the computational results. In the last section, the conclusion of the study and future directions are explained.


## 2. Literature Review

The origins of the WTA problem date back to the 1950s. Manne [1], Braford [2], and Day [3] developed and standardized the problem in the 1950s and 1960s. Since then, it has developed with the following several characteristics. First, it is NP-complete. Lloyd and Witsenhausen established the NP-completeness of the problem in 1986 [4]. So, as the instance size increases, the time of the algorithm for finding the optimal solution increases exponentially, and moreover, there is no optimal algorithm so far. Second, the WTA problem is a MINLP (mixed integer nonlinear programming). LP (linear programming) has an exact solution method like the simplex method, but MINLP does not. Thus, the relaxed method is usually used when solving the MINLP. Many researchers have suggested different formulations and devised exact algorithms to solve the problem. Relatively recently, there have been attempts by applying heuristic algorithms. Among them, an exact method had been researched for some of the following special cases.First of all, Manne [1] explained that it was possible to devise a LP formulation of the WTA problem under the assumption that all weapons are identical. And, denBroeder et al. [5] came up with an algorithm that provides an optimal solution to Manne’s problem. Day [3] researched the dimensionality and complexity of the problem. Day [3] considerably reduced the dimensions of the problem by decomposing the assignment problem into smaller targeting problems, and it provides information for solving larger targeting problems.Bracken and McCormick [6] handled two kinds of weapon allocation problems and suggested the LP approximation method. The first problem is to maximize target damage, and the second problem is to minimize cost with constraints that at least a specified target damage has to be needed. These traditional WTA problems, which is called static WTA problem, with the objective function of minimizing the threat or survival probability of the enemy were also suggested by Kolitz [7], Lemus and David [8], and Zu [9]. Chang et al. [10] suggested a nearly optimal algorithm for the large-scale WTA problem. Assuming that at least one weapon can be assigned to one target, a nearly optimal solution can be obtained. More recently, Ni et al. [11] devised the Lagrangian relaxation method and decomposed the Lagrangian function into smaller subproblems. And each problem can be solved more easily than the original problem. Andersen et al. [12] devised a specialized exact algorithm that is based on the piecewise linear approximation of the objective function which is computationally solvable.Meanwhile, dynamic WTA is a multistage problem, and it considers the weapons and targets to be attacked for subsequent stages. Soland et al. first suggested the dynamic version of the WTA problem with the assumption that there is a one kind of asset, and the assets are identical in each stage. The dynamic problem has been studied by Wacholder et al. [13], Christ et al. [14], and Hosein and Athans [15]. Especially, Hosein [16] devised the dynamic versions of the asset-based and target-based WTA. Hosein used the maximum marginal return algorithm and suboptimal algorithm, and it is efficient to obtain the bounds of the optimal solution. Silav et al. [17] developed biobjective model for dynamic WTA, and a linearization approach is applied to the objective function. And the result showed that the approach can be used for real-time constrained WTA problem and can assist decision-making in a relatively short time.Also, some of the heuristic method researches were suggested to solve the problem. Xin et al. [18] devised the Tabu search for finding a solution, Madni and Andrecut [19] used simulated annealing and threshold-accepting methods, and Xin et al. [20] suggested an efficient rule-based constructive heuristic method. Hong et al. [21] compared three heuristic algorithms, ACO (ant colony optimization), PSO (particle swarm optimization), and GA (genetic algorithm), and found out that ACO was outperforming in overall scenarios. And many researchers have tried to improve the original algorithms for better optimality and solving time. Hong et al. [21] proposed a population initialization algorithm to improve the ability of GA for solving WTA. When initializing the algorithm, the characteristics of the WTA are reflected, the dominant gene is inherited, and the search space is set wide so that a high-quality solution can be found efficiently. Li et al. [22] modified Pareto ACO, and in order to avoid defects in the original ACO, Li et al. tried to use an improved movement probability rule, a dynamic evaporation rate strategy, and a global updating rule of pheromone. Chang et al. [23] used ABC (artificial bee colony) algorithm for improving the timeliness of the dynamic WTA. Chang et al. progressed the original ABC algorithm, and the initial solution could be obtained more effectively. These heuristic methods may not be able to guarantee the optimal solution, but it is effective to find a local solution to the MINLP problem in a short time.So far, we reviewed representative WTA studies, and our research differs from other studies as follows. First, we propose a new mathematical model that has not been discussed before, reflecting the characteristics of ground combat described in Section 1. More specifically, our model also considers the survivability of weapons, which will become an important characteristic as battles between autonomous unmanned systems become more frequent in the future. Second, in terms of methodology for solving the problem, we have to admit that we borrowed the core idea of Ni et al. [11] that uses Lagrangian formulation and divides objective function into a few easily solvable subproblems. However, since the similar idea was applied to an entirely different problem, a distinct approach was needed in various aspects, such as the construction of Lagrangian functions and obtaining of solution for subproblems. Moreover, we described algorithms for finding Lagrangian multipliers with a comparative experiment, which was not covered in detail by Ni et al.


## 3. Methodology

### 3.1. The WTAG Problem

In this section, we define the WTAG problem and explain solving strategy. The notations and formulation of the problem are as shown in Table1.Table 1
Definitions and notations.
SetsW=0,1,2,…,mSet of friendly assets/units/weaponsT=0,1,2,…,nSet of enemy units/targetsParameterspijProbability of kill targetj by asset i∀i∈W,∀j∈TsijProbability of kill asseti by target j∀i∈W,∀j∈TViValue of asseti∀i∈WVjValue of targetj∀j∈TDecision variablexij=1,assetiisassignedtotargetj,0,o/w.∀i∈W,∀j∈TThen, our WTAG problem can be defined as follows:(1)min∑j∈TVj∏i∈W1−pijxij,(2)s.t.∑j∈Txij≤1∀i∈W,(3)∑i∈WVi∏j∈T1−sij∑i′∈Wxi′jxij≥β,(4)xij∈0,1∀i∈W,∀j∈T.The objective function is to minimize the expected target value, and inequality (2) means that all weapons can be assigned only once at a time. That is, weapons only can attack no more than one target. The inequality (3) is a single new constraint added to the original WTA problem. The left-hand side estimates the friendly assets’ remaining value after an engagement, and it is bounded by constant β in right-hand side. The constant β∈R is useful for implementing operation phases in military. For example, commanders in the security area may want to keep their combat power in reserve for the decisive operation. In this case, we may save the combat power of friendly forces by using large β. In contrast, in decisive operations, commanders will try to achieve their objective even if they consume almost all of their combat power. In this case, we may make an engagement actively using small β.There are two important assumptions in making the inequality (3). First, we determine that the enemy’s counter-attack is an immediate reaction. Figure 1(a) shows such a case. In this case, if W1 is assigned to T1, T1 also counter-attack W1. In contrast, Figure 1(b) shows nonimmediate reaction case. In this case, T1 can attack either W1 or W2, depending on the intention of the enemy. In this research, the enemy’s response is assumed to be (a).Figure 1
Enemy’s counter-attack cases. (a) Immediate reaction and (b) nonimmediate reaction.
(a)(b)Another important assumption is that each agent only can be assigned to no more than one enemy since an agent is considered as a weapon (e.g., unmanned ground vehicle, machine gunner) not a unit (e.g., squad, platoon). So, in Figure 2, the assignments of friendly assets are valid, but those of enemy targets are not. Since we assume that each enemy agent is also a weapon, not a unit, it can counterattack against only one friendly asset.Based on two assumptions, it is reasonable to divide parametersij by the number of assignments to target j in the sense that it is a kind of expected value. For instance, in Figure 2, the expected probability of kill asset W1 by target T1 can be estimated by sw1,t1/2 since T1 can either counter-attack to W1 or W2. In a similar way, probability of kill asset W2 by target T1 is sw2,t1/2. In this case, we assume that the probability that targets counter-attack to friendly assets follows a uniform distribution. However, because inequality (3) has the denominator that can have zero (i.e., none of the weapons are assigned to a certain target), it is not properly defined. We, hence, introduce a reformed formulation, WTAG_E (WTAG_Extended) which gives the same solution to WTAG but deals with the zero-divide issue using a surrogate modeling technique.(5)min∑j∈TVj∏i∈W1−pijxijs.t.∑j∈Txij≤1∀i∈W,(6)∑i∈Wxij≤M∗wj,∀j∈T,(7)wj≤∑i∈Wxij,∀j∈T,(8)tij≤sij+M∗wj,∀i∈W,∀j∈T,(9)tij≥sij−M∗wj,∀i∈W,∀j∈T,(10)tij≤sij∑i′∈Wxi′j+ϵ+M∗1−wj,∀i∈W,∀j∈T,(11)tij≥sij∑i′∈Wxi′j+ϵ−M∗1−wj,∀i∈W,∀j∈T,(12)∑i∈WVi∏j∈T1−tijxij≥β,xij∈0,1,∀i∈W,∀j∈T,wj∈0,1,∀j∈T.In this formulation,M is a big number, and ϵ is a small number. The decision variable tij∈R+ and wj∈0,1 is added in the formulation. The tij is a positive value that replaces the term sij/∑i′∈Wxi′j, and wj is a binary variable required to define tij. If ∑i′∈Wxi′j is zero, wj has to be zero by the inequality (8), and tij has to be sij by the inequality (9) and (10). At this moment, the inequality (11) and (12) does not force anything to the bound for tij. And if ∑i′∈Wxi′j≥1, wj has to be 1 by the inequality (7), and tij has to be sij/∑i′∈Wxi′j by the inequality (11) and (12). Similarly, the inequality (9) and (10) does not make meaningful bound for tij at this time.Unfortunately, both the WTAG and WTAG_E problem are NP-complete. Hence, it is almost not possible to solve the formulation efficiently. In the remainder of the section, we discuss an efficient algorithm for solving WTAG. As we mentioned before, we mainly refer to Ni et al. [11] for developing the idea.


### 3.2. Efficient Algorithm

The key idea of the algorithm is to reform the WTAG problem by using the Lagrangian relaxation method. And the relaxed objective function is divided into three subproblems.


#### 3.2.1. Lagrangian Relaxation

Before we construct the Lagrangian relaxation of the WTAG, we first transpose a term that has power using a natural log.(13)yj=∏i∈W1−pijxij∀j∈T,→lnyj=∑i∈Wln1−pij⋅xij∀j∈T.As the natural log is a nondecreasing monotonic function, it does not spoil the solution. Here,pij is the probability of kill target j by asset i, so yj means the survival probability of the target j. In this way, the survival probability of the friendly assets, zi, can be transposed and rearranged as follows:(14)zi=∏j∈T1−sij∑i′∈Wxi′jxij,∀i∈W,(15)→lnzi=∑j∈Tln1−sij∑i′∈Wxi′j⋅xij,∀i∈W,(16)⇔lnzi=∑j∈Tln1−siju⋅xij,∀i∈W,1≤u≤W,We substitute the denominator term in (19) using surrogate constant u where 1≤u≤W. Denominator term is a number of assignments that allocated to target j. As we discussed, such number can be zero. However, in this equation, as a term xij is multiplied, if there is no assignment, such term automatically has zero value. Hence, the bound of u is proper. The equations related to the transpositions are added to the WTAG formulation.(17)min∑j∈TVj⋅yj,s.t.∑j∈Txij≤1∀i∈W,(18)∑i∈WVi⋅zi≥β,(19)lnyj=∑i∈Wln1−pij⋅xij,∀j∈T,(20)lnzi=∑j∈Tln1−siju⋅xij∀i∈W,1≤u≤W,(21)xij∈0,1,∀i∈W,∀j∈T.In this formulation, theyj and zi are also decision variables, as well as xij. Variables yj and zi are nonnegative as it represents a probability. Now, the WTAG Lagrangian function is as follows:(22)minLx,y,z,λ=∑j∈TVj⋅yj−λβ∑i∈WVi⋅zi−β+∑j∈Tλjy∑i∈Wln1−pij⋅xij−lnyj+∑i∈Wλiz∑j∈Tln1−siju⋅xij−lnzi,(23)s.t.∑j∈Txij≤1∀i∈W,(24)λβ≥0,λiz,λjy−free∀i∈W,∀j∈T.The equations (23)–(25) are moved to the objective function with corresponding multipliers. Here, λβ is one dimensional vector, and λjy, λiz is T, W dimensional vector, respectively. So, integrated multiplier vector λ is a 1+T+W dimensional vector. The important point here, theoretically, is that if this problem is unbounded, it means that there are implicit restrictions on the dual variables, λ vector. It has to be applied when solving dual problem using subgradient method. Even we are given multiplier vector, it is still a mixed integer nonlinear problem. Accordingly, in the next section, we propose a solving method by decomposing the objective function into three subproblems.


#### 3.2.2. Subproblems

The Lagrangian function can be separated according to the decision variables,xij, yj, zi, and we define each as L1x, L2y, L3z.(25)L1x=min∑j∈Tλjy∑i∈Wln1−pij⋅xij+∑i∈Wλiz∑j∈Tln1−siju⋅xij+λβ⋅β,⇔∑j∈T∑i∈Wλjy⋅ln1−pij+λjz⋅ln1−siju⋅xij+λβ⋅β,(26)L2y=min∑j∈TVj⋅yj−λjy⋅lnyj,(27)L3z=min−∑i∈Wλβ⋅Vi⋅zi+λiz⋅lnzi.We now describe how to solve each subproblem.


#### 3.2.3. Solving SubproblemL1x

AsL1x term includes variable x, we consider constraint of xij.(28)L1x=min∑j∈T∑i∈Wλjyln1−pij+λizln1−sijuxij+λββ,s.t,∑j∈Txij≤1∀i∈W.Given multiplier, the objective function is a linear function ofxij, and the constraint matrix is a totally unimodular (TU). Therefore, we solve L1x as linear programming.


#### 3.2.4. Solving SubproblemL2y

L2y formulation, (28), can be decomposed into T one-dimension discrete problems. And we find a solution for each subproblem.(29)SL2jyj=minVj⋅yj−λjy⋅lnyj,∀j∈T.As discussed,yj is the survival probability of the target j. If the asset i that has the largest pij attacks target j, the survival probability of the target j will have smallest value, γj. In contrast, if the asset i that has the smallest pij attacks target j, the survival probability of the target j will have largest value, αj. We define bound of the yj as follows:(30)αj=1,∀j∈T(31)γj=1−max1≤i≤WpijW,∀j∈T(32)γj≤yj≤αj,∀j∈TSince there may not be any assets that attack targetj, the maximum value of the survival probability of the target j is 1. The solution yj can be obtained depending on the shape of the objective function graph within the bound, γj≤yj≤αj. The shape of the graph can be different from the value of λjy because the Vj and yj are positive. If λjy<0, the SL2jyj graph is monotonic increasing as yj increases, as seen in Figure 3(a). In this case, the minimum value within the range of y is obtained when yj=γj. If λjy=0, minimum value is obtained when yj=γj. Lastly, if λjy>0, SL2jyj is convex, and the minimum value is obtained at the critical point where the gradient is 0, as seen in Figure 3(b). And it is important to check whether the critical point is in the range of yj. Accordingly, the solution can be obtained by the conditions as follows:(33)yj∗=αj,αj≤λjyVj,λjVj,γj≤λjyVj≤αj,γj,λjyVj≤γj.Figure 2
Assignments of the enemy target.Figure 3
Dynamics ofSL2jyj with different λ values. (a) λjy<0, (b) λjy>0.
(a)(b)


#### 3.2.5. Solving SubproblemL3z

TheL3z, (28), can be decomposed into W one-dimension discrete problems.(34)SL3izi=min−λβ⋅Vi⋅zi−λiz⋅lnzi,∀i∈W.Aszi represents the survival probability of asset i, we define the bound of the zi as follows:(35)αi=1,∀i∈W,(36)γi=1−max1≤j≤TsijuT,1≤u≤W,∀i∈W,(37)γi≤zi≤αi,∀i∈W.The solution to subproblem (33) can be obtained in a similar way of L2y by observing the dynamics of L3z with respect to z.As shown in Figure4(a), when λiz<0, we have a concave function. However, different to the case of subproblem L2y, we solve minimization problem, and the critical point is no more optimal. And if λiz≥0, it shows monotonic nonincreasing, as seen in Figure 4(b). Therefore, we define the analytical solution of both cases as follows:(38)zi∗∈argminzL3αi,L3γi.However, γi in (35) has a variable u which is indeed u=∑i′∈Wxi′j. And this means that subproblem (33) is not separated from a variable x. Our heuristic remedy is to use u=W. Although this setting obviously will change the solution, it at least provides a feasible solution according to the following logic. Looking back on the previous discussion, zi is survival probability of weapon i. If u=W, zi has the greatest value. Let’s define such zi to be zmax. And we substitute zmax into (21), then following is satisfied.(39)∑i∈WVizumax≥∑i∈WVizopt≥β,where zopt is the solution where we set u=∑i′∈Wxi′jFigure 4
Dynamics ofSL3izi with different λ values. (a) λiz<0, (b) λiz>0.
(a)(b)


#### 3.2.6. Line Search Method

Finding exact multipliers in Lagrangian relaxation is another optimization problem, and we introduce two algorithms for solving the Lagrangian dual problem. The first approach is to find a localλ solution using iterative first-order line search, as described in Algorithm 1.Algorithm 1: Line search method: WTAG_LS
(1)
Set iteration numberk←1, initial λ1←1,1,…,1,1(2)
Solve 3 subproblems (solution:x^,y^,z^ )(3)
Computepk (step direction)            ⊳pk=∇λk(4)
Computeαk (step size)              ⊳αk=γ1/k+1λkmax/pkmax, 0≤γ≤1(5)
ifpk>ϵthen                ⊳ϵ = tolerance(6)
λk+1←λk+αkpk(7)
Go line 2 andk←k+1(8)
else ifpk≤ϵthen(9)
STOP(10)
end ifIn algorithm 1,λk=λjy,λiz,λβ,∀i∈W,∀j∈T. Through numerical experiments, we found that the choice of step length αk had the greatest effect on the performance of the algorithm. We empirically decided to adopt the ideas of Fan et al. [24] and Feng and Liao [25].(40)αk=γ1k+1λkmaxξkmax,0≤γ≤1,where λkmax is the maximum of multiplier for coupling constraint, and ξkmax is a maximum of subgradient for coupling constraint. In this algorithm, as described in line 8, it terminates the search if the gradient norm (i.e., search direction) is less than epsilon, 0.2 in this study. Moreover, as a result of experiments, we empirically suggest using one vector for the initial choice of λk.


#### 3.2.7. Subgradient Method

Another method for solving the Lagrangian dual problem is the subgradient method, and the pseudocode of the method for our problem is shown in Algorithm 2.In algorithm 2,λk=λjy,λiz,λβ,∀i∈W,∀j∈T at iteration step k. We define two initial lambda vectors nearly symmetric about the origin for the purpose of preventing the LP problem in line 5 is unbounded. The specific numbers in the initial lambda are determined based on parameters of instance which is described in Table 2. As a result of numerous experiments, we suggest to use minVj×max1−pij as the initial choice of lambda and −minVj×max1−pij for choice of second iterations. As explained in Section 3.2.1, since λβ≥0 has to be satisfied, we set the last element of the initial λ2 to zero [26]. As described in line 6, the algorithm stops if one of the two conditions is satisfied. One is to check if the value of the gradient norm (i.e., sk) gets closer to 0, and the second is to see if the L0 value doesn’t get any better for a certain period of iterations (i.e., L0k+1−L0∗/L0k−L0∗<ϵ where L0k is solution of k th iteration of LP problem in line 5 and L0∗ is best solution so far. In this study, ϵ=0.1) [27]. Our conjecture for both WTAG_LS and WTAG_SG is that first-order information (i.e., ∇λk) does not provide good information since we solve 3 subproblems in discrete manner.Algorithm 2: Subgradient Method: WTAG_SG
(1)
Set iteration numberk←1(2)
initialλ1←20,20,…,20,20, λ2←−20,−20,…,−20,0(3)
Solve 3 sub-problems (solution:x^,y^,z^ )(4)
Computesk (subgradient)               ⊳sk=∇λk(5)
Solve LP problem ofλ (solution: λ^ )max,L0,L0≤Lλk+skTλk+1−λk,(6)
ifsk≤ϵ or L0^ has no improvement then       ⊳ϵ = tolerance(7)
STOP(8)
else(9)
λk+1←λ^(10)
Go line 3,k←k+1, Add inequality to LP problem in line 5.(11)
end if


## 3.1. The WTAG Problem

In this section, we define the WTAG problem and explain solving strategy. The notations and formulation of the problem are as shown in Table1.Table 1
Definitions and notations.
SetsW=0,1,2,…,mSet of friendly assets/units/weaponsT=0,1,2,…,nSet of enemy units/targetsParameterspijProbability of kill targetj by asset i∀i∈W,∀j∈TsijProbability of kill asseti by target j∀i∈W,∀j∈TViValue of asseti∀i∈WVjValue of targetj∀j∈TDecision variablexij=1,assetiisassignedtotargetj,0,o/w.∀i∈W,∀j∈TThen, our WTAG problem can be defined as follows:(1)min∑j∈TVj∏i∈W1−pijxij,(2)s.t.∑j∈Txij≤1∀i∈W,(3)∑i∈WVi∏j∈T1−sij∑i′∈Wxi′jxij≥β,(4)xij∈0,1∀i∈W,∀j∈T.The objective function is to minimize the expected target value, and inequality (2) means that all weapons can be assigned only once at a time. That is, weapons only can attack no more than one target. The inequality (3) is a single new constraint added to the original WTA problem. The left-hand side estimates the friendly assets’ remaining value after an engagement, and it is bounded by constant β in right-hand side. The constant β∈R is useful for implementing operation phases in military. For example, commanders in the security area may want to keep their combat power in reserve for the decisive operation. In this case, we may save the combat power of friendly forces by using large β. In contrast, in decisive operations, commanders will try to achieve their objective even if they consume almost all of their combat power. In this case, we may make an engagement actively using small β.There are two important assumptions in making the inequality (3). First, we determine that the enemy’s counter-attack is an immediate reaction. Figure 1(a) shows such a case. In this case, if W1 is assigned to T1, T1 also counter-attack W1. In contrast, Figure 1(b) shows nonimmediate reaction case. In this case, T1 can attack either W1 or W2, depending on the intention of the enemy. In this research, the enemy’s response is assumed to be (a).Figure 1
Enemy’s counter-attack cases. (a) Immediate reaction and (b) nonimmediate reaction.
(a)(b)Another important assumption is that each agent only can be assigned to no more than one enemy since an agent is considered as a weapon (e.g., unmanned ground vehicle, machine gunner) not a unit (e.g., squad, platoon). So, in Figure 2, the assignments of friendly assets are valid, but those of enemy targets are not. Since we assume that each enemy agent is also a weapon, not a unit, it can counterattack against only one friendly asset.Based on two assumptions, it is reasonable to divide parametersij by the number of assignments to target j in the sense that it is a kind of expected value. For instance, in Figure 2, the expected probability of kill asset W1 by target T1 can be estimated by sw1,t1/2 since T1 can either counter-attack to W1 or W2. In a similar way, probability of kill asset W2 by target T1 is sw2,t1/2. In this case, we assume that the probability that targets counter-attack to friendly assets follows a uniform distribution. However, because inequality (3) has the denominator that can have zero (i.e., none of the weapons are assigned to a certain target), it is not properly defined. We, hence, introduce a reformed formulation, WTAG_E (WTAG_Extended) which gives the same solution to WTAG but deals with the zero-divide issue using a surrogate modeling technique.(5)min∑j∈TVj∏i∈W1−pijxijs.t.∑j∈Txij≤1∀i∈W,(6)∑i∈Wxij≤M∗wj,∀j∈T,(7)wj≤∑i∈Wxij,∀j∈T,(8)tij≤sij+M∗wj,∀i∈W,∀j∈T,(9)tij≥sij−M∗wj,∀i∈W,∀j∈T,(10)tij≤sij∑i′∈Wxi′j+ϵ+M∗1−wj,∀i∈W,∀j∈T,(11)tij≥sij∑i′∈Wxi′j+ϵ−M∗1−wj,∀i∈W,∀j∈T,(12)∑i∈WVi∏j∈T1−tijxij≥β,xij∈0,1,∀i∈W,∀j∈T,wj∈0,1,∀j∈T.In this formulation,M is a big number, and ϵ is a small number. The decision variable tij∈R+ and wj∈0,1 is added in the formulation. The tij is a positive value that replaces the term sij/∑i′∈Wxi′j, and wj is a binary variable required to define tij. If ∑i′∈Wxi′j is zero, wj has to be zero by the inequality (8), and tij has to be sij by the inequality (9) and (10). At this moment, the inequality (11) and (12) does not force anything to the bound for tij. And if ∑i′∈Wxi′j≥1, wj has to be 1 by the inequality (7), and tij has to be sij/∑i′∈Wxi′j by the inequality (11) and (12). Similarly, the inequality (9) and (10) does not make meaningful bound for tij at this time.Unfortunately, both the WTAG and WTAG_E problem are NP-complete. Hence, it is almost not possible to solve the formulation efficiently. In the remainder of the section, we discuss an efficient algorithm for solving WTAG. As we mentioned before, we mainly refer to Ni et al. [11] for developing the idea.


## 3.2. Efficient Algorithm

The key idea of the algorithm is to reform the WTAG problem by using the Lagrangian relaxation method. And the relaxed objective function is divided into three subproblems.


### 3.2.1. Lagrangian Relaxation

Before we construct the Lagrangian relaxation of the WTAG, we first transpose a term that has power using a natural log.(13)yj=∏i∈W1−pijxij∀j∈T,→lnyj=∑i∈Wln1−pij⋅xij∀j∈T.As the natural log is a nondecreasing monotonic function, it does not spoil the solution. Here,pij is the probability of kill target j by asset i, so yj means the survival probability of the target j. In this way, the survival probability of the friendly assets, zi, can be transposed and rearranged as follows:(14)zi=∏j∈T1−sij∑i′∈Wxi′jxij,∀i∈W,(15)→lnzi=∑j∈Tln1−sij∑i′∈Wxi′j⋅xij,∀i∈W,(16)⇔lnzi=∑j∈Tln1−siju⋅xij,∀i∈W,1≤u≤W,We substitute the denominator term in (19) using surrogate constant u where 1≤u≤W. Denominator term is a number of assignments that allocated to target j. As we discussed, such number can be zero. However, in this equation, as a term xij is multiplied, if there is no assignment, such term automatically has zero value. Hence, the bound of u is proper. The equations related to the transpositions are added to the WTAG formulation.(17)min∑j∈TVj⋅yj,s.t.∑j∈Txij≤1∀i∈W,(18)∑i∈WVi⋅zi≥β,(19)lnyj=∑i∈Wln1−pij⋅xij,∀j∈T,(20)lnzi=∑j∈Tln1−siju⋅xij∀i∈W,1≤u≤W,(21)xij∈0,1,∀i∈W,∀j∈T.In this formulation, theyj and zi are also decision variables, as well as xij. Variables yj and zi are nonnegative as it represents a probability. Now, the WTAG Lagrangian function is as follows:(22)minLx,y,z,λ=∑j∈TVj⋅yj−λβ∑i∈WVi⋅zi−β+∑j∈Tλjy∑i∈Wln1−pij⋅xij−lnyj+∑i∈Wλiz∑j∈Tln1−siju⋅xij−lnzi,(23)s.t.∑j∈Txij≤1∀i∈W,(24)λβ≥0,λiz,λjy−free∀i∈W,∀j∈T.The equations (23)–(25) are moved to the objective function with corresponding multipliers. Here, λβ is one dimensional vector, and λjy, λiz is T, W dimensional vector, respectively. So, integrated multiplier vector λ is a 1+T+W dimensional vector. The important point here, theoretically, is that if this problem is unbounded, it means that there are implicit restrictions on the dual variables, λ vector. It has to be applied when solving dual problem using subgradient method. Even we are given multiplier vector, it is still a mixed integer nonlinear problem. Accordingly, in the next section, we propose a solving method by decomposing the objective function into three subproblems.


### 3.2.2. Subproblems

The Lagrangian function can be separated according to the decision variables,xij, yj, zi, and we define each as L1x, L2y, L3z.(25)L1x=min∑j∈Tλjy∑i∈Wln1−pij⋅xij+∑i∈Wλiz∑j∈Tln1−siju⋅xij+λβ⋅β,⇔∑j∈T∑i∈Wλjy⋅ln1−pij+λjz⋅ln1−siju⋅xij+λβ⋅β,(26)L2y=min∑j∈TVj⋅yj−λjy⋅lnyj,(27)L3z=min−∑i∈Wλβ⋅Vi⋅zi+λiz⋅lnzi.We now describe how to solve each subproblem.


### 3.2.3. Solving SubproblemL1x

AsL1x term includes variable x, we consider constraint of xij.(28)L1x=min∑j∈T∑i∈Wλjyln1−pij+λizln1−sijuxij+λββ,s.t,∑j∈Txij≤1∀i∈W.Given multiplier, the objective function is a linear function ofxij, and the constraint matrix is a totally unimodular (TU). Therefore, we solve L1x as linear programming.


### 3.2.4. Solving SubproblemL2y

L2y formulation, (28), can be decomposed into T one-dimension discrete problems. And we find a solution for each subproblem.(29)SL2jyj=minVj⋅yj−λjy⋅lnyj,∀j∈T.As discussed,yj is the survival probability of the target j. If the asset i that has the largest pij attacks target j, the survival probability of the target j will have smallest value, γj. In contrast, if the asset i that has the smallest pij attacks target j, the survival probability of the target j will have largest value, αj. We define bound of the yj as follows:(30)αj=1,∀j∈T(31)γj=1−max1≤i≤WpijW,∀j∈T(32)γj≤yj≤αj,∀j∈TSince there may not be any assets that attack targetj, the maximum value of the survival probability of the target j is 1. The solution yj can be obtained depending on the shape of the objective function graph within the bound, γj≤yj≤αj. The shape of the graph can be different from the value of λjy because the Vj and yj are positive. If λjy<0, the SL2jyj graph is monotonic increasing as yj increases, as seen in Figure 3(a). In this case, the minimum value within the range of y is obtained when yj=γj. If λjy=0, minimum value is obtained when yj=γj. Lastly, if λjy>0, SL2jyj is convex, and the minimum value is obtained at the critical point where the gradient is 0, as seen in Figure 3(b). And it is important to check whether the critical point is in the range of yj. Accordingly, the solution can be obtained by the conditions as follows:(33)yj∗=αj,αj≤λjyVj,λjVj,γj≤λjyVj≤αj,γj,λjyVj≤γj.Figure 2
Assignments of the enemy target.Figure 3
Dynamics ofSL2jyj with different λ values. (a) λjy<0, (b) λjy>0.
(a)(b)


### 3.2.5. Solving SubproblemL3z

TheL3z, (28), can be decomposed into W one-dimension discrete problems.(34)SL3izi=min−λβ⋅Vi⋅zi−λiz⋅lnzi,∀i∈W.Aszi represents the survival probability of asset i, we define the bound of the zi as follows:(35)αi=1,∀i∈W,(36)γi=1−max1≤j≤TsijuT,1≤u≤W,∀i∈W,(37)γi≤zi≤αi,∀i∈W.The solution to subproblem (33) can be obtained in a similar way of L2y by observing the dynamics of L3z with respect to z.As shown in Figure4(a), when λiz<0, we have a concave function. However, different to the case of subproblem L2y, we solve minimization problem, and the critical point is no more optimal. And if λiz≥0, it shows monotonic nonincreasing, as seen in Figure 4(b). Therefore, we define the analytical solution of both cases as follows:(38)zi∗∈argminzL3αi,L3γi.However, γi in (35) has a variable u which is indeed u=∑i′∈Wxi′j. And this means that subproblem (33) is not separated from a variable x. Our heuristic remedy is to use u=W. Although this setting obviously will change the solution, it at least provides a feasible solution according to the following logic. Looking back on the previous discussion, zi is survival probability of weapon i. If u=W, zi has the greatest value. Let’s define such zi to be zmax. And we substitute zmax into (21), then following is satisfied.(39)∑i∈WVizumax≥∑i∈WVizopt≥β,where zopt is the solution where we set u=∑i′∈Wxi′jFigure 4
Dynamics ofSL3izi with different λ values. (a) λiz<0, (b) λiz>0.
(a)(b)


### 3.2.6. Line Search Method

Finding exact multipliers in Lagrangian relaxation is another optimization problem, and we introduce two algorithms for solving the Lagrangian dual problem. The first approach is to find a localλ solution using iterative first-order line search, as described in Algorithm 1.Algorithm 1: Line search method: WTAG_LS
(1)
Set iteration numberk←1, initial λ1←1,1,…,1,1(2)
Solve 3 subproblems (solution:x^,y^,z^ )(3)
Computepk (step direction)            ⊳pk=∇λk(4)
Computeαk (step size)              ⊳αk=γ1/k+1λkmax/pkmax, 0≤γ≤1(5)
ifpk>ϵthen                ⊳ϵ = tolerance(6)
λk+1←λk+αkpk(7)
Go line 2 andk←k+1(8)
else ifpk≤ϵthen(9)
STOP(10)
end ifIn algorithm 1,λk=λjy,λiz,λβ,∀i∈W,∀j∈T. Through numerical experiments, we found that the choice of step length αk had the greatest effect on the performance of the algorithm. We empirically decided to adopt the ideas of Fan et al. [24] and Feng and Liao [25].(40)αk=γ1k+1λkmaxξkmax,0≤γ≤1,where λkmax is the maximum of multiplier for coupling constraint, and ξkmax is a maximum of subgradient for coupling constraint. In this algorithm, as described in line 8, it terminates the search if the gradient norm (i.e., search direction) is less than epsilon, 0.2 in this study. Moreover, as a result of experiments, we empirically suggest using one vector for the initial choice of λk.


### 3.2.7. Subgradient Method

Another method for solving the Lagrangian dual problem is the subgradient method, and the pseudocode of the method for our problem is shown in Algorithm 2.In algorithm 2,λk=λjy,λiz,λβ,∀i∈W,∀j∈T at iteration step k. We define two initial lambda vectors nearly symmetric about the origin for the purpose of preventing the LP problem in line 5 is unbounded. The specific numbers in the initial lambda are determined based on parameters of instance which is described in Table 2. As a result of numerous experiments, we suggest to use minVj×max1−pij as the initial choice of lambda and −minVj×max1−pij for choice of second iterations. As explained in Section 3.2.1, since λβ≥0 has to be satisfied, we set the last element of the initial λ2 to zero [26]. As described in line 6, the algorithm stops if one of the two conditions is satisfied. One is to check if the value of the gradient norm (i.e., sk) gets closer to 0, and the second is to see if the L0 value doesn’t get any better for a certain period of iterations (i.e., L0k+1−L0∗/L0k−L0∗<ϵ where L0k is solution of k th iteration of LP problem in line 5 and L0∗ is best solution so far. In this study, ϵ=0.1) [27]. Our conjecture for both WTAG_LS and WTAG_SG is that first-order information (i.e., ∇λk) does not provide good information since we solve 3 subproblems in discrete manner.Algorithm 2: Subgradient Method: WTAG_SG
(1)
Set iteration numberk←1(2)
initialλ1←20,20,…,20,20, λ2←−20,−20,…,−20,0(3)
Solve 3 sub-problems (solution:x^,y^,z^ )(4)
Computesk (subgradient)               ⊳sk=∇λk(5)
Solve LP problem ofλ (solution: λ^ )max,L0,L0≤Lλk+skTλk+1−λk,(6)
ifsk≤ϵ or L0^ has no improvement then       ⊳ϵ = tolerance(7)
STOP(8)
else(9)
λk+1←λ^(10)
Go line 3,k←k+1, Add inequality to LP problem in line 5.(11)
end if


## 3.2.1. Lagrangian Relaxation

Before we construct the Lagrangian relaxation of the WTAG, we first transpose a term that has power using a natural log.(13)yj=∏i∈W1−pijxij∀j∈T,→lnyj=∑i∈Wln1−pij⋅xij∀j∈T.As the natural log is a nondecreasing monotonic function, it does not spoil the solution. Here,pij is the probability of kill target j by asset i, so yj means the survival probability of the target j. In this way, the survival probability of the friendly assets, zi, can be transposed and rearranged as follows:(14)zi=∏j∈T1−sij∑i′∈Wxi′jxij,∀i∈W,(15)→lnzi=∑j∈Tln1−sij∑i′∈Wxi′j⋅xij,∀i∈W,(16)⇔lnzi=∑j∈Tln1−siju⋅xij,∀i∈W,1≤u≤W,We substitute the denominator term in (19) using surrogate constant u where 1≤u≤W. Denominator term is a number of assignments that allocated to target j. As we discussed, such number can be zero. However, in this equation, as a term xij is multiplied, if there is no assignment, such term automatically has zero value. Hence, the bound of u is proper. The equations related to the transpositions are added to the WTAG formulation.(17)min∑j∈TVj⋅yj,s.t.∑j∈Txij≤1∀i∈W,(18)∑i∈WVi⋅zi≥β,(19)lnyj=∑i∈Wln1−pij⋅xij,∀j∈T,(20)lnzi=∑j∈Tln1−siju⋅xij∀i∈W,1≤u≤W,(21)xij∈0,1,∀i∈W,∀j∈T.In this formulation, theyj and zi are also decision variables, as well as xij. Variables yj and zi are nonnegative as it represents a probability. Now, the WTAG Lagrangian function is as follows:(22)minLx,y,z,λ=∑j∈TVj⋅yj−λβ∑i∈WVi⋅zi−β+∑j∈Tλjy∑i∈Wln1−pij⋅xij−lnyj+∑i∈Wλiz∑j∈Tln1−siju⋅xij−lnzi,(23)s.t.∑j∈Txij≤1∀i∈W,(24)λβ≥0,λiz,λjy−free∀i∈W,∀j∈T.The equations (23)–(25) are moved to the objective function with corresponding multipliers. Here, λβ is one dimensional vector, and λjy, λiz is T, W dimensional vector, respectively. So, integrated multiplier vector λ is a 1+T+W dimensional vector. The important point here, theoretically, is that if this problem is unbounded, it means that there are implicit restrictions on the dual variables, λ vector. It has to be applied when solving dual problem using subgradient method. Even we are given multiplier vector, it is still a mixed integer nonlinear problem. Accordingly, in the next section, we propose a solving method by decomposing the objective function into three subproblems.


## 3.2.2. Subproblems

The Lagrangian function can be separated according to the decision variables,xij, yj, zi, and we define each as L1x, L2y, L3z.(25)L1x=min∑j∈Tλjy∑i∈Wln1−pij⋅xij+∑i∈Wλiz∑j∈Tln1−siju⋅xij+λβ⋅β,⇔∑j∈T∑i∈Wλjy⋅ln1−pij+λjz⋅ln1−siju⋅xij+λβ⋅β,(26)L2y=min∑j∈TVj⋅yj−λjy⋅lnyj,(27)L3z=min−∑i∈Wλβ⋅Vi⋅zi+λiz⋅lnzi.We now describe how to solve each subproblem.


## 3.2.3. Solving SubproblemL1x

AsL1x term includes variable x, we consider constraint of xij.(28)L1x=min∑j∈T∑i∈Wλjyln1−pij+λizln1−sijuxij+λββ,s.t,∑j∈Txij≤1∀i∈W.Given multiplier, the objective function is a linear function ofxij, and the constraint matrix is a totally unimodular (TU). Therefore, we solve L1x as linear programming.


## 3.2.4. Solving SubproblemL2y

L2y formulation, (28), can be decomposed into T one-dimension discrete problems. And we find a solution for each subproblem.(29)SL2jyj=minVj⋅yj−λjy⋅lnyj,∀j∈T.As discussed,yj is the survival probability of the target j. If the asset i that has the largest pij attacks target j, the survival probability of the target j will have smallest value, γj. In contrast, if the asset i that has the smallest pij attacks target j, the survival probability of the target j will have largest value, αj. We define bound of the yj as follows:(30)αj=1,∀j∈T(31)γj=1−max1≤i≤WpijW,∀j∈T(32)γj≤yj≤αj,∀j∈TSince there may not be any assets that attack targetj, the maximum value of the survival probability of the target j is 1. The solution yj can be obtained depending on the shape of the objective function graph within the bound, γj≤yj≤αj. The shape of the graph can be different from the value of λjy because the Vj and yj are positive. If λjy<0, the SL2jyj graph is monotonic increasing as yj increases, as seen in Figure 3(a). In this case, the minimum value within the range of y is obtained when yj=γj. If λjy=0, minimum value is obtained when yj=γj. Lastly, if λjy>0, SL2jyj is convex, and the minimum value is obtained at the critical point where the gradient is 0, as seen in Figure 3(b). And it is important to check whether the critical point is in the range of yj. Accordingly, the solution can be obtained by the conditions as follows:(33)yj∗=αj,αj≤λjyVj,λjVj,γj≤λjyVj≤αj,γj,λjyVj≤γj.Figure 2
Assignments of the enemy target.Figure 3
Dynamics ofSL2jyj with different λ values. (a) λjy<0, (b) λjy>0.
(a)(b)


## 3.2.5. Solving SubproblemL3z

TheL3z, (28), can be decomposed into W one-dimension discrete problems.(34)SL3izi=min−λβ⋅Vi⋅zi−λiz⋅lnzi,∀i∈W.Aszi represents the survival probability of asset i, we define the bound of the zi as follows:(35)αi=1,∀i∈W,(36)γi=1−max1≤j≤TsijuT,1≤u≤W,∀i∈W,(37)γi≤zi≤αi,∀i∈W.The solution to subproblem (33) can be obtained in a similar way of L2y by observing the dynamics of L3z with respect to z.As shown in Figure4(a), when λiz<0, we have a concave function. However, different to the case of subproblem L2y, we solve minimization problem, and the critical point is no more optimal. And if λiz≥0, it shows monotonic nonincreasing, as seen in Figure 4(b). Therefore, we define the analytical solution of both cases as follows:(38)zi∗∈argminzL3αi,L3γi.However, γi in (35) has a variable u which is indeed u=∑i′∈Wxi′j. And this means that subproblem (33) is not separated from a variable x. Our heuristic remedy is to use u=W. Although this setting obviously will change the solution, it at least provides a feasible solution according to the following logic. Looking back on the previous discussion, zi is survival probability of weapon i. If u=W, zi has the greatest value. Let’s define such zi to be zmax. And we substitute zmax into (21), then following is satisfied.(39)∑i∈WVizumax≥∑i∈WVizopt≥β,where zopt is the solution where we set u=∑i′∈Wxi′jFigure 4
Dynamics ofSL3izi with different λ values. (a) λiz<0, (b) λiz>0.
(a)(b)


## 3.2.6. Line Search Method

Finding exact multipliers in Lagrangian relaxation is another optimization problem, and we introduce two algorithms for solving the Lagrangian dual problem. The first approach is to find a localλ solution using iterative first-order line search, as described in Algorithm 1.Algorithm 1: Line search method: WTAG_LS
(1)
Set iteration numberk←1, initial λ1←1,1,…,1,1(2)
Solve 3 subproblems (solution:x^,y^,z^ )(3)
Computepk (step direction)            ⊳pk=∇λk(4)
Computeαk (step size)              ⊳αk=γ1/k+1λkmax/pkmax, 0≤γ≤1(5)
ifpk>ϵthen                ⊳ϵ = tolerance(6)
λk+1←λk+αkpk(7)
Go line 2 andk←k+1(8)
else ifpk≤ϵthen(9)
STOP(10)
end ifIn algorithm 1,λk=λjy,λiz,λβ,∀i∈W,∀j∈T. Through numerical experiments, we found that the choice of step length αk had the greatest effect on the performance of the algorithm. We empirically decided to adopt the ideas of Fan et al. [24] and Feng and Liao [25].(40)αk=γ1k+1λkmaxξkmax,0≤γ≤1,where λkmax is the maximum of multiplier for coupling constraint, and ξkmax is a maximum of subgradient for coupling constraint. In this algorithm, as described in line 8, it terminates the search if the gradient norm (i.e., search direction) is less than epsilon, 0.2 in this study. Moreover, as a result of experiments, we empirically suggest using one vector for the initial choice of λk.


## 3.2.7. Subgradient Method

Another method for solving the Lagrangian dual problem is the subgradient method, and the pseudocode of the method for our problem is shown in Algorithm 2.In algorithm 2,λk=λjy,λiz,λβ,∀i∈W,∀j∈T at iteration step k. We define two initial lambda vectors nearly symmetric about the origin for the purpose of preventing the LP problem in line 5 is unbounded. The specific numbers in the initial lambda are determined based on parameters of instance which is described in Table 2. As a result of numerous experiments, we suggest to use minVj×max1−pij as the initial choice of lambda and −minVj×max1−pij for choice of second iterations. As explained in Section 3.2.1, since λβ≥0 has to be satisfied, we set the last element of the initial λ2 to zero [26]. As described in line 6, the algorithm stops if one of the two conditions is satisfied. One is to check if the value of the gradient norm (i.e., sk) gets closer to 0, and the second is to see if the L0 value doesn’t get any better for a certain period of iterations (i.e., L0k+1−L0∗/L0k−L0∗<ϵ where L0k is solution of k th iteration of LP problem in line 5 and L0∗ is best solution so far. In this study, ϵ=0.1) [27]. Our conjecture for both WTAG_LS and WTAG_SG is that first-order information (i.e., ∇λk) does not provide good information since we solve 3 subproblems in discrete manner.Algorithm 2: Subgradient Method: WTAG_SG
(1)
Set iteration numberk←1(2)
initialλ1←20,20,…,20,20, λ2←−20,−20,…,−20,0(3)
Solve 3 sub-problems (solution:x^,y^,z^ )(4)
Computesk (subgradient)               ⊳sk=∇λk(5)
Solve LP problem ofλ (solution: λ^ )max,L0,L0≤Lλk+skTλk+1−λk,(6)
ifsk≤ϵ or L0^ has no improvement then       ⊳ϵ = tolerance(7)
STOP(8)
else(9)
λk+1←λ^(10)
Go line 3,k←k+1, Add inequality to LP problem in line 5.(11)
end if


## 4. Computational Results

In order to verify the effectiveness of the algorithm, we conduct various experiments. Again, so far, we proposed three different variations of WTAG. All these models have in common that they are for solving WTAG, but WTAG_E is MINLP formulation, and both WTAG_LS and WTAG_SG are Lagrangian-based algorithms. We use BARON (version 20.4.14) with default setting (optimality gap 99% and CPU time limit 4000 sec) under GAMS (version 32.2) environment for solving WTAG_E. Two Lagrangian algorithms (i.e., WTAG_LS and WTAG_SG) are implemented using MATLAB (R2020b 9.9.0). All experiments are conducted on a CPU (Intel i5-10210U) with RAM (16 GB).


### 4.1. Instances

This section describes the instance used in the experiments. First, we create small-size random instances where parameters are uniformly distributed on a given bound described in Table2. We choose uniform distribution because military data should not be disclosed to the public, and uniform sampling is preferable for validating the robustness of the algorithm by avoiding being biased induced by data selection. The bound in Table 2 is also an arbitrary number, but we set a lower bound of pij,sij to 0.4 since the situation where probability of kill is very small is unrealistic. Examples of instances created by these rules are presented in Table 3.Table 2
Instance bound of experiments.
pij,sijVi,Vjβ0.4≤pij,sij≤160≤Vi,Vj≤1001/5∑j∈TVj≤β≤3/5∑j∈TVjTable 3
Instances in 2 weapons and 2 targets.
Instancespw1t1pw2t1pw1t2pw2t2sw1t1sw2t1sw1t2sw2t2Vw1Vw2Vt1Vt2βPar-2 : 2-10.70.40.60.70.70.40.60.88090807090Par-2 : 2-20.90.60.70.50.50.50.50.4100100100100100Par-2 : 2-30.40.50.90.80.80.80.70.96060907060Par-2 : 2-40.70.70.50.50.50.70.70.67090909060Par-2 : 2-50.40.80.80.50.50.50.50.58080808080Par-2 : 2-60.60.80.50.70.70.70.80.58095859090Par-2 : 2-70.80.90.70.60.70.50.50.4100100100100100Par-2 : 2-80.40.40.40.50.50.40.80.410060707080Par-2 : 2-90.80.70.40.90.40.60.50.57080908080Par-2 : 2-100.40.70.40.80.60.50.80.5100100100100100


### 4.2. Result of Small Instances

Computational results of 2 weapons and 2 targets instances are described in Table4. Table 4 shows the objective function value, solution, and CPU elapsed time of three different models. Because WTAG_E is a nonlinear formulation, it does not guarantee a global optimality. However, in small-size instances, as described in Table 4, it may provide an optimal solution. All WTAG_E solutions in Table 4 are turned out to be optimal (i.e., confirmed by enumerating all possible solutions). Besides, WTAG_E was the fastest in time. On the other hand, the two heuristic methods do not outperform in terms of optimality and time. The optimality of WTAG_LS and WTAG_SG was distributed between approximately 43% and 99.9%. In experiment of small size instances, we find that Lagrangian algorithms provide solid lower bound. However, it did not help much in terms of time or optimality.Table 4
Results of 2 weapons and 2 targets.
InstancesWTAG_EWTAG_LSWTAG_SGObjxTime (sec)ObjxTime (sec)Opt (%)ObjxTime (sec)Opt (%)Par-2 : 2-184.41, 1, 0, 00.1637.641, 0, 0, 13.0244.637.971, 0, 0, 16.2944.9Par-2 : 2-220.01, 0, 0, 10.1019.721, 0, 0, 13.5398.619.971, 0, 0, 17.8299.9Par-2 : 2-394.20, 0, 1, 10.1550.721, 1, 0, 02.9253.851.291, 1, 0, 05.7854.4Par-2 : 2-498.11, 1, 0, 00.2144.441, 0, 0, 13.1745.344.741, 1, 0, 06.2545.6Par-2 : 2-532.00, 1, 1, 00.1331.830, 1, 1, 03.1899.531.990, 1, 1, 07.6499.9Par-2 : 2-689.50,0,1,10.1238.111, 1, 0, 02.5642.638.621, 0, 0, 14.7843.2Par-2 : 2-740.00, 1, 1, 00.1038.210, 1, 1, 03.5895.539.840, 1, 1, 08.2099.6Par-2 : 2-877.01, 0, 0, 10.0974.241, 0, 0, 12.8096.476.991,0,0,17.3499.9Par-2 : 2-922.01, 0, 0, 10.1421.761, 0, 0, 12.5898.921.981, 0, 0, 16.7299.9Par-2 : 2-101120, 0, 1, 10.1065.611, 1, 0, 02.8158.571.151, 0, 0, 15.9266.2


### 4.3. Results of Large Instances

Results of large-size instances are described in Table5. Up to 6 weapons and 6 targets, the WTAG_E provided a local solution in the shortest CPU time. In this table, the objective value in the WTAG_E column means an integer feasible solution value. For WTAG_LS and WTAG_SG, it is a lower bound for true optimal value. As instance size increases, WTAG_E is not able to obtain a local solution whereas two heuristic methods provide a bound in a relatively shorter time.Table 5
Results of large-size instances 1.
InstancesSizeWTAG_EWTAG_LSWTAG_SGObjTime (sec)ObjTime (sec)ObjTime (sec)Par-3 : 33 : 332.020.5026.0321.1725.899.05Par-4 : 44 : 480.862.1620.7024.0235.2623.25Par-5 : 55 : 5111.818.9240.7433.5155.3535.80Par-6 : 66 : 690.5615.5326.7836.6132.8565.00Par-7 : 77 : 7104.6486.8126.7046.1940.00106.22Par-8 : 88 : 8—≥400013.6435.1919.47135.39Par-9 : 99 : 9—≥400023.2245.4539.66230.23In Table6, we describe the result of WTAG_LS and WTAG_SG for larger instances, and those are not solvable using WTAG_E. In this table, the objective value is the Lagrangian function value, and ∑j∈TVjyj in third column is objective value without the Lagrangian function which is actually the original objective function. The optimality gap in the fourth column calculates the gap between these two values. The reason why the value of ∑j∈TVjyj is shown in this table is that, if inequality (23) and (25) is satisfied, such value becomes the feasible bound of this problem. For example, with instance par-9:11-12, a solution from WTAG_LS satisfies inequality (23), and we can use 59.97 as a feasible bound for the original problem.Table 6
Results of large-size 2.
InstancesWTAG_LSWTAG_SGObjTime (sec)∑j∈TVjyjOpt (%)ObjTime (sec)∑j∈TVjyjOpt (%)Par-8 : 718.1031.9122.1781.625.87105.3426.5197.2Par-8 : 813.6435.1917.1779.419.47135.3919.7798.4Par-8 : 1059.2442.1160.5097.976.09199.8976.4099.6Par-9 : 820.6642.4822.1493.323.87127.0624.6396.9Par-9 : 923.2245.4533.8868.539.66230.2339.7099.4Par-9 : 1138.2443.9459.9763.869.85237.5073.9094.5Par-10 : 916.7446.4718.9788.222.20199.3025.0588.7Par-10 : 1034.6750.7348.8571.047.35241.2548.5697.5Par-10 : 1128.6258.5532.9686.839.87267.7143.3891.9Par-11 : 1123.0963.3929.0879.431.73409.8332.0099.1Par-12 : 1231.1281.0636.5085.334.00515.2536.0394.4Figure5 shows the dynamics of function value according to the CPU time of WTAG_LS and WTAG_SG. The WTAG_LS obtains solution in a relatively short time compared to WTAG_SG in all cases. Comparing the two algorithms, our conclusion is WTAG_LS which is superior in terms of time.Figure 5
Dynamics of function value with respect to time of two algorithms with large instances. (a) par-5.5, (b) par-6.6, (c) par-7:7, (d) par-8:8, (e) par-9:9.
(a)(b)(c)(d)(e)Figure6 shows examples of assignments in 2-dimensional spaces. The coordinates of the weapons and targets in the figure are given arbitrarily. And through visualization, we find some considerations that occur in real-world ground combat and that we need to reflect on in the model. For instance, since every weapon has range, assigning a target at a long distance (e.g., t2 and w6 in Figures 6(a) and 6(b)) is not allowed. And we can deal with these issues by adding or modifying the equations in the math model. For example, the rule of not making assignments that are too far apart can be implemented using a set of inequalities xij≤dij, ∀i∈W,j∈T where dij is predefined binary parameters that takes 1 if the distance between weapon i and target j is not too far. From this point of view, optimization model for solving weapon-target assignment is useful despite the difficulties it contains. In the practical aspect of the result, since a small ground unit consists of about 10 different types of weapons, the methods are usable referring to the WTAG_LS algorithm which can find the solution in about 1 minute or less.Figure 6
Weapon target assignment example. (a) par-7:7 (solved using WTAG LS), (b) par-7:7 (solved using WTAG SG).
(a)(b)


## 4.1. Instances

This section describes the instance used in the experiments. First, we create small-size random instances where parameters are uniformly distributed on a given bound described in Table2. We choose uniform distribution because military data should not be disclosed to the public, and uniform sampling is preferable for validating the robustness of the algorithm by avoiding being biased induced by data selection. The bound in Table 2 is also an arbitrary number, but we set a lower bound of pij,sij to 0.4 since the situation where probability of kill is very small is unrealistic. Examples of instances created by these rules are presented in Table 3.Table 2
Instance bound of experiments.
pij,sijVi,Vjβ0.4≤pij,sij≤160≤Vi,Vj≤1001/5∑j∈TVj≤β≤3/5∑j∈TVjTable 3
Instances in 2 weapons and 2 targets.
Instancespw1t1pw2t1pw1t2pw2t2sw1t1sw2t1sw1t2sw2t2Vw1Vw2Vt1Vt2βPar-2 : 2-10.70.40.60.70.70.40.60.88090807090Par-2 : 2-20.90.60.70.50.50.50.50.4100100100100100Par-2 : 2-30.40.50.90.80.80.80.70.96060907060Par-2 : 2-40.70.70.50.50.50.70.70.67090909060Par-2 : 2-50.40.80.80.50.50.50.50.58080808080Par-2 : 2-60.60.80.50.70.70.70.80.58095859090Par-2 : 2-70.80.90.70.60.70.50.50.4100100100100100Par-2 : 2-80.40.40.40.50.50.40.80.410060707080Par-2 : 2-90.80.70.40.90.40.60.50.57080908080Par-2 : 2-100.40.70.40.80.60.50.80.5100100100100100


## 4.2. Result of Small Instances

Computational results of 2 weapons and 2 targets instances are described in Table4. Table 4 shows the objective function value, solution, and CPU elapsed time of three different models. Because WTAG_E is a nonlinear formulation, it does not guarantee a global optimality. However, in small-size instances, as described in Table 4, it may provide an optimal solution. All WTAG_E solutions in Table 4 are turned out to be optimal (i.e., confirmed by enumerating all possible solutions). Besides, WTAG_E was the fastest in time. On the other hand, the two heuristic methods do not outperform in terms of optimality and time. The optimality of WTAG_LS and WTAG_SG was distributed between approximately 43% and 99.9%. In experiment of small size instances, we find that Lagrangian algorithms provide solid lower bound. However, it did not help much in terms of time or optimality.Table 4
Results of 2 weapons and 2 targets.
InstancesWTAG_EWTAG_LSWTAG_SGObjxTime (sec)ObjxTime (sec)Opt (%)ObjxTime (sec)Opt (%)Par-2 : 2-184.41, 1, 0, 00.1637.641, 0, 0, 13.0244.637.971, 0, 0, 16.2944.9Par-2 : 2-220.01, 0, 0, 10.1019.721, 0, 0, 13.5398.619.971, 0, 0, 17.8299.9Par-2 : 2-394.20, 0, 1, 10.1550.721, 1, 0, 02.9253.851.291, 1, 0, 05.7854.4Par-2 : 2-498.11, 1, 0, 00.2144.441, 0, 0, 13.1745.344.741, 1, 0, 06.2545.6Par-2 : 2-532.00, 1, 1, 00.1331.830, 1, 1, 03.1899.531.990, 1, 1, 07.6499.9Par-2 : 2-689.50,0,1,10.1238.111, 1, 0, 02.5642.638.621, 0, 0, 14.7843.2Par-2 : 2-740.00, 1, 1, 00.1038.210, 1, 1, 03.5895.539.840, 1, 1, 08.2099.6Par-2 : 2-877.01, 0, 0, 10.0974.241, 0, 0, 12.8096.476.991,0,0,17.3499.9Par-2 : 2-922.01, 0, 0, 10.1421.761, 0, 0, 12.5898.921.981, 0, 0, 16.7299.9Par-2 : 2-101120, 0, 1, 10.1065.611, 1, 0, 02.8158.571.151, 0, 0, 15.9266.2


## 4.3. Results of Large Instances

Results of large-size instances are described in Table5. Up to 6 weapons and 6 targets, the WTAG_E provided a local solution in the shortest CPU time. In this table, the objective value in the WTAG_E column means an integer feasible solution value. For WTAG_LS and WTAG_SG, it is a lower bound for true optimal value. As instance size increases, WTAG_E is not able to obtain a local solution whereas two heuristic methods provide a bound in a relatively shorter time.Table 5
Results of large-size instances 1.
InstancesSizeWTAG_EWTAG_LSWTAG_SGObjTime (sec)ObjTime (sec)ObjTime (sec)Par-3 : 33 : 332.020.5026.0321.1725.899.05Par-4 : 44 : 480.862.1620.7024.0235.2623.25Par-5 : 55 : 5111.818.9240.7433.5155.3535.80Par-6 : 66 : 690.5615.5326.7836.6132.8565.00Par-7 : 77 : 7104.6486.8126.7046.1940.00106.22Par-8 : 88 : 8—≥400013.6435.1919.47135.39Par-9 : 99 : 9—≥400023.2245.4539.66230.23In Table6, we describe the result of WTAG_LS and WTAG_SG for larger instances, and those are not solvable using WTAG_E. In this table, the objective value is the Lagrangian function value, and ∑j∈TVjyj in third column is objective value without the Lagrangian function which is actually the original objective function. The optimality gap in the fourth column calculates the gap between these two values. The reason why the value of ∑j∈TVjyj is shown in this table is that, if inequality (23) and (25) is satisfied, such value becomes the feasible bound of this problem. For example, with instance par-9:11-12, a solution from WTAG_LS satisfies inequality (23), and we can use 59.97 as a feasible bound for the original problem.Table 6
Results of large-size 2.
InstancesWTAG_LSWTAG_SGObjTime (sec)∑j∈TVjyjOpt (%)ObjTime (sec)∑j∈TVjyjOpt (%)Par-8 : 718.1031.9122.1781.625.87105.3426.5197.2Par-8 : 813.6435.1917.1779.419.47135.3919.7798.4Par-8 : 1059.2442.1160.5097.976.09199.8976.4099.6Par-9 : 820.6642.4822.1493.323.87127.0624.6396.9Par-9 : 923.2245.4533.8868.539.66230.2339.7099.4Par-9 : 1138.2443.9459.9763.869.85237.5073.9094.5Par-10 : 916.7446.4718.9788.222.20199.3025.0588.7Par-10 : 1034.6750.7348.8571.047.35241.2548.5697.5Par-10 : 1128.6258.5532.9686.839.87267.7143.3891.9Par-11 : 1123.0963.3929.0879.431.73409.8332.0099.1Par-12 : 1231.1281.0636.5085.334.00515.2536.0394.4Figure5 shows the dynamics of function value according to the CPU time of WTAG_LS and WTAG_SG. The WTAG_LS obtains solution in a relatively short time compared to WTAG_SG in all cases. Comparing the two algorithms, our conclusion is WTAG_LS which is superior in terms of time.Figure 5
Dynamics of function value with respect to time of two algorithms with large instances. (a) par-5.5, (b) par-6.6, (c) par-7:7, (d) par-8:8, (e) par-9:9.
(a)(b)(c)(d)(e)Figure6 shows examples of assignments in 2-dimensional spaces. The coordinates of the weapons and targets in the figure are given arbitrarily. And through visualization, we find some considerations that occur in real-world ground combat and that we need to reflect on in the model. For instance, since every weapon has range, assigning a target at a long distance (e.g., t2 and w6 in Figures 6(a) and 6(b)) is not allowed. And we can deal with these issues by adding or modifying the equations in the math model. For example, the rule of not making assignments that are too far apart can be implemented using a set of inequalities xij≤dij, ∀i∈W,j∈T where dij is predefined binary parameters that takes 1 if the distance between weapon i and target j is not too far. From this point of view, optimization model for solving weapon-target assignment is useful despite the difficulties it contains. In the practical aspect of the result, since a small ground unit consists of about 10 different types of weapons, the methods are usable referring to the WTAG_LS algorithm which can find the solution in about 1 minute or less.Figure 6
Weapon target assignment example. (a) par-7:7 (solved using WTAG LS), (b) par-7:7 (solved using WTAG SG).
(a)(b)


## 5. Conclusion

The WTA problem has been studied for a long time in the optimization field, but there have been few attempts to apply it to the ground battle. The purpose of this study is that (1) the new WTAG, WTA for Ground, the problem is formulated by considering the characteristics of ground combat, and (2) the most efficient way to solve this problem was presented.First, the WTAG problem reflected the important characteristic that engagement should be considered when assigning a target, and this was reflected through a parameter calledβ that the commander can select according to the operational situations. Second, the Lagrangian relaxation is used to solve the problem efficiently, and then, we divided the problem into three subproblems. To obtain the Lagrangian multipliers, the line search method and the subgradient method are implicated. There are three computational methods proposed to solve this problem. WTAG_E is MINLP, and as the size of the problem increases, it is difficult to obtain a solution. WTAG_LS and WTAG_SG can both provide a dual bound of the problem, even a primal bound under specific conditions. However, not all methods provide global solutions.Future research directions are as follows. As discussed a bit in Section 4, there are many characteristics of small unit ground combat that need to be further considered in order for the allocation derived from this problem to be put to practical use (e.g., range of the weapon, line of sight, closeness, cross-shoot). Lastly, the Lagrangian-based method can find the solutions that provide bound, but it is not fast enough to find the solution in real-time meaning that the performance of the algorithm in terms of time needs to be improved.



---

*Source: 9228993-2022-12-24.xml*