手把手教你用Gurobi求解一个数学模型_gurobi r语言 多元方程

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

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

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

对于Vrptw问题来说，数学模型主要由以下部分组成。首先我们定义一些相关参数，一个图可以表示为$G(V,A)$，其中 V = { 0 , 1 , . . . , n , n + 1 } V= \{ 0,1,...,n,n+1 \} 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之外的点)：
∑ ( 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：
$$\sum _{(i,n+1)\in A}{x}_{i,n+1}^k=1,\forall k\in K\text{\hspace{1em}(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_{ij}+ser{v}_i-M(1-x_{ij})\le s_j,\forall (i,j)\in A\text{\hspace{1em}(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)
约束七保证不违反车辆的载重约束：
$$\sum _{(i,j)\in A}{x}_{ij}^k{q}_i\le Q_k,\forall k\in K\text{\hspace{1em}(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   = [[]]   
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】获取。