手把手教你用Gurobi求解一个数学模型

在接触Gurobi之前，一直使用Java语言调用cplex求解数学模型，这段时间在师兄的指点下，学习了使用python调用Gurobi的一些基础操作，感叹实在是太简易了。
在此分享一个求解Vrptw问题的小例子。

带时间窗的车辆路径规划问题(Vrptw)

对于Vrptw问题来说，数学模型主要由以下部分组成。首先我们定义一些相关参数，一个图可以表示为$G(V,A)$，其中$V=\{0,1,...,n,n+1\}$为图中所有点的集合，$A$为图中所有弧的集合，有$(i,j)\in$A，$\forall i,j\in V,i\ne j$。弧$(i,j)$的单位运输费用为$c_{ij}$，运输时间为$t_{ij}$，每个客户点的需求为$q_{ij}$，可服务的时间窗为$[e_i,l_i]$，服务时长为$ser{v}_i$。令车辆的集合为$K$，每辆车的最大载重为$Q_k,\forall k\in K$。决策变量为$x_{ij}^k$，代表第$k$辆车是否服务了弧$(i,j)$；$s_i$为客户点$i$开始被服务的时间。

接下来构建数学模型。
目标函数为最小化运输成本：
$$min\sum _{k\in K}\sum _{(i,j)\in A}{c}_{ij}{x}_{ij}\text{\hspace{1em}(1)}$$
约束一让车辆驶出仓库(depot)：
$$\sum _{(0,j)\in A}{x}_{0j}^k=1,\forall k\in K\text{\hspace{1em}(2)}$$
约束二为流平衡(除去depot之外的点)：
$$\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\}\text{\hspace{1em}(3)}$$
约束三让车辆驶回depot：
$$\sum _{(i,n+1)\in A}{x}_{i,n+1}^k=1,\forall k\in K\text{\hspace{1em}(4)}$$
约束四保证每个客户点都被服务：
$$\sum _{k\in K}\sum _{(i,j)\in A}{x}_{i,j}^k=1,\forall i\in V\setminus \{s,t\}\text{\hspace{1em}(5)}$$
约束五保证被服务的相邻节点开始服务时间的大小关系(去回路)：
$$s_i+t_{ij}+ser{v}_i-M(1-x_{ij})\le s_j,\forall (i,j)\in A\text{\hspace{1em}(6)}$$
约束六保证不违反客户的时间窗：
$$e_i\le s_i\le l_i,\forall i\in V\setminus \{s,t\}\text{\hspace{1em}(7)}$$
约束七保证不违反车辆的载重约束：
$$\sum _{(i,j)\in A}{x}_{ij}^k{q}_i\le Q_k,\forall k\in K\text{\hspace{1em}(8)}$$
最后是决策变量的取值约束：
$$\begin{gathered} 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\}\text{\hspace{1em}(9)} \end{gathered}$$
那么(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   = [[]] 
  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
 

定义一个读取数据的函数，并对节点之间的距离进行计算：

# 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

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
 

读取数据，并定义一些参数：

data = Data()
path = 'c101.txt' //读取Solomon数据集
customerNum = 20  //设置客户数量
readData(data, path, customerNum)
data.vehicleNum = 10  //设置车辆数
printData(data, customerNum)
BigM = 100000 //定义一个极大值

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
 

调用gurobi进行模型的建立与求解：

x = {}
s = {}  //定义字典用来存放决策变量

  

 

  
1
  
2
 

根据式(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
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
 

根据式(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)  

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
 

根据式(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)) 

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
 

约束二（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()

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
 

约束三（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))

  

 

  
1
  
2
  
3
  
4
  
5
  
6
 

约束四（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)) 

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
 

约束五（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_')

  

 

  
1
  
2
  
3
  
4
  
5
 

约束六（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')

  

 

  
1
  
2
  
3
  
4
 

约束七（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))

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
 

求解模型（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)

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
 

导出模型：

model.write('VRPTW.lp')

  

 

  
1
 

求解结果（部分）：

这样就完成了，感谢大家的阅读。

作者：夏旸，清华大学，工业工程系/深圳国际研究生院 (硕士在读)
刘兴禄，清华大学，清华伯克利深圳学院（博士在读）

完整Project文件请关注运小筹公众号，回复【Vrptw】获取。

文章来源: blog.csdn.net，作者：旸1111，版权归原作者所有，如需转载，请联系作者。

原文链接：blog.csdn.net/qq_35008055/article/details/110201571