【CVRP&VRPTW】Python+Gurobi求解运输问题建模实践二

采用Python+Gurobi实现车辆路径规划问题建模求解【CVRP&VRPTW】

今天以经典CVRP和VRPTW为例，讲解如何采用python+gurobi编程求解。个人认为采用求解器求解数学模型的关键在于拿到一个正确的数学模型！本人在求解CVRP时，通过文献库搜索了很多文章，限于个人能力，无法采用gurobi对其模型进行求解。最后借用了微信博文中提到的数学模型。（文末给出链接）

往期资料

• python实现6种智能算法求解CVRP问题
• python实现7种智能算法求解MDVRP问题
• python实现7种智能算法求解MDVRPTW问题
• Python版MDHFVRPTW问题智能求解算法代码【TS算法】
• Python版MDHFVRPTW问题智能求解算法代码【SA算法】
• Python版MDHFVRPTW问题智能求解算法代码【GA算法】
• Python版MDHFVRPTW问题智能求解算法代码【DPSO算法】
• Python版MDHFVRPTW问题智能求解算法代码【DE算法】
• Python版MDHFVRPTW问题智能求解算法代码【ACO算法】

1. CVRP

1.1 适用场景

• 求解CVRP
• python+gurobi

1.2 数学模型

后续代码所对应的数学模型源自以下文章：

https://mp.weixin.qq.com/s/DYh-5WkrYxk1gCKo8ZjvAw

1.3 数据结构

采用.xlsx文件储存数据文件，共包含两个sheet，命名为：node，link。
在 ‘node’ sheet中内容如下：

在 ‘link’ sheet中内容如下：

1.4 编程实现

（1）读取文件

def readXlsxFile(filename):
    df=pd.read_excel(filename,sheet_name=None)
    "读取节点信息"
    df_node=df['node']
    depot=df_node['id'][0]
    N = [depot]
    C = []
    Q = {}
    for i in range(df_node.shape[0]):
        id=df_node['id'][i]
        demand=df_node['demand'][i]
        N.append(id)
        C.append(id)
        Q[id]=int(demand)
    N = tuplelist(N)
    C = tuplelist(C)
    Q=tupledict(Q)

    "读取网络弧信息"
    Cost={}
    df_link=df['link']
    for i in range(df_link.shape[0]):
        from_node_id=df_link['from_node_id'][i]
        to_node_id=df_link['to_node_id'][i]
        cost=df_link['link_cost'][i]
        Cost[from_node_id,to_node_id]=cost
    Cost=tupledict(Cost)
    return depot,C,N,Q,Cost

（2）模型构建与求解

def solveCVRPModel(depot,C,N,Q,Cost,n_vehicles=20,CAP=100):
    """
    :param depot: 车场索引
    :param C: 需求点集合
    :param N: 所有点集合
    :param Q: 需求集合
    :param Cost: 弧运输成本集合
    :param n_vehicles: 最大车辆数量
    :param CAP: 车辆容量
    :return:
    """
    "构建车队"
    K=tuplelist([f'v'+str(k) for k in range(n_vehicles)])
    "实例化模型"
    model=Model('cvrp')
    "添加变量"
    X=model.addVars( N,N,K,vtype=GRB.BINARY,name='X[i,j,k]')
    Y=model.addVars( K,N,vtype=GRB.BINARY,name='Y[k,i]')
    U=model.addVars( K,N,vtype=GRB.CONTINUOUS,name='U[k,i]')
    "目标函数：最小化路径成本"
    z1=quicksum( Cost[i,j]*X[i,j,k] for i in N for j in N for k in K if i!=j)
    model.setObjective(z1,GRB.MINIMIZE)
    "约束：需求覆盖约束"
    model.addConstrs( quicksum( Y[k,i] for k in K ) ==1 for i in C )
    "约束：车辆数量约束"
    model.addConstr( quicksum( Y[k,depot] for k in K) == n_vehicles )
    "约束：流平衡"
    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 )
    "约束：决策变量关联"
    model.addConstrs( quicksum( X[i,j,k] for j in N ) == Y[k,i] for i in N for k in K )
    "约束：容量限制"
    model.addConstrs( quicksum( Q[i]*Y[k,i] for i in C ) <= CAP for k in K )
    "约束：破除子环"
    model.addConstrs( U[k,i]-U[k,j]+CAP*X[i,j,k] <= CAP-Q[i] for i in C for j in C for k in K )
    model.addConstrs( Q[i] <= U[k,i] for k in K for i in C )
    model.addConstrs( U[k,i] <= CAP for k in K for i in C)
    "求解"
    model.Params.TimeLimit = 300 # 设置求解时间上限
    model.optimize()
    if model.status == GRB.Status.OPTIMAL or model.status == GRB.Status.TIME_LIMIT:
        for k in K:
            for i in N:
                for j in N:
                    if i!=j :
                        if X[i,j,k].x>0:
                            print("X[{},{},{}]=1".format(i,j,k))
        print("obj:{}".format(model.objVal))
    else:
        print("no solution")

2. VRPTW

2.1 适用场景

• 求解VRPTW
• python+gurobi

2.2 数学模型

后续代码所对应的数学模型源自以下文章：

https://mp.weixin.qq.com/s/tF-ayzjpZfuZvelvItuecw

2.3 数据结构

采用.csv文件格式储存数据文件，分别为：node.csv 和link.csv
在 ‘node.csv’ 中内容如下：

在 ‘link.csv’ 中内容如下：

2.4 编程实现

（1）读取文件

def readCsvFile(node_file,link_file):
    """
    :param node_file: 节点文件，
    :param link_file: 网络弧文件
    :return:
    """
    N=[] #所有节点
    Q={} #节点需求
    TT={} #节点旅行时间
    ET={} #节点最早开始服务时间
    LT={} #节点最晚结束服务时间
    ST={} #节点服务时间
    Cost={}
    with open(node_file,'r') as f:
        node_reader=csv.DictReader(f)
        for row in node_reader:
            node_id = row['id']
            demand = float(row['demand'])
            start_time = float(row['start_time'])
            end_time = float(row['end_time'])
            service_time = float(row['service_time'])
            N.append(node_id)
            Q[node_id] = demand
            ET[node_id] = start_time
            LT[node_id] = end_time
            ST[node_id] = service_time
    with open(link_file,'r') as f:
        link_reader = csv.DictReader(f)
        for row in link_reader:
            from_node_id = row['from_node_id']
            to_node_id = row['to_node_id']
            travel_time = float(row['travel_time'])
            travel_cost = float(row['link_cost'])
            TT[from_node_id,to_node_id] = travel_time
            Cost[from_node_id,to_node_id] = travel_cost
    return N,Q,TT,ET,LT,ST,Cost

（2）模型构建与求解

def solveVRPTWModel(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=tuplelist(N[1:]) #需求节点
    N=tuplelist(N)
    Q=tupledict(Q)
    TT=tupledict(TT)
    E=tupledict(ET)
    L=tupledict(LT)
    S=tupledict(ST)
    K=tuplelist([f'v'+str(i) for i in range(K)])
    M=10**5
    depot = N[0]
    "创建模型"
    model=Model()
    "添加变量"
    X=model.addVars(N,N,K,vtype=GRB.BINARY,name='X(i,j,k)')
    T=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)
    model.setObjective(z1,GRB.MINIMIZE)
    "车辆起点约束"
    model.addConstrs(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
    "车辆路径连续约束"
    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)
    "车辆终点约束"
    model.addConstrs(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
    "需求服务约束"
    model.addConstrs( quicksum(X[i,j,k] for k in K for j in N if j!=i)==1 for i in C)
    "车辆容量约束"
    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 )
    "时间窗约束"
    model.addConstrs( T[i,k]+S[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 )
    model.addConstrs( T[i,k] >= E[i] for i in N for k in K)
    model.addConstrs( T[i,k] <= L[i] for i in N for k in K)
    "设置模型参数"
    model.Params.TimeLimit = 300 # 规模较大时可设置求解时间限制
    # model.setParam('OutputFlag', 0) #是否输出求解日志
    "模型求解"
    model.optimize()
    "判断求解状态"
    if model.status == GRB.Status.OPTIMAL or model.status == GRB.Status.TIME_LIMIT:
        print('obj={}'.format(model.objVal))
        res=[]
        for k in K:
            for i in N:
                for j in N:
                    if i!=j:
                        if X[i,j,k].x>0:
                            print("X[{},{},{}]=1".format(i,j,k))
                            res.append([i,j,k,T[i,k].x,T[j,k].x])
        saveFile(res)
    else:
        print("no solution")

（3）保存结果

def saveFile(data):
    outfile = open('results.csv', 'w', newline='', errors='ignore')
    write = csv.writer(outfile)
    write.writerow(['from_node_id', 'to_node_id', 'vehicle','Ti','Tj'])
    for v in data:
        write.writerow(v)
    outfile.close()

3. 完整代码

如有错误，欢迎交流。
代码和数据文件可从github主页免费获取：

https://github.com/PariseC/modeling_examples_using_gurobi_in_python

可关注微信公众号，及时获取新内容（推送内容正在准备中）：

Python助力交通

参考

1. https://mp.weixin.qq.com/s/DYh-5WkrYxk1gCKo8ZjvAw
2. https://mp.weixin.qq.com/s/tF-ayzjpZfuZvelvItuecw
3. http://www.gurobi.cn/pic.asp?bigclassname=%D1%A7%CF%B0%D7%CA%C1%CF