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

手把手教你用Gurobi求解一个数学模型
- 带时间窗的车辆路径规划问题(Vrptw)
- python调用Gurobi求解Vrptw
- - 首先我们定义一下需要用到的参数：
  - 定义一个读取数据的函数，并对节点之间的距离进行计算：
  - 读取数据，并定义一些参数：
  - 调用gurobi进行模型的建立与求解：
  - - 根据式(9)定义决策变量，并加入模型当中：
    - 根据式(1)定义目标函数，并加入模型当中：
    - 根据式(2)~(8)定义决策变量，并加入模型当中：

手把手教你用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}为图中所有点的集合，AAA为图中所有弧的集合，有(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)的单位运输费用为cijc_{ij}cij，运输时间为tijt_{ij}tij，每个客户点的需求为qijq_{ij}qij，可服务的时间窗为[ei,li][e_i,l_i][ei,li]，服务时长为serviserv_iservi。令车辆的集合为KKK，每辆车的最大载重为Qk,∀k∈KQ_k,\forall k\in KQk,∀k∈K。决策变量为xijkx_{ij}^kxijk，代表第kkk辆车是否服务了弧(i,j)(i,j)(i,j)；sis_isi为客户点iii开始被服务的时间。

接下来构建数学模型。
目标函数为最小化运输成本：
min∑k∈K∑(i,j)∈Acijxij(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)∈Ax0jk=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)∈Axijk−∑(j,i)∈Axjik=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)∈Axi,n+1k=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)∈Axi,jk=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)
约束五保证被服务的相邻节点开始服务时间的大小关系(去回路)：
si+tij+servi−M(1−xij)≤sj,∀(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)
约束六保证不违反客户的时间窗：
ei≤si≤li,∀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)∈Axijkqi≤Qk,∀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)
最后是决策变量的取值约束：
xijk∈{0,1}∀(i,j)∈A,k∈Ksi≥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 = 0nodeNum     = 0vehicleNum  = 0capacity    = 0cor_X       = []cor_Y       = []demand      = []serviceTime = []readyTime   = []dueTime     = []disMatrix   = [[]]

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

# function to read data from .txt files
def readData(data, path, customerNum):data.customerNum = customerNumdata.nodeNum = customerNum + 2f = open(path, 'r')lines = f.readlines()count = 0# read the infofor line in lines:count = count + 1if (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 matrixdata.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]) ** 2data.disMatrix[i][j] = math.sqrt(temp)temp = 0return 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】获取。