软件工具 | 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数学模型

4.2 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个节点）：

VRPTW问题数据结构如下（保留100个节点）：

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求解结果汇总

对于CVRP问题，在相同的求解时间下，Gurobi在求解质量方面有最好的表现，COPT比SCIP略差。

9.2 VRPTW求解结果汇总

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

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

从上界下降曲线可以看出，Gurobi找到的初始可行解的质量最好，且收敛性最好，而COPT和SCIP的初始解质量相差不大，收敛性不如Gurobi。

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

c101-31（Gurobi）

c201-31（Gurobi）

r101-31（Gurobi）

10. 小节

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