软件工具 | Python调用运筹优化求解器（二）：以VRPCB为例

import copyimport csvimport mathimport xlsxwriterimport matplotlib.pyplot as pltfrom gurobipy import GRB,Model,quicksum
# 读取文件def read_data(filename):    Q = {}    L = []    B = []    V = []    XY = {}    depot = None    # 读取网络节点、需求    with open(filename, 'r') as f:        node_reader = csv.DictReader(f)        for row in node_reader:            if float(row['linehaul']) == 0 and  float(row['backhaul']) == 0:                depot = row['id']                Q[depot] = 0            elif float(row['linehaul']) > 0 and  float(row['backhaul']) == 0:                L.append(row['id'])                Q[row['id']] = float(row['linehaul'])            elif float(row['linehaul']) == 0 and  float(row['backhaul']) > 0:                B.append(row['id'])                Q[row['id']] = float(row['backhaul'])            V.append(row['id'])            XY[row['id']] = (float(row['x_coord']), float(row['y_coord']))    # 计算网络弧    L0 = L + [depot]    B0 = B + [depot]    AL = [ (i, j) for i in L0 for j in L if i != j ]    AB = [ (i, j) for i in B for j in B0 if i != j]    AC = [ (i, j) for i in L for j in B0 if i != j]    # 计算网络弧距离    Cost = {}    for i in V:        x1, y1 = XY[i][0], XY[i][0]        for j in V:            x2, y2 = XY[j][0], XY[j][0]            Cost[i,j] = math.sqrt( (x1-x2)**2 + (y1-y2)**2 )    return depot,L,L0,B,B0,AL,AB,AC,Q,Cost,XY# 提取结果，形成车辆路径def extract_routes(depot,L,B,B0,X,Y,Z):    L = copy.deepcopy(L)    B = copy.deepcopy(B)    B0 = copy.deepcopy(B0)    route_list = []    V = []    while len(L):        # 提取 派送阶段路径        route = [depot]        cur_node = depot        for j in L:            if X[depot, j].x > 0:                cur_node = j                route.append(j)                L.remove(j)                break        stop = True        while len(L) > 0:            for j in L:                if X[cur_node, j].x > 0:                    cur_node = j                    route.append(j)                    if j != depot:                        L.remove(j)                    stop = False                    break            if stop:                break            else:                stop = True        # 提取 取货阶段路径        for j in B:            if Z[cur_node,j].x > 0:                cur_node = j                route.append(j)                B0.remove(j)                break        if cur_node in B:            while cur_node != depot:                for j in B0:                    if Y[cur_node, j].x > 0:                        cur_node = j                        route.append(j)                        if j != depot:                            B0.remove(j)                        break        else:            route.append(depot)        route_list.append(route)        V.extend(route[1:-1])    print(len(V))    return route_list# 绘制车辆路径def draw_routes(route_list,XY,L,B):    for route in route_list:        path_x = []        path_y = []        for n in route:            path_x.append(XY[n][0])            path_y.append(XY[n][1])        plt.plot(path_x, path_y, linewidth=0.5, ms=5,color='black')        linehual_point_x = [XY[n][0] for n in L]        linehual_point_y = [XY[n][1] for n in L]        backhual_point_x = [XY[n][0] for n in B]        backhual_point_y = [XY[n][1] for n in B]        plt.scatter(linehual_point_x, linehual_point_y, marker='s', c='b', s=5)        plt.scatter(backhual_point_x, backhual_point_y, marker='o', c='r', s=5)        plt.show()# 保存结果def save_file(route_list,total_cost,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,'距离')    for row,route in enumerate(route_list):        route_str = [str(i) for i in route]        dist = sum(Cost[route[i], route[i + 1]] for i in range(len(route) - 1))        ws.write(row + 2, 0, f'{row + 1}')        ws.write(row+2,1,'-'.join(route_str))        ws.write(row + 2, 2, dist)        row += 1    wb.close()# 建模和求解def solve_model(depot,L,L0,B,B0,AL,AB,AC,Q,Cost,K,CAP,XY):    """
    :param depot:车场id    :param L:linehaul节点集合    :param B:backhaul节点集合    :param AL:linehaul节点衔接弧集合    :param AB:backhaul节点衔接弧集合    :param AC:linehaul节点和backhaul节点衔接弧集合    :param Q:节点需求集合    :param Cost:网络弧费用    :return:    """    model = Model()    # 添加变量    X = model.addVars(AL,vtype=GRB.BINARY,name='X[i,j]')    Y = model.addVars(AB,vtype=GRB.BINARY,name='Y[i,j]')    Z = model.addVars(AC,vtype=GRB.BINARY,name='Z[i,j]')    U1 = model.addVars(L0,vtype=GRB.CONTINUOUS,name='U1[i]')    U2 = model.addVars(B0, vtype=GRB.CONTINUOUS, name='U2[i]')    # 目标函数    obj = (quicksum(X[i,j]*Cost[i,j] for i,j in AL) + quicksum(Y[i,j]*Cost[i,j] for i,j in AB) +           quicksum(Z[i,j]*Cost[i,j] for i,j in AC))    model.setObjective(obj,GRB.MINIMIZE)    # linebaul 相关约束    model.addConstr( quicksum(X[depot,j] for j in L) == K ) # 车辆数约束    model.addConstrs( (quicksum(X[i,j] for i in L0 if i != j) == 1 for j in L) ) # 派送需求约束    model.addConstrs( (U1[i] - U1[j] + CAP*X[i,j] <= CAP - Q[j] for i,j in AL) ) # 破圈约束    # backbaul 相关约束    model.addConstr( quicksum([Y[j,depot] for j in B]) + quicksum(Z[i,depot] for i in L) == K ) # 车辆数约束    model.addConstrs( (quicksum(Y[i,j] for j in B0 if i != j) == 1 for i in B) ) # 取货需求约束    model.addConstrs( (U2[i] - U2[j] + CAP*Y[i,j] <= CAP - Q[j] for i,j in AB) ) # 破圈约束    # connection 相关约束    model.addConstr( quicksum(Z[i,j] for i,j in AC) == K ) # 车辆数约束    model.addConstrs( quicksum(X[i,j] for j in L if i != j ) + quicksum(Z[i,j] for j in B0) <= 1 for i in L ) # 接续约束    model.addConstrs( quicksum(Y[i,j] for i in B if i != j) + quicksum(Z[i,j] for i in L) <= 1 for j in B) # 接续约束    # 模型求解    model.Params.TimeLimit = 300  # 规模较大时可设置求解时间限制    model.optimize()    # 判断求解状态    if model.status == GRB.Status.OPTIMAL or model.status == GRB.Status.TIME_LIMIT:        route_list = extract_routes(depot,L,B,B0,X,Y,Z)        draw_routes(route_list, XY, L, B)        save_file(route_list, model.objVal,Cost)if __name__=='__main__':    filename=r'demand-X-n120-50-k3.csv'    depot,L,L0,B,B0,AL,AB,AC,Q,Cost, XY = read_data(filename)    solve_model(depot=depot,L=L,L0=L0,B=B,B0=B0,AL=AL,AB=AB,AC=AC,Q=Q,Cost=Cost,K=3,CAP=21,XY=XY)