1. 引言

優(yōu)化求解器是解決復雜工程問題不可或缺的工具，可以幫助我們驗證模型的正確性、理解決策變量的耦合關系、獲取最優(yōu)決策方案（合適規(guī)模條件下）。小編搜羅了網上關于各類常見（其實并不常見）的優(yōu)化求解器介紹的帖子：

• 優(yōu)化求解器資源盤點
  • 干貨 | 運籌學、數學規(guī)劃、離散優(yōu)化求解器大PK，總有一款適合你
  • 【學界】運籌學數學規(guī)劃|離散優(yōu)化求解器大搜羅
  • 除了Gurobi /SCIP，國內外還有哪些優(yōu)化求解器？
  • 開源線性規(guī)劃求解器（Linear Programming solver）LP_Solve和CLP的PK
  • 有哪些簡便好用的解凸優(yōu)化的工具箱或者包？
  • 開源建?？蚣?開源求解器 | 使用pyomo建?？蚣軐崿F(xiàn)交通物流優(yōu)化問題的求解

除了以上求解器外，還有一些針對特定問題而量身定制高效率求解器：

• VRP問題求解器
  • 基于求解器的路徑規(guī)劃算法實現(xiàn)及性能分析
  • 智能優(yōu)化算法工具包scikit-opt介紹及vrp問題求解評測（一）
  • Python主要智能優(yōu)化算法庫匯總

今天的主題是以CVRP和VRPTW問題為例，分享Python調用運籌優(yōu)化求解器（Gurobi、COPT、SCIP）的教程。

2. 求解器介紹

（1）Gurobi
Gurobi是由美國 Gurobi Optimization 公司開發(fā)新一代大規(guī)模優(yōu)化器，提供 C++, Java, Python, .Net, Matlab 和R等多種接口，支持求解以下模型：
（1）連續(xù)和混合整數線性問題
（2）凸目標或約束連續(xù)和混合整數二次問題
（3）非凸目標或約束連續(xù)和混合整數二次問題
（4）含有對數、指數、三角函數、高階多項式目標或約束，以及任何形式的分段約束的非線性問題
（5）含有絕對值、最大值、最小值、邏輯與或非目標或約束的非線性問題

（2）COPT
杉數求解器COPT（Cardinal Optimizer）是杉數自主研發(fā)的針對大規(guī)模優(yōu)化問題的高效數學規(guī)劃求解器套件，也是支撐杉數端到端供應鏈平臺的核心組件，是目前同時具備大規(guī)?；旌险麛狄?guī)劃、線性規(guī)劃（單純形法和內點法）、半定規(guī)劃、（混合整數）二階錐規(guī)劃以及（混合整數）凸二次規(guī)劃和（混合整數）凸二次約束規(guī)劃問題求解能力的綜合性能數學規(guī)劃求解器，為企業(yè)應對高性能求解的需求提供了更多選擇。COPT支持所有主流編程接口：C、C++、C#、Python、Julia、Java、AMPL、GAMS、Pyomo、PuLP、CVXPY。

（3）SCIP
SCIP起源于ZIB(Zuse Institute Berlin)，由Tobias Achterberg奠定整個框架，是目前用于混合整數規(guī)劃（MIP）和混合整數非線性規(guī)劃（MINLP）的最快的非商業(yè)解算器之一，它允許用戶對求解過程進行全面控制，支持自定義搜索樹中的各個模塊，在分支限界(Branch and Bound)過程中添加變量等功能。SCIP支持Python Java AMPL GAMS MATLAB等編程語言。

3. 基礎語言

3.1 創(chuàng)建模型

# 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]) }       # 節(jié)點橫坐標
    y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) }       # 節(jié)點縱坐標
    demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) }         # 節(jié)點需求
    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: 網絡節(jié)點集合
    :param K: 車隊集合
    :param depot: 配送中心id
    :param CAP: 車輛容量
    :param cost: 網絡弧費用集合
    :param demand: 網絡節(jié)點需求集合
    :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("最優(yōu)路徑距離:", cvrp_model.objVal)
    print("最優(yōu)路徑使用車輛數:", 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 = [] #所有節(jié)點
    Q = {} #節(jié)點需求
    TT = {} #節(jié)點旅行時間
    ET = {} #節(jié)點最早開始服務時間
    LT = {} #節(jié)點最晚結束服務時間
    ST = {} #節(jié)點服務時間
    x_coord = {} # 節(jié)點橫坐標
    y_coord = {} # 節(jié)點縱坐標
    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: 所有節(jié)點集合，其中N[0]為車場
    :param Q: 節(jié)點需求集合
    :param TT: 旅行時間
    :param ET: 節(jié)點最早開始服務時間
    :param LT：節(jié)點最晚結束服務時間
    :param ST: 節(jié)點服務時間
    :param CAP: 車輛容量
    :param Cost: 旅行費用
    :param K: 車隊
    :return:
    """
    C = N[1:] #需求節(jié)點
    M=10**5
    depot = N[0]
    # 創(chuàng)建模型
    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)
    # 車輛路徑連續(xù)約束
    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("最優(yōu)路徑距離:", vrptw_model.objVal)
    print("最優(yōu)路徑使用車輛數:", 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]) }       # 節(jié)點橫坐標
    y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) }       # 節(jié)點縱坐標
    demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) }         # 節(jié)點需求
    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):
    # 創(chuàng)建求解環(huán)境
    env = Envr()
    # 創(chuàng)建求解模型
    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("最優(yōu)路徑距離:", cvrp_model.objval)
    print("最優(yōu)路徑使用車輛數:", 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 = [] #所有節(jié)點
    Q = {} #節(jié)點需求
    TT = {} #節(jié)點旅行時間
    ET = {} #節(jié)點最早開始服務時間
    LT = {} #節(jié)點最晚結束服務時間
    ST = {} #節(jié)點服務時間
    x_coord = {} # 節(jié)點橫坐標
    y_coord = {} # 節(jié)點縱坐標
    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: 所有節(jié)點集合，其中N[0]為車場
    :param Q: 節(jié)點需求集合
    :param TT: 旅行時間
    :param ET: 節(jié)點最早開始服務時間
    :param LT：節(jié)點最晚結束服務時間
    :param ST: 節(jié)點服務時間
    :param CAP: 車輛容量
    :param Cost: 旅行費用
    :param K: 車隊
    :return:
    """
    C = N[1:] #需求節(jié)點
    M=10**5
    depot = N[0]
    # 創(chuàng)建求解環(huán)境
    env = Envr()
    # 創(chuàng)建求解模型
    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)
    # 車輛路徑連續(xù)約束
    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("最優(yōu)路徑距離:", vrptw_model.objval)
    print("最優(yōu)路徑使用車輛數:", 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]) }       # 節(jié)點橫坐標
    y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) }       # 節(jié)點縱坐標
    demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) }         # 節(jié)點需求
    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: 網絡節(jié)點集合
    :param K: 車隊集合
    :param depot: 配送中心id
    :param CAP: 車輛容量
    :param cost: 網絡弧費用集合
    :param demand: 網絡節(jié)點需求集合
    :return:
    """
    cvrp_model = Model('cvrp')
    # 添加變量
    X = {}
    Y = {}
    U = {}
    # 創(chuàng)建變量 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})")
    # 創(chuàng)建變量 Y[i,k]
    for i in N:
        for k in K:
            Y[i,k] = cvrp_model.addVar(vtype="B", name=f"y({i},{k})")
    #  創(chuàng)建輔助變量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("最優(yōu)路徑距離:", cvrp_model.getObjVal())
    print("最優(yōu)路徑使用車輛數:", 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 = [] #所有節(jié)點
    Q = {} #節(jié)點需求
    TT = {} #節(jié)點旅行時間
    ET = {} #節(jié)點最早開始服務時間
    LT = {} #節(jié)點最晚結束服務時間
    ST = {} #節(jié)點服務時間
    x_coord = {} # 節(jié)點橫坐標
    y_coord = {} # 節(jié)點縱坐標
    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: 所有節(jié)點集合，其中N[0]為車場
    :param Q: 節(jié)點需求集合
    :param TT: 旅行時間
    :param ET: 節(jié)點最早開始服務時間
    :param LT：節(jié)點最晚結束服務時間
    :param ST: 節(jié)點服務時間
    :param CAP: 車輛容量
    :param Cost: 旅行費用
    :param K: 車隊
    :return:
    """
    C = N[1:] #需求節(jié)點
    M=10**5
    depot = N[0]
    # 創(chuàng)建模型
    vrptw_model=Model()
    X = {}
    T = {}
    # 創(chuàng)建變量 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})")
    # 創(chuàng)建變量 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)
    # 車輛路徑連續(xù)約束
    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("最優(yōu)路徑距離:", vrptw_model.getObjVal())
    print("最優(yōu)路徑使用車輛數:", 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個節(jié)點）：

VRPTW問題數據結構如下（保留100個節(jié)點）：

8. 測試參數

軟件版本號：

Gurbi：9.5.1
SCIP：8.0.0
COPT：5.0.1

CVRP相關參數設置：

車輛容量：80
車隊規(guī)模：10
客戶節(jié)點：30
求解時間：1200s

VRPTW相關參數設置：

車輛容量：80
車隊規(guī)模：30
客戶節(jié)點：100
求解時間：1500s

9. 測試結果

9.1 CVRP求解結果匯總

對于CVRP問題，在相同的求解時間下，Gurobi在求解質量方面有最好的表現(xiàn)，COPT比SCIP略差。

9.2 VRPTW求解結果匯總

對于VRPTW問題，在相同的求解時間下，Gurobi可以找到不錯的優(yōu)化解，而COPT和SCIP并未找到初始解。

9.3 上界下降曲線對比（以CVRP為例）

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

9.5 車輛路徑可視化（以CVRP為例）

c101-31（Gurobi）

c201-31（Gurobi）

r101-31（Gurobi）

10. 小節(jié)

根據以上測試結果，Gurobi的求解性能遠優(yōu)于SCIP和COPT，作為后起之秀COPT已經接近SCIP水平（最新版的COPT尚未測試，也許性能有更大提升）。
雖然求解器可以直接對數學模型進行求解，但僅在小規(guī)模案例上具有較好的性能，對于大規(guī)模問題仍需要設計更高效的算法。
小編后續(xù)將繼續(xù)分享VRP問題系列求解算法，并將其求解效果與上述求解器進行比較。文章來源地址http://www.zghlxwxcb.cn/news/detail-498736.html

到了這里，關于軟件工具 | Python調用運籌優(yōu)化求解器（一）：以CVRP&VRPTW為例的文章就介紹完了。如果您還想了解更多內容，請在右上角搜索TOY模板網以前的文章或繼續(xù)瀏覽下面的相關文章，希望大家以后多多支持TOY模板網！