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

此文转载自：https://blog.csdn.net/qq_35008055/article/details/110201571

手把手教你用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,jin V, i e 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 kin 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_{kin K} sum_{(i,j)in A}c_{ij}x_{ij} ag{1} mink∈K∑​(i,j)∈A∑​cij​xij​(1)
约束一让车辆驶出仓库(depot)：
∑ ( 0 , j ) ∈ A x 0 j k = 1 ,   ∀ k ∈ K (2) sum_{(0,j)in A}x_{0j}^k =1, forall kin K ag{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 kin K ,i in Vsetminus {s,t} ag{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 kin K ag{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 Vsetminus {s,t} ag{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 ag{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 Vsetminus {s,t} ag{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 ag{8} (i,j)∈A∑​xijk​qi​≤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}^kin {0,1} forall (i,j) in A,k in K \ s_i ge 0, forall i in Vsetminus {s,t} ag{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("

-----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】获取。
在这里插入图片描述

   

更多内容详见微信公众号：Python测试和开发

Python测试和开发

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/phyger/p/14054715.html