前言

靜態(tài)武器目標(biāo)分配（static weapon-target assignment,，WTA）問(wèn)題是軍事領(lǐng)域的著名問(wèn)題,，也是分配問(wèn)題中的研究熱點(diǎn)之一。SWTA問(wèn)題易于描述但是難以求解,，自1958年提出以來(lái)，大量學(xué)者對(duì)其進(jìn)行了研究,，但是對(duì)于大規(guī)模的SWTA問(wèn)題求解直至近兩年才有較大突破,，本文將借助Gurobi精確求解小規(guī)模的SWTA問(wèn)題，借此進(jìn)一步熟悉Gurobi建模過(guò)程,。

SWTA問(wèn)題

靜態(tài)武器目標(biāo)的描述及其簡(jiǎn)單非線性模型如下：

非線性模型難以求解,，因此首先對(duì)其進(jìn)行線性化處理，轉(zhuǎn)化成如下線性模型：

其中s_j為m位的二進(jìn)制數(shù),，表示m個(gè)武器在目標(biāo)j上的分配情況,，假如我方有3個(gè)武器，那么[0,0,0]表示沒(méi)有武器分配給目標(biāo)j,，[0,1,0]表示僅有第二個(gè)武器分配給目標(biāo)j,，[1,1,1]表示三個(gè)武器都分配給目標(biāo)j。

建模與求解過(guò)程

建模和求解通過(guò)Python調(diào)用Gurobi實(shí)現(xiàn),。首先是求解過(guò)程中用到的包,，主要是gurobipy。

import pandas as pd
from gurobi import *
import random
import copy

第二部分是數(shù)據(jù)讀取,，數(shù)據(jù)來(lái)源于txt文件,，以一個(gè)武器數(shù)為5，目標(biāo)數(shù)為3的小規(guī)模問(wèn)題為例,，數(shù)據(jù)如下：

#data.txt
5
3
0.3   0.6   0.5
0.441 0.424 0.401
0.468 0.449 0.422
0.462 0.445 0.420
0.433 0.420 0.400
0.467 0.452 0.428

其中第一行為武器數(shù),，第二行為目標(biāo)數(shù),第三行為目標(biāo)威脅度，第四至八行為各武器對(duì)目標(biāo)的毀傷概率,，數(shù)據(jù)處理的代碼如下：

data=pd.read_table('data/5_3.txt',header=None)
m=int(data.iloc[0,0])
n=int(data.iloc[1,0])
v=[float(item) for item in data.iloc[2,0].split('   ')]
p=[[float(item[j]) for j in range(len(item))] for item in [item.split(' ') for item in data.iloc[3:,0].values]]

為了后續(xù)建模方便,，接下來(lái)生成所有可能的分配狀態(tài)s

def enumeration(m,n,b,bits):
    if len(b)==m:
        bits.append(copy.deepcopy(b))
        return
    for i in range(2):
        b.append(i)
        enumeration(m,n,b,bits)
        b.pop()
def get_bit(m,n):
    b=[]
    bits=[]
    enumeration(m,n,b,bits)
    return bits

def get_destroy(j,bit):
    survive=v[j]
    for i in range(len(bit)):
        if bit[i]==0:continue
        survive*=(1-p[i][j])
    return v[j]-survive
    
bits=get_bit(m,n)

其中bits的取值如下：

[[0, 0, 0, 0, 0],
 [0, 0, 0, 0, 1],
 [0, 0, 0, 1, 0],
 [0, 0, 0, 1, 1],
 [0, 0, 1, 0, 0],
 [0, 0, 1, 0, 1],
 [0, 0, 1, 1, 0],
 [0, 0, 1, 1, 1],
 [0, 1, 0, 0, 0],
 [0, 1, 0, 0, 1],
 [0, 1, 0, 1, 0],
 [0, 1, 0, 1, 1],
 [0, 1, 1, 0, 0],
 [0, 1, 1, 0, 1],
 [0, 1, 1, 1, 0],
 [0, 1, 1, 1, 1],
 [1, 0, 0, 0, 0],
 [1, 0, 0, 0, 1],
 [1, 0, 0, 1, 0],
 [1, 0, 0, 1, 1],
 [1, 0, 1, 0, 0],
 [1, 0, 1, 0, 1],
 [1, 0, 1, 1, 0],
 [1, 0, 1, 1, 1],
 [1, 1, 0, 0, 0],
 [1, 1, 0, 0, 1],
 [1, 1, 0, 1, 0],
 [1, 1, 0, 1, 1],
 [1, 1, 1, 0, 0],
 [1, 1, 1, 0, 1],
 [1, 1, 1, 1, 0],
 [1, 1, 1, 1, 1]]

數(shù)據(jù)準(zhǔn)備完畢，接下來(lái)依據(jù)之前給出來(lái)的線性模型建模即可,，建模過(guò)程如下：

# 模型
model=Model('WTA')

# 變量
y=model.addVars([(j,s) for j in range(n) for s in range(2**m)],name='y')

#目標(biāo)函數(shù)
model.setObjective(quicksum(quicksum(get_destroy(j,bits[s])*y[j,s] for s in range(2**m))                            for j in range(n)),GRB.MAXIMIZE)

#約束
#(1)對(duì)每個(gè)目標(biāo)j,，只有一個(gè)狀態(tài)
model.addConstrs(((quicksum(y[j,s] for s in range(2**m))==1) for j in range(n)),name='state constraint')
#(2)對(duì)每個(gè)武器i，只能分配給一個(gè)目標(biāo)
model.addConstrs(((quicksum(quicksum(y[j,s]*bits[s][i] for j in range(n)) for s in range(2**m))<=1)                  for i in range(m)),name='assign constraint')

#求解
model.params.OutputFlag=0
model.optimize()

#輸出
print('Obj={}'.format(model.objVal))
for j in range(n):
    for s in range(2**m):
        if y[j,s].x:
            weapon=[]
            for i in range(len(bits[s])):
                if bits[s][i]:
                    weapon.append(i+1)
            print('weapon {} assigned to target [{}]'.format(weapon,j+1))

完整代碼與求解結(jié)果

將上述過(guò)程整合，完整代碼如下：

import pandas as pd
from gurobi import *
import random
import copy


data=pd.read_table('data/5_3.txt',header=None)
m=int(data.iloc[0,0])
n=int(data.iloc[1,0])
v=[float(item) for item in data.iloc[2,0].split('   ')]
p=[[float(item[j]) for j in range(len(item))] for item in [item.split(' ') for item in data.iloc[3:,0].values]]


def enumeration(m,n,b,bits):
    if len(b)==m:
        bits.append(copy.deepcopy(b))
        return
    for i in range(2):
        b.append(i)
        enumeration(m,n,b,bits)
        b.pop()
def get_bit(m,n):
    b=[]
    bits=[]
    enumeration(m,n,b,bits)
    return bits

def get_destroy(j,bit):
    survive=v[j]
    for i in range(len(bit)):
        if bit[i]==0:continue
        survive*=(1-p[i][j])
    return v[j]-survive
    
bits=get_bit(m,n)


# 模型
model=Model('WTA')

# 變量
y=model.addVars([(j,s) for j in range(n) for s in range(2**m)],name='y')

#目標(biāo)函數(shù)
model.setObjective(quicksum(quicksum(get_destroy(j,bits[s])*y[j,s] for s in range(2**m))                            for j in range(n)),GRB.MAXIMIZE)

#約束
#(1)對(duì)每個(gè)目標(biāo)j,，只有一個(gè)狀態(tài)
model.addConstrs(((quicksum(y[j,s] for s in range(2**m))==1) for j in range(n)),name='state constraint')
#(2)對(duì)每個(gè)武器i,，只能分配給一個(gè)目標(biāo)
model.addConstrs(((quicksum(quicksum(y[j,s]*bits[s][i] for j in range(n)) for s in range(2**m))<=1)                  for i in range(m)),name='assign constraint')

#求解
model.params.OutputFlag=0
model.optimize()

#輸出
print('Obj={}'.format(model.objVal))
for j in range(n):
    for s in range(2**m):
        if y[j,s].x:
            weapon=[]
            for i in range(len(bits[s])):
                if bits[s][i]:
                    weapon.append(i+1)
            print('weapon {} assigned to target [{}]'.format(weapon,j+1))

模型求解結(jié)果如下，給出了我方的最大毀傷概率0.878216,，毀傷比例為0.878216/(0.3+0.5+0.6)=62.73%,，對(duì)應(yīng)的分配方案也已給出。

Obj=0.878216
weapon [2] assigned to target [1]
weapon [3, 5] assigned to target [2]
weapon [1, 4] assigned to target [3]

一類武器有多個(gè)的情況

上一節(jié)中求解的問(wèn)題是,，每個(gè)武器當(dāng)成一類,，即每類武器的數(shù)量均為1，在實(shí)際情況中,，武器種類往往是不多,，每一類武器會(huì)有多個(gè)，接下來(lái)將求解一類武器有多個(gè)的WTA問(wèn)題,，數(shù)據(jù)與上一節(jié)同,，將每一類武器的數(shù)量從1變?yōu)?，建模求解過(guò)程如下：

import pandas as pd
import numpy as np
from gurobipy import *
import copy

#原始數(shù)據(jù)
data=pd.read_table('data/5_3.txt',header=None)
m=int(data.iloc[0,0])
weapon_num=[2 for i in range(m)]
n=int(data.iloc[1,0])
v=[float(item) for item in data.iloc[2,0].split('   ')]
p=[[float(item[j]) for j in range(len(item))] for item in [item.split(' ') for item in data.iloc[3:,0].values]]

#中間數(shù)據(jù)
bits=[]
def next_bit(bit):
    for i in range(m):
        if bit[i]==weapon_num[i]:
            continue
        bit[i]+=1
        for j in range(i):
            bit[j]=0
        break
    return bit

def get_destroy(j,bit):
    survive=1
    for i in range(len(bit)):
        survive*=pow(1-p[i][j],bit[i])
    return v[j]*(1-survive)

S=np.prod([i+1 for i in weapon_num])
bit=[0 for i in range(m)]
bits=[bit]
while bit!=next_bit(copy.deepcopy(bit)):
    bit=next_bit(copy.deepcopy(bit))
    bits.append(bit)

#模型
model=Model('WTA')

#變量
y=model.addVars([(j,s) for j in range(n) for s in range(S)],                obj=[get_destroy(j,bits[s]) for j in range(n) for s in range(S)],vtype='B',name='y')

#目標(biāo)
model.setAttr(GRB.Attr.ModelSense,-1)

#約束
#(1)對(duì)每個(gè)目標(biāo)j,，只有一個(gè)狀態(tài)
model.addConstrs((quicksum(y[j,s] for s in range(S))==1 for j in range(n)),name='state constraint')
#(2)對(duì)每個(gè)武器i,，只能分配給一個(gè)目標(biāo)
model.addConstrs((quicksum(quicksum(y[j,s]*bits[s][i] for j in range(n))                          for s in range(S))<=weapon_num[i] for i in range(m)),name='assign constraint')
#求解
model.setParam(GRB.Param.OutputFlag,0)
model.optimize()

#輸出
print('毀傷值為{}'.format(model.objVal))
for j in range(n):
    for s in range(S):
        if y[j,s].x:
            weapon=[]
            for i in range(len(bits[s])):
                if bits[s][i]:
                    print('武器{}分配給目標(biāo){}的數(shù)量為{}'.format(i+1,j+1,bits[s][i]))

結(jié)果如下：

毀傷值為1.19508327488
武器2分配給目標(biāo)1的數(shù)量為2
武器1分配給目標(biāo)2的數(shù)量為1
武器3分配給目標(biāo)2的數(shù)量為2
武器5分配給目標(biāo)2的數(shù)量為1
武器1分配給目標(biāo)3的數(shù)量為1
武器4分配給目標(biāo)3的數(shù)量為2
武器5分配給目標(biāo)3的數(shù)量為1

總結(jié)

對(duì)于小規(guī)模的SWTA問(wèn)題（20武器20目標(biāo)以內(nèi)），我們基本都可以使用求解器快速求解,，但是對(duì)于大規(guī)模的問(wèn)題,，這里仍然存在指數(shù)爆炸的問(wèn)題，變量數(shù)n*2^m隨著武器數(shù)量的增加呈指數(shù)中增長(zhǎng),，試想一下,，100武器100目標(biāo)的SWTA問(wèn)題，變量數(shù)將達(dá)到天文數(shù)字,，那么我們?nèi)绾吻蠼?？是否還能求解,？答案是肯定的,，具體細(xì)節(jié)見(jiàn)論文A new exact algorithm for the Weapon-Target Assignment problem，論文精確求解了400武器400目標(biāo)的超大規(guī)模SWTA問(wèn)題,。