1. Toy模板网首页
  2. > Toy博客

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

10月前 作者：Better.C 分类：Toy博客 阅读(65) 违法举报

这篇具有很好参考价值的文章主要介绍了软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。


欢迎关注个人微信公众号：Python助力交通

1. 引言

优化求解器是解决复杂工程问题不可或缺的工具，可以帮助我们验证模型的正确性、理解决策变量的耦合关系、获取最优决策方案（合适规模条件下）。小编搜罗了网上关于各类常见（其实并不常见）的优化求解器介绍的帖子：

  • 优化求解器资源盘点
    • 干货 | 运筹学、数学规划、离散优化求解器大PK，总有一款适合你
    • 【学界】运筹学数学规划|离散优化求解器大搜罗
    • 除了Gurobi /SCIP，国内外还有哪些优化求解器？
    • 开源线性规划求解器（Linear Programming solver）LP_Solve和CLP的PK
    • 有哪些简便好用的解凸优化的工具箱或者包？
    • 开源建模框架+开源求解器 | 使用pyomo建模框架实现交通物流优化问题的求解

除了以上求解器外，还有一些针对特定问题而量身定制高效率求解器：

  • VRP问题求解器
    • 基于求解器的路径规划算法实现及性能分析
    • 智能优化算法工具包scikit-opt介绍及vrp问题求解评测（一）
    • Python主要智能优化算法库汇总

今天的主题是以CVRP和VRPTW问题为例，分享Python调用运筹优化求解器（Gurobi、COPT、SCIP）的教程。

2. 求解器介绍

（1）Gurobi
Gurobi是由美国 Gurobi Optimization 公司开发新一代大规模优化器，提供 C++, Java, Python, .Net, Matlab 和R等多种接口，支持求解以下模型：
（1）连续和混合整数线性问题
（2）凸目标或约束连续和混合整数二次问题
（3）非凸目标或约束连续和混合整数二次问题
（4）含有对数、指数、三角函数、高阶多项式目标或约束，以及任何形式的分段约束的非线性问题
（5）含有绝对值、最大值、最小值、逻辑与或非目标或约束的非线性问题

（2）COPT
杉数求解器COPT（Cardinal Optimizer）是杉数自主研发的针对大规模优化问题的高效数学规划求解器套件，也是支撑杉数端到端供应链平台的核心组件，是目前同时具备大规模混合整数规划、线性规划（单纯形法和内点法）、半定规划、（混合整数）二阶锥规划以及（混合整数）凸二次规划和（混合整数）凸二次约束规划问题求解能力的综合性能数学规划求解器，为企业应对高性能求解的需求提供了更多选择。COPT支持所有主流编程接口：C、C++、C#、Python、Julia、Java、AMPL、GAMS、Pyomo、PuLP、CVXPY。

（3）SCIP
SCIP起源于ZIB(Zuse Institute Berlin)，由Tobias Achterberg奠定整个框架，是目前用于混合整数规划（MIP）和混合整数非线性规划（MINLP）的最快的非商业解算器之一，它允许用户对求解过程进行全面控制，支持自定义搜索树中的各个模块，在分支限界(Branch and Bound)过程中添加变量等功能。SCIP支持Python Java AMPL GAMS MATLAB等编程语言。

3. 基础语言

3.1 创建模型

# gurobi
cvrp_model = Model('cvrp')
# copt
env = Envr()
cvrp_model = env.createModel('cvrp')
# scip
cvrp_model = Model('cvrp')

3.2 添加变量

# gurobi
x = cvrp_model.addVar(vtype=GRB.BINARY, name='x') # 单个变量
X = cvrp_model.addVars(N, N, K, vtype=GRB.BINARY, name='X[i,j,k]') # 多个变量
# copt
x = cvrp_model.addVars(vtype=COPT.BINARY, name='x') # 单个变量
X = cvrp_model.addVars(N, N, K, vtype=COPT.BINARY, nameprefix='X[i,j,k]') # 多个变量
# scip
X[i,j,k] = cvrp_model.addVar(vtype="B", name=f"x({i},{j},{k})") # 单个变量(貌似只能添加单个变量)

3.3 添加目标函数

# gurobi
cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K), GRB.MINIMIZE)
# copt
cvrp_model.setObjective(quicksum(X[i, j, k] * cost[i, j] for i in N for j in N for k in K), sense=COPT.MINIMIZE)
# scip
cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K),'minimize' )

3.4 添加约束


# gurobi
cvrp_model.addConstr() # 单个约束
cvrp_model.addConstrs() # 多个约束
# copt
cvrp_model.addConstr() # 单个约束
cvrp_model.addConstrs() # 多个约束
# scip
cvrp_model.addCons() # 单个约束
cvrp_model.addConss()  # 多个约束

3.5 设置参数

# gurobi
cvrp_model.setParam(GRB.Param.LogFile, './gurobi_r101-31.log') # 保存日志
cvrp_model.Params.TimeLimit = 1200 # 设置求解时间
# copt
cvrp_model.setLogFile('./copt_r101-31.log')  # 保存日志
cvrp_model.param.timelimit = 1200 # 设置求解时间
# scip
cvrp_model.setLogfile('./scip_r101-31.log') # 保存日志
cvrp_model.setRealParam('limits/time', 1200)  # 设置求解时间

3.6 求解

# gurobi
cvrp_model.optimize()
# copt
cvrp_model.solve()
# scip
cvrp_model.optimize()

4. 数学模型

4.1 CVRP数学模型

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

4.2 VRPTW数学模型

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

5. 完整代码

5.1 Python调用Gurobi求解CVRP

import xlsxwriter
import math
import pandas as pd
import matplotlib.pyplot as plt
from gurobipy import Model,quicksum,GRB
def read_input(filename):
    """
    :param filename: 数据文件路径
    :return:
    """
    df = pd.read_csv(filename)
    x_coord = { df['id'][i]:df['x_coord'][i] for i in range(df.shape[0]) }       # 节点横坐标
    y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) }       # 节点纵坐标
    demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) }         # 节点需求
    cost = {}
    N = df['id'].tolist()
    for f_n in N:
        for t_n in N:
            dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
            cost[f_n,t_n] = dist
    return N,cost,demand,x_coord,y_coord

def build_model(N,K,depot,CAP,cost,demand):
    """
    :param N: 网络节点集合
    :param K: 车队集合
    :param depot: 配送中心id
    :param CAP: 车辆容量
    :param cost: 网络弧费用集合
    :param demand: 网络节点需求集合
    :return:
    """
    cvrp_model = Model('cvrp')
    # 添加变量
    X = cvrp_model.addVars(N, N, K, vtype=GRB.BINARY, name='X[i,j,k]')
    Y = cvrp_model.addVars(N, K, vtype=GRB.BINARY, name='Y[i,k]')
    U = cvrp_model.addVars(N, K, vtype=GRB.INTEGER, name='U[i,k]')
    # 设置目标函数
    cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K), GRB.MINIMIZE)
    # 添加约束
    #  需求覆盖约束
    cvrp_model.addConstrs( quicksum(Y[i,k] for k in K) == 1 for i in N[1:] )
    #  车辆启用约束
    cvrp_model.addConstr( quicksum(Y[depot,k] for k in K) == len(K) )
    #  车辆流平衡约束
    cvrp_model.addConstrs( quicksum(X[i,j,k] for j in N ) == quicksum(X[j,i,k] for j in N ) for i in N for k in K )
    #  车辆路径限制
    cvrp_model.addConstrs( quicksum(X[i,j,k] for j in N) == Y[i,k] for i in N for k in K )
    #  车辆容量约束
    cvrp_model.addConstrs( quicksum(Y[i,k] * demand[i] for i in N) <= CAP for k in K )
    #  破圈约束
    cvrp_model.addConstrs( U[i,k] - U[j,k] + CAP * X[i,j,k] <= CAP - demand[i] for i in N[1:] for j in N[1:] for k in K )
    cvrp_model.addConstrs( U[i,k] <=  CAP for i in N[1:] for k in K)
    cvrp_model.addConstrs( U[i,k] >= demand[i] for i in N[1:] for k in K )
    return cvrp_model,X,Y,U

def draw_routes(route_list,x_coord,y_coord):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(x_coord[n])
            path_y.append(y_coord[n])
        plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
    plt.show()

def save_file(route_list,total_cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    row = 2
    for route in route_list:
        ws.write(row,0,route[0])
        route_str = [str(i) for i in route[1:]]
        ws.write(row,1,'-'.join(route_str))
        row += 1
    wb.close()

if __name__ == '__main__':
    filename = './datasets/CVRP/r101-31.csv'
    N, cost, demand, x_coord, y_coord = read_input(filename)
    depot = N[0]
    K = list(range(1,10))
    CAP = 80
    cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)
    cvrp_model.setParam(GRB.Param.LogFile, './gurobi_r101-31.log')
    cvrp_model.Params.TimeLimit = 1200
    cvrp_model.optimize()
    route_list = []
    for k in K:
        route = [depot]
        cur_node = depot
        cur_k = None
        for j in N[1:]:
            if X[depot, j, k].x > 0:
                cur_node = j
                cur_k = k
                route.append(j)
                N.remove(j)
                break
        if cur_k is None:
            continue
        while cur_node != depot:
            for j in N:
                if X[cur_node, j, cur_k].x > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        N.remove(j)
                    break
        route_list.append(route)
    print("最优路径距离:", cvrp_model.objVal)
    print("最优路径使用车辆数:", len(route_list))
    draw_routes(route_list, x_coord, y_coord)
    save_file(route_list, cvrp_model.objVal)

5.2 Python调用Gurobi求解VRPTW

import math
import pandas as pd
import matplotlib.pyplot as plt
import xlsxwriter
from gurobipy import GRB,Model,quicksum,tupledict,tuplelist

def read_input(filename):
    """
    :param filename: 数据文件路径
    :return:
    """
    N = [] #所有节点
    Q = {} #节点需求
    TT = {} #节点旅行时间
    ET = {} #节点最早开始服务时间
    LT = {} #节点最晚结束服务时间
    ST = {} #节点服务时间
    x_coord = {} # 节点横坐标
    y_coord = {} # 节点纵坐标
    Cost={}
    df = pd.read_csv(filename)
    for i in range(df.shape[0]):
        id = df['id'][i]
        N.append(id)
        x_coord[id] = df['x_coord'][i]
        y_coord[id] = df['y_coord'][i]
        Q[id] = df['demand'][i]
        ET[id] = df['start_time'][i]
        LT[id] = df['end_time'][i]
        ST[id] = df['service_time'][i]
    for f_n in N:
        for t_n in N:
            dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
            Cost[f_n,t_n] = dist
            TT[f_n,t_n] = dist/1 # 假设速度为1
    return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coord

def build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K):
    """
    :param N: 所有节点集合，其中N[0]为车场
    :param Q: 节点需求集合
    :param TT: 旅行时间
    :param ET: 节点最早开始服务时间
    :param LT：节点最晚结束服务时间
    :param ST: 节点服务时间
    :param CAP: 车辆容量
    :param Cost: 旅行费用
    :param K: 车队
    :return:
    """
    C = N[1:] #需求节点
    M=10**5
    depot = N[0]
    # 创建模型
    vrptw_model=Model()
    # 添加变量
    X = vrptw_model.addVars(N,N,K,vtype=GRB.BINARY,name='X(i,j,k)')
    T = vrptw_model.addVars( N,K,vtype=GRB.CONTINUOUS,lb=0,name='T[i,k]')
    # 设置目标函数
    z1 = quicksum( Cost[i,j]*X[i,j,k] for i in N for j in N for k in K if i!=j)
    vrptw_model.setObjective(z1,GRB.MINIMIZE)
    # 车辆起点约束
    vrptw_model.addConstrs(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
    # 车辆路径连续约束
    vrptw_model.addConstrs( quicksum(X[i,j,k] for j in N if j!=i)==quicksum(X[j,i,k] for j in N if j!=i) for i in C for k in K)
    # 车辆终点约束
    vrptw_model.addConstrs(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
    # 需求服务约束
    vrptw_model.addConstrs( quicksum(X[i,j,k] for k in K for j in N if j!=i)==1 for i in C)
    # 车辆容量约束
    vrptw_model.addConstrs( quicksum(Q[i]*X[i,j,k] for i in C for j in N if i!=j)<=CAP for k in K )
    # 时间窗约束
    vrptw_model.addConstrs( T[i,k]+ST[i]+TT[i,j]-(1-X[i,j,k])*M<=T[j,k] for i in C for j in C for k in K if i!=j )
    vrptw_model.addConstrs( T[i,k] >= ET[i] for i in N for k in K)
    vrptw_model.addConstrs( T[i,k] <= LT[i] for i in N for k in K)
    return vrptw_model,X,T
def draw_routes(route_list,x_coord,y_coord):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(x_coord[n])
            path_y.append(y_coord[n])
        plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
    plt.show()
def save_file(route_list,route_timetable,total_cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    ws.write(1,2,'时刻')
    row = 2
    for id,route in enumerate(route_list):
        ws.write(row,0,route[0])
        route_str = [str(i) for i in route[1:]]
        ws.write(row,1,'-'.join(route_str))
        timetable_str = [str(i) for i in route_timetable[id]]
        ws.write(row,2,'-'.join(timetable_str))
        row += 1
    wb.close()
if __name__=='__main__':
    filename = './data/r101.csv'
    N, Q, TT, ET, LT, ST, Cost,x_coord,y_coord = read_input(filename)
    depot = N[0]
    K = list(range(1, 30))
    CAP = 80
    vrptw_model,X,T = build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K)
    vrptw_model.setParam(GRB.Param.LogFile, './log/gurobi_r101.log')
    vrptw_model.Params.TimeLimit = 1500
    vrptw_model.optimize()
    route_list = []
    route_timetable = []
    # 提取车辆路径
    for k in K:
        route = [depot]
        cur_node = depot
        cur_k = None
        for j in N[1:]:
            if X[depot, j, k].x > 0:
                cur_node = j
                cur_k = k
                route.append(j)
                N.remove(j)
                break
        if cur_k is None:
            continue
        while cur_node != depot:
            for j in N:
                if X[cur_node, j, cur_k].x > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        N.remove(j)
                    break
        route.insert(0,k)
        route_list.append(route)
    # 提取车辆的时刻点
    for route in route_list:
        k = route[0]
        timetable = []
        for i in route[1:]:
            timetable.append(T[i,k].x)
        route_timetable.append(timetable)

    print("最优路径距离:", vrptw_model.objVal)
    print("最优路径使用车辆数:", len(route_list))
    draw_routes(route_list, x_coord, y_coord)
    save_file(route_list, route_timetable, vrptw_model.objVal)

5.3 Python调用COPT求解CVRP


import xlsxwriter
import math
import pandas as pd
import matplotlib.pyplot as plt
from coptpy import *
def read_input(filename):
    """
    :param filename: 数据文件路径
    :return:
    """
    df = pd.read_csv(filename)
    x_coord = { df['id'][i]:df['x_coord'][i] for i in range(df.shape[0]) }       # 节点横坐标
    y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) }       # 节点纵坐标
    demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) }         # 节点需求
    cost = {}
    N = df['id'].tolist()
    for f_n in N:
        for t_n in N:
            dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
            cost[f_n,t_n] = dist
    return N,cost,demand,x_coord,y_coord
def build_model(N,K,depot,CAP,cost,demand):
    # 创建求解环境
    env = Envr()
    # 创建求解模型
    cvrp_model = env.createModel('cvrp')
    # 添加决策变量
    X = cvrp_model.addVars(N, N, K, vtype=COPT.BINARY, nameprefix='X[i,j,k]')
    Y = cvrp_model.addVars(N, K, vtype=COPT.BINARY, nameprefix='Y[i,k]')
    U = cvrp_model.addVars(N, K, vtype=COPT.INTEGER, nameprefix='U[i,k]')
    # 设置目标函数
    cvrp_model.setObjective(quicksum(X[i, j, k] * cost[i, j] for i in N for j in N for k in K), sense=COPT.MINIMIZE)
    # 添加约束
    #  需求覆盖约束
    cvrp_model.addConstrs(quicksum(Y[i, k] for k in K) == 1 for i in N[1:])
    #  车辆启用约束
    cvrp_model.addConstr(quicksum(Y[depot, k] for k in K) == len(K))
    #  车辆流平衡约束
    cvrp_model.addConstrs(quicksum(X[i, j, k] for j in N) == quicksum(X[j, i, k] for j in N) for i in N for k in K)
    #  车辆路径限制
    cvrp_model.addConstrs(quicksum(X[i, j, k] for j in N) == Y[i, k] for i in N for k in K)
    #  车辆容量约束
    cvrp_model.addConstrs(quicksum(Y[i, k] * demand[i] for i in N) <= CAP for k in K)
    #  破圈约束
    cvrp_model.addConstrs(
        U[i, k] - U[j, k] + CAP * X[i, j, k] <= CAP - demand[i] for i in N[1:] for j in N[1:] for k in K)
    cvrp_model.addConstrs(U[i, k] <= CAP for i in N[1:] for k in K)
    cvrp_model.addConstrs(U[i, k] >= demand[i] for i in N[1:] for k in K)
    return cvrp_model, X, Y, U
def draw_routes(route_list,x_coord,y_coord):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(x_coord[n])
            path_y.append(y_coord[n])
        plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
    plt.show()
def save_file(route_list,total_cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    row = 2
    for route in route_list:
        ws.write(row,0,route[0])
        route_str = [str(i) for i in route[1:]]
        ws.write(row,1,'-'.join(route_str))
        row += 1
    wb.close()
if __name__ == '__main__':
    filename = './datasets/CVRP/r101-31.csv'
    N, cost, demand, x_coord,y_coord = read_input(filename)
    depot = N[0]
    K = list(range(1,10))
    CAP = 80
    cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)
    cvrp_model.param.timelimit = 1200
    cvrp_model.setLogFile('./copt_r101-31.log')
    cvrp_model.solve()
    route_list = []
    for k in K:
        route = [depot]
        cur_node = depot
        cur_k = None
        for j in N[1:]:
            if X[depot, j, k].x > 0:
                cur_node = j
                cur_k = k
                route.append(j)
                N.remove(j)
                break
        if cur_k is None:
            continue
        while cur_node != depot:
            for j in N:
                if X[cur_node, j, cur_k].x > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        N.remove(j)
                    break
        route_list.append(route)
    print("最优路径距离:", cvrp_model.objval)
    print("最优路径使用车辆数:", len(route_list))
    draw_routes(route_list, x_coord, y_coord)
    save_file(route_list, cvrp_model.objval)

5.4 Python调用COPT求解VRPTW


import math
import pandas as pd
import matplotlib.pyplot as plt
import xlsxwriter
from coptpy import Envr,COPT,quicksum

def read_input(filename):
    """
    :param filename: 数据文件路径
    :return:
    """
    N = [] #所有节点
    Q = {} #节点需求
    TT = {} #节点旅行时间
    ET = {} #节点最早开始服务时间
    LT = {} #节点最晚结束服务时间
    ST = {} #节点服务时间
    x_coord = {} # 节点横坐标
    y_coord = {} # 节点纵坐标
    Cost={}
    df = pd.read_csv(filename)
    for i in range(df.shape[0]):
        id = df['id'][i]
        N.append(id)
        x_coord[id] = df['x_coord'][i]
        y_coord[id] = df['y_coord'][i]
        Q[id] = df['demand'][i]
        ET[id] = df['start_time'][i]
        LT[id] = df['end_time'][i]
        ST[id] = df['service_time'][i]
    for f_n in N:
        for t_n in N:
            dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
            Cost[f_n,t_n] = dist
            TT[f_n,t_n] = dist/1 # 假设速度为1
    return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coord
def build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K):
    """
    :param N: 所有节点集合，其中N[0]为车场
    :param Q: 节点需求集合
    :param TT: 旅行时间
    :param ET: 节点最早开始服务时间
    :param LT：节点最晚结束服务时间
    :param ST: 节点服务时间
    :param CAP: 车辆容量
    :param Cost: 旅行费用
    :param K: 车队
    :return:
    """
    C = N[1:] #需求节点
    M=10**5
    depot = N[0]
    # 创建求解环境
    env = Envr()
    # 创建求解模型
    vrptw_model = env.createModel('vrptw')
    # 添加决策变量
    X = vrptw_model.addVars(N, N, K, vtype=COPT.BINARY, nameprefix='X[i,j,k]')
    T = vrptw_model.addVars( N,K,vtype=COPT.CONTINUOUS,lb=0,nameprefix='T[i,k]')
    # 设置目标函数
    vrptw_model.setObjective(quicksum(X[i, j, k] * Cost[i, j] for i in N for j in N for k in K), sense=COPT.MINIMIZE)
    # 车辆起点约束
    vrptw_model.addConstrs(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
    # 车辆路径连续约束
    vrptw_model.addConstrs(
        quicksum(X[i, j, k] for j in N if j != i) == quicksum(X[j, i, k] for j in N if j != i) for i in C for k in K)
    # 车辆终点约束
    vrptw_model.addConstrs(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
    # 需求服务约束
    vrptw_model.addConstrs(quicksum(X[i, j, k] for k in K for j in N if j != i) == 1 for i in C)
    # 车辆容量约束
    vrptw_model.addConstrs(quicksum(Q[i] * X[i, j, k] for i in C for j in N if i != j) <= CAP for k in K)
    # 时间窗约束
    vrptw_model.addConstrs(
        T[i, k] + ST[i] + TT[i, j] - (1 - X[i, j, k]) * M <= T[j, k] for i in C for j in C for k in K if i != j)
    vrptw_model.addConstrs(T[i, k] >= ET[i] for i in N for k in K)
    vrptw_model.addConstrs(T[i, k] <= LT[i] for i in N for k in K)
    return vrptw_model, X, T
def draw_routes(route_list,x_coord,y_coord):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(x_coord[n])
            path_y.append(y_coord[n])
        plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
    plt.show()
def save_file(route_list,route_timetable,total_cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    ws.write(1,2,'时刻')
    row = 2
    for id,route in enumerate(route_list):
        ws.write(row,0,route[0])
        route_str = [str(i) for i in route[1:]]
        ws.write(row,1,'-'.join(route_str))
        timetable_str = [str(i) for i in route_timetable[id]]
        ws.write(row,2,'-'.join(timetable_str))
        row += 1
    wb.close()
if __name__=='__main__':
    filename = './data/r101.csv'
    N, Q, TT, ET, LT, ST, Cost,x_coord,y_coord = read_input(filename)
    depot = N[0]
    K = list(range(1, 30))
    CAP = 80
    vrptw_model,X,T = build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K)
    vrptw_model.setLogFile('./log/copt_c201.log')
    vrptw_model.param.timelimit = 1500
    vrptw_model.solve()
    route_list = []
    route_timetable = []
    # 提取车辆路径
    for k in K:
        route = [depot]
        cur_node = depot
        cur_k = None
        for j in N[1:]:
            if X[depot, j, k].x > 0:
                cur_node = j
                cur_k = k
                route.append(j)
                N.remove(j)
                break
        if cur_k is None:
            continue
        while cur_node != depot:
            for j in N:
                if X[cur_node, j, cur_k].x > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        N.remove(j)
                    break
        route.insert(0,k)
        route_list.append(route)
    # 提取车辆的时刻点
    for route in route_list:
        k = route[0]
        timetable = []
        for i in route[1:]:
            timetable.append(T[i,k].x)
        route_timetable.append(timetable)

    print("最优路径距离:", vrptw_model.objval)
    print("最优路径使用车辆数:", len(route_list))
    draw_routes(route_list, x_coord, y_coord)
    save_file(route_list, route_timetable, vrptw_model.objval)

5.5 Python调用SCIP求解CVRP


import xlsxwriter
import math
import pandas as pd
import matplotlib.pyplot as plt
from pyscipopt import Model, quicksum, multidict

def read_input(filename):
    """
    :param filename: 数据文件路径
    :return:
    """
    df = pd.read_csv(filename)
    x_coord = { df['id'][i]:df['x_coord'][i] for i in range(df.shape[0]) }       # 节点横坐标
    y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) }       # 节点纵坐标
    demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) }         # 节点需求
    cost = {}
    N = df['id'].tolist()
    for f_n in N:
        for t_n in N:
            dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
            cost[f_n,t_n] = dist
    return N,cost,demand,x_coord,y_coord

def build_model(N,K,depot,CAP,cost,demand):
    """
    :param N: 网络节点集合
    :param K: 车队集合
    :param depot: 配送中心id
    :param CAP: 车辆容量
    :param cost: 网络弧费用集合
    :param demand: 网络节点需求集合
    :return:
    """
    cvrp_model = Model('cvrp')
    # 添加变量
    X = {}
    Y = {}
    U = {}
    # 创建变量 X[i,j,k]
    for i in N:
        for j in N:
            for k in K:
                X[i,j,k] = cvrp_model.addVar(vtype="B", name=f"x({i},{j},{k})")
    # 创建变量 Y[i,k]
    for i in N:
        for k in K:
            Y[i,k] = cvrp_model.addVar(vtype="B", name=f"y({i},{k})")
    #  创建辅助变量U[i,k]
    for i in N:
        for k in K:
            U[i,k] = cvrp_model.addVar(vtype="I", name=f"u({i},{k})")
    # 设置目标函数
    cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K),'minimize' )
    # 添加约束
    #  需求覆盖约束
    cvrp_model.addConss( quicksum(Y[i,k] for k in K) == 1 for i in N[1:] )
    #  车辆启用约束
    cvrp_model.addCons( quicksum(Y[depot,k] for k in K) == len(K) )
    #  车辆流平衡约束
    cvrp_model.addConss( quicksum(X[i,j,k] for j in N ) == quicksum(X[j,i,k] for j in N ) for i in N for k in K )
    #  车辆路径限制
    cvrp_model.addConss( quicksum(X[i,j,k] for j in N) == Y[i,k] for i in N for k in K )
    #  车辆容量约束
    cvrp_model.addConss( quicksum(Y[i,k] * demand[i] for i in N) <= CAP for k in K )
    #  破圈约束
    cvrp_model.addConss( U[i,k] - U[j,k] + CAP * X[i,j,k] <= CAP - demand[i] for i in N[1:] for j in N[1:] for k in K )
    cvrp_model.addConss( U[i,k] <=  CAP for i in N[1:] for k in K)
    cvrp_model.addConss( U[i,k] >= demand[i] for i in N[1:] for k in K )
    return cvrp_model, X, Y, U

def draw_routes(route_list,x_coord,y_coord):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(x_coord[n])
            path_y.append(y_coord[n])
        plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
    plt.show()
def save_file(route_list,total_cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    row = 2
    for route in route_list:
        ws.write(row,0,route[0])
        route_str = [str(i) for i in route[1:]]
        ws.write(row,1,'-'.join(route_str))
        row += 1
    wb.close()
if __name__ == '__main__':
    filename = './datasets/CVRP/r101-31.csv'
    N, cost, demand, x_coord, y_coord = read_input(filename)
    depot = N[0]
    K = list(range(1,10))
    CAP = 80
    cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)
    # cvrp_model.setRealParam('limits/gap', 0.1)  # 求解gap为10%
    cvrp_model.setRealParam('limits/time', 1200)  # 求解时间为1800s
    cvrp_model.setLogfile('./scip_r101-31.log')
    cvrp_model.optimize()
    route_list = []
    for k in K:
        route = [depot]
        cur_node = depot
        cur_k = None
        for j in N[1:]:
            if cvrp_model.getVal(X[depot, j, k]) > 0:
                cur_node = j
                cur_k = k
                route.append(j)
                N.remove(j)
                break
        if cur_k is None:
            continue
        while cur_node != depot:
            for j in N:
                if cvrp_model.getVal(X[cur_node, j, cur_k]) > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        N.remove(j)
                    break
        route_list.append(route)
    print("最优路径距离:", cvrp_model.getObjVal())
    print("最优路径使用车辆数:", len(route_list))
    draw_routes(route_list, x_coord, y_coord)
    save_file(route_list, cvrp_model.getObjVal())

5.6 Python调用SCIP求解VRPTW


import math
import pandas as pd
import matplotlib.pyplot as plt
import xlsxwriter
from pyscipopt import Model, quicksum

def read_input(filename):
    """
    :param filename: 数据文件路径
    :return:
    """
    N = [] #所有节点
    Q = {} #节点需求
    TT = {} #节点旅行时间
    ET = {} #节点最早开始服务时间
    LT = {} #节点最晚结束服务时间
    ST = {} #节点服务时间
    x_coord = {} # 节点横坐标
    y_coord = {} # 节点纵坐标
    Cost={}
    df = pd.read_csv(filename)
    for i in range(df.shape[0]):
        id = df['id'][i]
        N.append(id)
        x_coord[id] = df['x_coord'][i]
        y_coord[id] = df['y_coord'][i]
        Q[id] = df['demand'][i]
        ET[id] = df['start_time'][i]
        LT[id] = df['end_time'][i]
        ST[id] = df['service_time'][i]
    for f_n in N:
        for t_n in N:
            dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
            Cost[f_n,t_n] = dist
            TT[f_n,t_n] = dist/1 # 假设速度为1
    return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coord
def build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K):
    """
    :param N: 所有节点集合，其中N[0]为车场
    :param Q: 节点需求集合
    :param TT: 旅行时间
    :param ET: 节点最早开始服务时间
    :param LT：节点最晚结束服务时间
    :param ST: 节点服务时间
    :param CAP: 车辆容量
    :param Cost: 旅行费用
    :param K: 车队
    :return:
    """
    C = N[1:] #需求节点
    M=10**5
    depot = N[0]
    # 创建模型
    vrptw_model=Model()
    X = {}
    T = {}
    # 创建变量 X[i,j,k]
    for i in N:
        for j in N:
            for k in K:
                X[i, j, k] = vrptw_model.addVar(vtype="B", name=f"x({i},{j},{k})")
    # 创建变量 T[i,k]
    for i in N:
        for k in K:
            T[i, k] = vrptw_model.addVar(vtype="I", name=f"t({i},{k})")
    # 设置目标函数
    vrptw_model.setObjective(quicksum(X[i, j, k] * Cost[i, j] for i in N for j in N for k in K), 'minimize')
    # 车辆起点约束
    # 车辆起点约束
    vrptw_model.addConss(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
    # 车辆路径连续约束
    vrptw_model.addConss(
        quicksum(X[i, j, k] for j in N if j != i) == quicksum(X[j, i, k] for j in N if j != i) for i in C for k in K)
    # 车辆终点约束
    vrptw_model.addConss(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
    # 需求服务约束
    vrptw_model.addConss(quicksum(X[i, j, k] for k in K for j in N if j != i) == 1 for i in C)
    # 车辆容量约束
    vrptw_model.addConss(quicksum(Q[i] * X[i, j, k] for i in C for j in N if i != j) <= CAP for k in K)
    # 时间窗约束
    vrptw_model.addConss(
        T[i, k] + ST[i] + TT[i, j] - (1 - X[i, j, k]) * M <= T[j, k] for i in C for j in C for k in K if i != j)
    vrptw_model.addConss(T[i, k] >= ET[i] for i in N for k in K)
    vrptw_model.addConss(T[i, k] <= LT[i] for i in N for k in K)
    return vrptw_model, X, T
def draw_routes(route_list,x_coord,y_coord):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(x_coord[n])
            path_y.append(y_coord[n])
        plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
    plt.show()
def save_file(route_list,route_timetable,total_cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    ws.write(1,2,'时刻')
    row = 2
    for id,route in enumerate(route_list):
        ws.write(row,0,route[0])
        route_str = [str(i) for i in route[1:]]
        ws.write(row,1,'-'.join(route_str))
        timetable_str = [str(i) for i in route_timetable[id]]
        ws.write(row,2,'-'.join(timetable_str))
        row += 1
    wb.close()
if __name__=='__main__':
    filename = './data/c201.csv'
    N, Q, TT, ET, LT, ST, Cost,x_coord,y_coord = read_input(filename)
    depot = N[0]
    K = list(range(1, 30))
    CAP = 80
    vrptw_model,X,T = build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K)
    vrptw_model.setLogfile('./log/scip_c101.log')
    vrptw_model.setRealParam('limits/time', 2000)  # 求解时间为1800s
    vrptw_model.optimize()
    route_list = []
    route_timetable = []
    # 提取车辆路径
    for k in K:
        route = [depot]
        cur_node = depot
        cur_k = None
        for j in N[1:]:
            if vrptw_model.getVal(X[depot, j, k]) > 0:
                cur_node = j
                cur_k = k
                route.append(j)
                N.remove(j)
                break
        if cur_k is None:
            continue
        while cur_node != depot:
            for j in N:
                if vrptw_model.getVal(X[cur_node, j, cur_k]) > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        N.remove(j)
                    break
        route.insert(0,k)
        route_list.append(route)
    # 提取车辆的时刻点
    for route in route_list:
        k = route[0]
        timetable = []
        for i in route[1:]:
            timetable.append(vrptw_model.getVal(T[i,k]))
        route_timetable.append(timetable)

    print("最优路径距离:", vrptw_model.getObjVal())
    print("最优路径使用车辆数:", len(route_list))
    draw_routes(route_list, x_coord, y_coord)
    save_file(route_list, route_timetable, vrptw_model.getObjVal())

7. 测试案例

本期从solomon数据集中选择了c101，c201，cr101三个算例用于测评。
CVRP问题数据结构如下（保留30个节点）：
软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例
VRPTW问题数据结构如下（保留100个节点）：
软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

8. 测试参数

软件版本号：

Gurbi：9.5.1
SCIP：8.0.0
COPT：5.0.1

CVRP相关参数设置：

车辆容量：80
车队规模：10
客户节点：30
求解时间：1200s

VRPTW相关参数设置：

车辆容量：80
车队规模：30
客户节点：100
求解时间：1500s

9. 测试结果

9.1 CVRP求解结果汇总

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例
对于CVRP问题，在相同的求解时间下，Gurobi在求解质量方面有最好的表现，COPT比SCIP略差。

9.2 VRPTW求解结果汇总

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

对于VRPTW问题，在相同的求解时间下，Gurobi可以找到不错的优化解，而COPT和SCIP并未找到初始解。

9.3 上界下降曲线对比（以CVRP为例）

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例
软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例
软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例
从上界下降曲线可以看出，Gurobi找到的初始可行解的质量最好，且收敛性最好，而COPT和SCIP的初始解质量相差不大，收敛性不如Gurobi。

9.5 车辆路径可视化（以CVRP为例）

c101-31（Gurobi）

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

c201-31（Gurobi）

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

r101-31（Gurobi）

软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例

10. 小节

根据以上测试结果，Gurobi的求解性能远优于SCIP和COPT，作为后起之秀COPT已经接近SCIP水平（最新版的COPT尚未测试，也许性能有更大提升）。
虽然求解器可以直接对数学模型进行求解，但仅在小规模案例上具有较好的性能，对于大规模问题仍需要设计更高效的算法。
小编后续将继续分享VRP问题系列求解算法，并将其求解效果与上述求解器进行比较。文章来源地址https://www.toymoban.com/news/detail-498736.html

到了这里，关于软件工具 | Python调用运筹优化求解器（一）：以CVRP&VRPTW为例的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处： 如若内容造成侵权/违法违规/事实不符，请点击违法举报进行投诉反馈，一经查实，立即删除！

分享到：
领支付宝红包 赞助服务器费用

相关文章

  • 基于GA-PSO遗传粒子群混合优化算法的VRPTW问题求解matlab仿真

    基于GA-PSO遗传粒子群混合优化算法的VRPTW问题求解matlab仿真

    目录 1.程序功能描述 2.测试软件版本以及运行结果展示 3.核心程序 4.本算法原理 4.1 遗传算法（GA）基本原理 4.2 粒子群优化（PSO）基本原理 4.3 算法优化策略 5.完整程序        VRPTW是车辆路径问题（VRP）的一个扩展，它在基本的车辆路径问题上增加了对客户服务时间窗的考虑

    2024年02月02日
    浏览(77)
  • 【运筹优化】ALNS自适应大领域搜索算法求解TSP问题 + Java代码实现

    【运筹优化】ALNS自适应大领域搜索算法求解TSP问题 + Java代码实现

    旅行推销员问题(TSP)提出以下问题：“给定 n n n 个城市的列表，其中有一个起始城市，以及每对城市之间的距离，访问每个城市一次并返回起始城市的最短可能路线是什么？”。 这又是一个重要的NP-hard组合优化，特别是在运筹学和理论计算机科学领域。这个问题最早是在1

    2024年02月13日
    浏览(51)
  • 【运筹优化】子集和问题（Subset Sum Problems , SSP）介绍 + 动态规划求解 + Java代码实现

    【运筹优化】子集和问题（Subset Sum Problems , SSP）介绍 + 动态规划求解 + Java代码实现

    子集和问题（Subset Sum Problems , SSP），它是复杂性理论中最重要的问题之一。 SSP会给定一组整数 a 1 , a 2 , . . . . , a n a_1,a_2,....,a_n a 1 ​ , a 2 ​ , .... , a n ​ ，最多 n n n 个整数，我们需要判断是否存在一个非空子集，使得子集的总和为 M M M 整数？如果存在则需要输出该子集。

    2024年01月17日
    浏览(50)
  • 运筹学—例题求解

    运筹学—例题求解

    作答如下：      图解法验证：  由图可得在点x1=20,x2=24取到最大值 Z =4080； 作答如下： 解： (1)设 xij 为从产地Ai运往销地Bj的运输量，得到下列运输量表设 xij 为从产地Ai运往销地Bj的运输量，得到下列运输量表   B1 B2 B3 产量 A1 x 11 x 12 x 13 200 A2 x 21 x 22 x 23 230 销量 100 150 180

    2024年02月04日
    浏览(59)
  • 运筹系列82：使用动态规划求解TSP问题

    运筹系列82：使用动态规划求解TSP问题

    定义 c ( s , k ) c(s,k) c ( s , k ) 为当前在 k k k ，待访问点的集合 s s s ，最后返回城市0的最短路径，那么Bellman方程为： c ( s , k ) = min ⁡ i ∈ s { c ( s − { i } , i ) + d i , k } c(s,k)=min_{i in s}{c(s-{i},i)+d_{i,k}} c ( s , k ) = min i ∈ s ​ { c ( s − { i } , i ) + d i , k ​ } 为了使用方便，这里

    2024年02月06日
    浏览(44)
  • 运筹系列87：julia求解随机动态规划问题入门

    运筹系列87：julia求解随机动态规划问题入门

    随机动态规划问题的特点是： 有多个阶段，每个阶段的随机性互不相关，且有有限个实现值 (finite realizations) 具有马尔可夫性质，即每个阶段只受上一个阶段影响，可以用状态转移方程来描述阶段与阶段之间的变化过程。 我们使用julia的SDDP算法包来求解随机动态规划问题。

    2024年01月16日
    浏览(43)

  • 带约束条件的运筹规划问题求解（模拟退火算法实现）

    超级简单的模拟退火算法实现ε٩(๑ ₃ )۶з搭配最简单的线性规划模型进行讲解！但是如果需要的话，可以直接修改程序求解非线性问题哦(´つヮ⊂︎) [max,f(x)=10x_1+9x_2] (s.t.) [6x_1+5x_2leq{60}tag{1}] [10x_1+20x_2leq{150}tag{2}] [0leq{x_1}leq{8}tag{3}] [0leq{x_2}leq{8}tag{4}] 对约束

    2023年04月18日
    浏览(46)
  • 运筹说 第67期 | 动态规划模型的建立与求解

    运筹说 第67期 | 动态规划模型的建立与求解

    通过前一期的学习，我们已经学会了动态规划的基本概念和基本原理。本期小编带大家学习动态规划模型的建立与求解。 动态规划模型的建立 一   概述 建立动态规划的模型，就是分析问题并建立问题的动态规划基本方程。 成功地应用动态规划方法的关键，在于识别问题的

    2024年01月16日
    浏览(36)
  • 数学建模1：lingo软件求解优化模型

    数学建模1：lingo软件求解优化模型

    本次数学建模学习笔记系列，以代码学习为主，附带建模及论文亮点记录 由于队友为两位经济学小伙伴，因此以大数据类型题目为主要学习方向 注：论文代码资料来源网络 1、结构清晰（后附该论文前两问的目录结构） 2、lingo求解优化模型，涉及函数循环与求和 3、表格很好

    2024年02月08日
    浏览(65)
  • python求解带约束的优化问题

    python求解带约束的优化问题

    带约束的优化问题可被定义为： 在python中，可以使用 scipy 的 optimize 包进行求解，具体求解函数为 linprog ，下面举例说明求解方法： 假设问题被定义为： 首先，求解最大值问题，我们可以通过取负转换为求解最小值问题，包括不等式约束也是如此，那么该问题的python求解代

    2024年02月08日
    浏览(43)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包，优惠每天领

二维码1

领取红包

二维码2

领红包