混合线性规划-求解VRPTW问题（pyomo-glpk）

1、开始-准备工作

求解器--glpk

conda安装方法：

conda install -c conda-forge glpk

建模工具pyomo

conda安装方法：

conda install -c conda-forge pyomo.extras

建模文档：（这里的时间窗约束为软时间窗，早到可以等待不惩罚，晚到惩罚-R）

 若存在商品不能混装问题：

 简单解释一下模型中的常量

• 决策变量

  ​x_{ij}^k:binary，表示车辆k从i访问到j的弧是否被选择；

  ​s_{ki}：表示车辆k在顾客i处的到达时间；

  ​R_i：晚到具体时间；

• 常数项

  ​d_i：客户i的需求量

  lbtw_i​：最早到达时间

  ​ubtw_i：最晚到达时间

  ​M：常数项(较大数)，详细的可以搜索大M约束

代码实现：

1、数据处理-不返回（若是车辆需返回，在x_loc和y_loc中更改None为原点坐标，距离中注销掉if条件即可）

import numpy as np
# 这里原始数据不能发出来，就简单做几个模拟数据大家可以试运行一下；
# 数据处理-获取数据
d=[0,10,11,13,1,0]    # 需求量，包含起始点
x_loc = [0,1,2,3,4,None]   # 横坐标,最后模拟坐标不画出来，设为none
y_loc = [0,2,1,4,5,None]   # 纵坐标
t_ser = [0,0,0,0,0,0]   # 服务时间
lbtw = [0,0,0,0,0,0]    # 最早到达时间
ubtw = [100,100,120,150,140,100]  # 最晚到达时间

# 数据处理-处理为自己想要的数据格式
distance = [
    [0,1,2,3,4,0],
    [1,0,4,7,8,0],
    [2,4,0,6,2,0],
    [3,7,6,0,9,0],
    [4,8,2,9,0,0],
    [0,0,0,0,0,0]
]

v = 1 # m/s
C = np.array(distance)
t = C
M = 100
c_num = len(d) - 2   # 顾客数
v_num = 3  #卡车数，数量视情况而定
q = (0, 50)    # 卡车的容量

CM = [i for i in range (1,c_num+1)]  # 顾客列表
V = [i for i in range(1,v_num+1)]   # 卡车列表
N = [0]+CM+[c_num+1]   # 结点列表

## 若运输商品中有冲突的货物，此刻不注销
# Mod = a  # 货物A的客户数量      
# Col = b  # 货物B的客户数量
# Mo = [i for i in range(1, Mod+1)]
# Co = [i for i in range(Mod+1, Mod+Col+1)]
c1 = 5    # 惩罚系数

A = [(i,j) for i in N for j in N if i != j]
A_car = [(k,i,j) for k in V for i,j in A]
S_time = [(k,i) for k in V for i in N]

 2、建立模型

'''
此模型时间窗为软约束:
1、车辆早到需等待
2、不能出现晚到情况
3、不能处理较大的数据量数据
4、晚到时间不能太长<0.5
'''
from pyomo.environ import *
import pandas as pd
import numpy as np
import time


# 构建模型
model = ConcreteModel()

model.x = Var(A_car, within=Binary)
model.s = Var(N, within=NonNegativeReals)
model.R = Var(N, bounds=(0,0.5), within=NonNegativeReals)

# 目标函数
model.Obj = Objective(expr=sum(model.x[k,i,j] * C[i,j] * c1 for k,i,j in A_car) + sum(model.R[i] for i in N), sense=minimize)

# 约束1-每个客户都需要被服务到
def cons1(model, i):
    return sum(model.x[k,i,j] for k in V for j in N if j != i) == 1
model.cons_1 = Constraint(CM, rule=cons1)

# 约束2-不能超过总承重
def cons2_overweight(model, k):
    return sum(d[i]* sum(model.x[k,i,j] for j in N if j != i) for i in CM) <= q[1]
model.cons_2overweight = Constraint(V, rule=cons2_overweight)

# 约束3-车辆从起点出发
def cons3(model, k):
    return sum(model.x[k,0,j] for j in N if j != 0) == 1
model.cons_3 = Constraint(V, rule=cons3)

# 约束4-进出平衡
def cons4(model, k, h):
    return sum(model.x[k,i,h] for i in N if i != h) - sum(model.x[k,h,j] for j in N if j != h) == 0
model.cons_4 = Constraint(V, CM, rule=cons4)

# 约束5.1-回到原点
def cons5_1(model, k):
    return sum(model.x[k,i,c_num+1] for i in N if i != (c_num+1)) == 1    # 视情况需要，可以选择返回或者不返回，返回原点或者其他点，在开始数据处理的时候将返回点的坐标放上即可，不反悔的话让任何点到次点的距离为0即可
model.cons_5_1 = Constraint(V, rule=cons5_1)

# 约束5.2-禁止从虚拟节点到客户
def cons5_2(model, k):
    return sum(model.x[k,c_num+1,i] for i in N if i != (c_num+1)) == 0


# 约束6-发车数量小于总车辆
model.cons_6 = Constraint(expr=sum(model.x[k,0,j] for k in V for j in N if j != 0) <= v_num)

# 约束7-时间窗约束
model.cons_7 = ConstraintList()
for i,j in A:
    for k in V:
        model.cons_7.add(model.s[i] + t[i,j] - M* (1- model.x[k,i,j]) + R[i]<= model.s[j])
        model.cons_7.add(model.s[i] + t[i,j] + M* (1- model.x[k,i,j]) + R[i]>= model.s[j])
        model.cons_7.add(model.s[i] <= ubtw[i])
        model.cons_7.add(model.s[i] >= lbtw[i])

# 回到虚拟节点后禁止出来
def cons8(model, k):
    return (sum(model.x[k,c_num+1,i] for i in N if i != (c_num+1)) == 0)
model.cons_8 = Constraint(V, rule=cons8)

# 混装约束
# if Mod != 0 and Col != 0):    # 判断是否有混装商品，有的话进行约束
#     model.cons_8 = ConstraintList()
#     for i in range(Mod):
#         for j in range(Mod, Mod+Col):
#             for v in V:
#                 model.cons_8.add(sum(model.x[k,h,i] for h in self.N if h != i) + sum(model.x[k,h,j] for h in self.N if h != j) <= 1)
        
a = time.time()
# opt = SolverFactory('glpk').solve(model, options=glpk_opt)    # 当客户数量较多时可能短时间内给不出最优解，此时可以适当降低精度，减少求解时间
opt = SolverFactory('glpk').solve(model)

# 输出结果

for i in N:

    print(model.x[i]()) 

大致就是这样，注意仅适用于样本量较少的情况，若是样本量较大那么求解器的效率就不会太高，此时可能需要启发式的一些算法来解决问题；

update time --2022.8.22:

更新内容：晚到约束更改；数据处理，更新数据；

update time --2022.9.28:

更新内容：不能混装约束修改，优化之前存在的漏洞；修改时间窗约束，早到进行的等待（规定时间窗内）并进行惩罚；可行数据；增加约束，回到终点后禁止再出来；

update time --2022.10.17:

更新内容：禁止从虚拟节点到其他任何节点；