此文轉載自：https://blog.csdn.net/qq_35008055/article/details/110201571#commentBox
手把手教你用Gurobi求解一個數學模型
手把手教你用Gurobi求解一個數學模型
在接觸Gurobi之前，一直使用Java語言調用cplex求解數學模型，這段時間在師兄的指點下，學習了使用python調用Gurobi的一些基礎操作，感嘆實在是太簡易了。
在此分享一個求解Vrptw問題的小例子。
帶時間窗的車輛路徑規划問題(Vrptw)
對於Vrptw問題來說，數學模型主要由以下部分組成。首先我們定義一些相關參數，一個圖可以表示為 G ( V , A ) G(V,A) G(V,A)，其中 V = { 0 , 1 , . . . , n , n + 1 } V= \{ 0,1,...,n,n+1 \} V={0,1,...,n,n+1}為圖中所有點的集合， A A A為圖中所有弧的集合，有 ( i , j ) ∈ (i,j)\in (i,j)∈A， ∀ i , j ∈ V , i ≠ j \forall i,j\in V, i\ne j ∀i,j∈V,i=j。弧 ( i , j ) (i,j) (i,j)的單位運輸費用為 c i j c_{ij} cij，運輸時間為 t i j t_{ij} tij，每個客戶點的需求為 q i j q_{ij} qij，可服務的時間窗為 [ e i , l i ] [e_i,l_i] [ei,li]，服務時長為 s e r v i serv_i servi。令車輛的集合為 K K K，每輛車的最大載重為 Q k , ∀ k ∈ K Q_k,\forall k\in K Qk,∀k∈K。決策變量為 x i j k x_{ij}^k xijk，代表第 k k k輛車是否服務了弧 ( i , j ) (i,j) (i,j)； s i s_i si為客戶點 i i i開始被服務的時間。
接下來構建數學模型。
目標函數為最小化運輸成本：
m i n ∑ k ∈ K ∑ ( i , j ) ∈ A c i j x i j (1) min \sum_{k\in K} \sum_{(i,j)\in A}c_{ij}x_{ij} \tag{1} mink∈K∑(i,j)∈A∑cijxij(1)
約束一讓車輛駛出倉庫(depot)：
∑ ( 0 , j ) ∈ A x 0 j k = 1 , ∀ k ∈ K (2) \sum_{(0,j)\in A}x_{0j}^k =1,\ \forall k\in K \tag{2} (0,j)∈A∑x0jk=1, ∀k∈K(2)
約束二為流平衡(除去depot之外的點)：
∑ ( i , j ) ∈ A x i j k − ∑ ( j , i ) ∈ A x j i k = 0 , ∀ k ∈ K , i ∈ V ∖ { s , t } (3) \sum_{(i,j)\in A}x_{ij}^k - \sum_{(j,i)\in A}x_{ji}^k = 0,\ \forall k\in K ,i \in V\setminus \{s,t\} \tag{3} (i,j)∈A∑xijk−(j,i)∈A∑xjik=0, ∀k∈K,i∈V∖{s,t}(3)
約束三讓車輛駛回depot：
∑ ( i , n + 1 ) ∈ A x i , n + 1 k = 1 , ∀ k ∈ K (4) \sum_{(i,n+1)\in A}x_{i,n+1}^k =1,\ \forall k\in K \tag{4} (i,n+1)∈A∑xi,n+1k=1, ∀k∈K(4)
約束四保證每個客戶點都被服務：
∑ k ∈ K ∑ ( i , j ) ∈ A x i , j k = 1 , ∀ i ∈ V ∖ { s , t } (5) \sum_{k \in K}\sum_{(i,j)\in A}x_{i,j}^k =1,\ \forall i \in V\setminus \{s,t\} \tag{5} k∈K∑(i,j)∈A∑xi,jk=1, ∀i∈V∖{s,t}(5)
約束五保證被服務的相鄰節點開始服務時間的大小關系(去回路)：
s i + t i j + s e r v i − M ( 1 − x i j ) ≤ s j , ∀ ( i , j ) ∈ A (6) s_i+t_{ij}+serv_i-M(1-x_{ij})\le s_j,\ \forall (i,j) \in A \tag{6} si+tij+servi−M(1−xij)≤sj, ∀(i,j)∈A(6)
約束六保證不違反客戶的時間窗：
e i ≤ s i ≤ l i , ∀ i ∈ V ∖ { s , t } (7) e_i \le s_i \le l_i ,\ \forall i \in V\setminus \{s,t\} \tag{7} ei≤si≤li, ∀i∈V∖{s,t}(7)
約束七保證不違反車輛的載重約束：
∑ ( i , j ) ∈ A x i j k q i ≤ Q k , ∀ k ∈ K (8) \sum_{(i,j)\in A}x_{ij}^kq_i \le Q_k ,\ \forall k \in K \tag{8} (i,j)∈A∑xijkqi≤Qk, ∀k∈K(8)
最后是決策變量的取值約束：
x i j k ∈ { 0 , 1 } ∀ ( i , j ) ∈ A , k ∈ K s i ≥ 0 , ∀ i ∈ V ∖ { s , t } (9) x_{ij}^k\in \{0,1\} \ \forall (i,j) \in A,k \in K \\ s_i \ge 0,\ \forall i \in V\setminus \{s,t\} \tag{9} xijk∈{0,1} ∀(i,j)∈A,k∈Ksi≥0, ∀i∈V∖{s,t}(9)
那么(1)~(9)式就組成了Vrptw問題的數學模型。
python調用Gurobi求解Vrptw
首先我們定義一下需要用到的參數：
class Data:
customerNum = 0
nodeNum = 0
vehicleNum = 0
capacity = 0
cor_X = []
cor_Y = []
demand = []
serviceTime = []
readyTime = []
dueTime = []
disMatrix = [[]]
定義一個讀取數據的函數，並對節點之間的距離進行計算：
# function to read data from .txt files
def readData(data, path, customerNum):
data.customerNum = customerNum
data.nodeNum = customerNum + 2
f = open(path, 'r')
lines = f.readlines()
count = 0
# read the info
for line in lines:
count = count + 1
if (count == 5):
line = line[:-1].strip()
str = re.split(r" +", line)
data.vehicleNum = int(str[0])
data.capacity = float(str[1])
elif (count >= 10 and count <= 10 + customerNum):
line = line[:-1]
str = re.split(r" +", line)
data.cor_X.append(float(str[2]))
data.cor_Y.append(float(str[3]))
data.demand.append(float(str[4]))
data.readyTime.append(float(str[5]))
data.dueTime.append(float(str[6]))
data.serviceTime.append(float(str[7]))
data.cor_X.append(data.cor_X[0])
data.cor_Y.append(data.cor_Y[0])
data.demand.append(data.demand[0])
data.readyTime.append(data.readyTime[0])
data.dueTime.append(data.dueTime[0])
data.serviceTime.append(data.serviceTime[0])
# compute the distance matrix
data.disMatrix = [([0] * data.nodeNum) for p in range(data.nodeNum)] # 初始化距離矩陣的維度,防止淺拷貝
# data.disMatrix = [[0] * nodeNum] * nodeNum]; 這個是淺拷貝，容易重復
for i in range(0, data.nodeNum):
for j in range(0, data.nodeNum):
temp = (data.cor_X[i] - data.cor_X[j]) ** 2 + (data.cor_Y[i] - data.cor_Y[j]) ** 2
data.disMatrix[i][j] = math.sqrt(temp)
temp = 0
return data
讀取數據，並定義一些參數：
data = Data()
path = 'c101.txt' //讀取Solomon數據集
customerNum = 20 //設置客戶數量
readData(data, path, customerNum)
data.vehicleNum = 10 //設置車輛數
printData(data, customerNum)
BigM = 100000 //定義一個極大值
調用gurobi進行模型的建立與求解：
x = {}
s = {} //定義字典用來存放決策變量
根據式(9)定義決策變量，並加入模型當中：
for i in range(data.nodeNum):
for k in range(data.vehicleNum):
name = 's_' + str(i) + '_' + str(k)
s[i,k] = model.addVar(0
, 1500
, vtype= GRB.CONTINUOUS
, name= name) //定義訪問時間為連續變量
for j in range(data.nodeNum):
if(i != j):
name = 'x_' + str(i) + '_' + str(j) + '_' + str(k)
x[i,j,k] = model.addVar(0
, 1
, vtype= GRB.BINARY
, name= name) //定義是否服務為0-1變量
根據式(1)定義目標函數，並加入模型當中：
//首先定義一個線性表達式
obj = LinExpr(0)
for i in range(data.nodeNum):
for k in range(data.vehicleNum):
for j in range(data.nodeNum):
if(i != j):
//將目標函數系數與決策變量相乘，並進行連加
obj.addTerms(data.disMatrix[i][j], x[i,j,k])
//將表示目標函數的線性表達式加入模型，並定義為求解最小化問題
model.setObjective(obj, GRB.MINIMIZE)
根據式(2)~(8)定義決策變量，並加入模型當中：
約束一（vehicle_depart）：
for k in range(data.vehicleNum):
//同樣先定義一個線性表達式
lhs = LinExpr(0)
for j in range(data.nodeNum):
if(j != 0):
//約束系數與決策變量相乘
lhs.addTerms(1, x[0,j,k])
//將約束加入模型
model.addConstr(lhs == 1, name= 'vehicle_depart_' + str(k))
約束二（flow_conservation）：
for k in range(data.vehicleNum):
for h in range(1, data.nodeNum - 1):
expr1 = LinExpr(0)
expr2 = LinExpr(0)
for i in range(data.nodeNum):
if (h != i):
expr1.addTerms(1, x[i,h,k])
for j in range(data.nodeNum):
if (h != j):
expr2.addTerms(1, x[h,j,k])
model.addConstr(expr1 == expr2, name= 'flow_conservation_' + str(i))
expr1.clear()
expr2.clear()
約束三（vehicle_enter）：
for k in range(data.vehicleNum):
lhs = LinExpr(0)
for j in range(data.nodeNum - 1):
if(j != 0):
lhs.addTerms(1, x[j, data.nodeNum-1, k])
model.addConstr(lhs == 1, name= 'vehicle_enter_' + str(k))
約束四（customer_visit）：
for i in range(1, data.nodeNum - 1):
lhs = LinExpr(0)
for k in range(data.vehicleNum):
for j in range(1, data.nodeNum):
if(i != j):
lhs.addTerms(1, x[i,j,k])
model.addConstr(lhs == 1, name= 'customer_visit_' + str(i))
約束五（time_windows）：
for k in range(data.vehicleNum):
for i in range(data.nodeNum):
for j in range(data.nodeNum):
if(i != j):
model.addConstr(s[i,k] + data.disMatrix[i][j] + data.serviceTime[i] - s[j,k]- BigM + BigM * x[i,j,k] <= 0 , name= 'time_windows_')
約束六（ready_time and due_time）：
for i in range(1,data.nodeNum-1):
for k in range(data.vehicleNum):
model.addConstr(data.readyTime[i] <= s[i,k], name= 'ready_time')
model.addConstr(s[i,k] <= data.dueTime[i], name= 'due_time')
約束七（capacity_vehicle）：
for k in range(data.vehicleNum):
lhs = LinExpr(0)
for i in range(1, data.nodeNum - 1):
for j in range(data.nodeNum):
if(i != j):
lhs.addTerms(data.demand[i], x[i,j,k])
model.addConstr(lhs <= data.capacity, name= 'capacity_vehicle' + str(k))
求解模型（solve）並輸出解：
model.optimize()
print("\n\n-----optimal value-----")
print(model.ObjVal)
for key in x.keys():
if(x[key].x > 0 ):
print(x[key].VarName + ' = ', x[key].x)
導出模型：
model.write('VRPTW.lp')
求解結果（部分）：
這樣就完成了，感謝大家的閱讀。
作者：夏暘，清華大學，工業工程系/深圳國際研究生院 (碩士在讀)
劉興祿，清華大學，清華伯克利深圳學院（博士在讀）
完整Project文件請關注運小籌公眾號，回復【Vrptw】獲取。
