Solving the CVRP variant problem

ユーザーの入力を待っています。

2023年11月11日 14:10

I want to implement a variant of the CVRP problem using gurobi. Specifically, he has n customer nodes (1..n) and a warehouse node depot=0. The capacity of the warehouse is set to Q=1. And assume that the distance between two points is known, and the distance between any two points has two values d^{+} and d^{-}. The mathematical model is modeled as:

the data is:

n:5
distance:(i,j,d_ij^{+},d_ij^{-})
0 1 0.74 0.37
0 2 0.68 0.44
0 3 0.26 0.2
0 4 0.7 0.49
0 5 0.64 0.28
1 2 0.87 0.36
1 3 0.68 0.17
1 4 0.86 0.69
1 5 0.68 0.48
2 3 0.42 0.27
2 4 0.88 0.36
2 5 0.51 0.25
3 4 0.66 0.63
3 5 0.66 0.48
4 5 0.57 0.21
node_demand:1..n) depot_demand=0
0.4
0.2
0.1
0.6
0.1

the code is as follows, But the current situation is that even for 5 points, the optimal solution cannot be found.

import gurobipy as gp
from gurobipy import GRB
import time
EPS = 0.0001


def get_wst_scenario(n, sol, d_down, d_up):
    wst = dict()
    for i in range(n+1):
        for j in range(i+1,n+1):
            if sol[(i,j)]>0.5 or sol[(j,i)]>0.5:
                wst[(i, j)] = wst[(j, i)] = d_up[(i,j)]
            else:
                wst[(i, j)] = wst[(j, i)] = d_down[(i,j)]

    return wst


def get_regret(N, Q, q, d_down, d_up, sol):
    d_wst = get_wst_scenario(n=len(N), sol=sol, d_down=d_down, d_up=d_up)
    y_val, y_sol = solve_cvrp_bigM(N=N, E=d_down.keys(), d=d_wst, Q=Q, q=q)
    cost_x = sum(d_up[e] * sol[e] for e in sol)
    regret = sum(d_up[e] * sol[e] for e in sol) - y_val

    return regret, cost_x, y_val, y_sol


#---V：include depot node 0  N：all customer node(i..n)  E：all edges
#   d: cost  q:demand  Q:capacity limitation for the vehicle---#
def solve_cvrp_bigM(N, E, d, Q, q, time_limit=None):
    V = [0] + N
    model = gp.Model("CVRP")
    ### Variable
    # x 0-1
    x = {e: model.addVar(vtype=GRB.BINARY, name="x[{}]".format(e)) for e in E}

    u = model.addVars(V,vtype=GRB.CONTINUOUS)
    u[0] = 0  # 在仓库节点的已服务量为 0
    model.update()

    ### Objective
    model.setObjective(gp.quicksum(d[e] * x[e] for e in E), GRB.MINIMIZE)

    ### Constraints
    #
    model.addConstrs(gp.quicksum(x[i,j] for j in V if i != j) == 1 for i in N)
    model.addConstrs(gp.quicksum(x[i,j] for i in V if i != j) == 1 for j in N)
    #
    model.addConstrs( ( (u[i] + q[j]) <= (u[j] + Q*(1-x[i,j])) ) for i in V for j in N if i!=j)
    model.addConstrs( ( (u[i] + q[j])
                      >= ( u[j] - (Q - q[i] - q[j])*(1-x[i,j]) ) ) for i in V for j in N if i!=j)


    model.Params.outputFlag = False
    model.Params.threads = 1
    model.Params.MIPGap = 0.0
    model.Params.lazyConstraints = 1
    if time_limit is not None:
        model.Params.timeLimit = time_limit
    model.optimize()
    if model.SolCount <= 0:
        return None, None
    x_sol = {e: round(x[e].x) for e in x}
    return (model.ObjVal), x_sol



def set_bd_model(N, E, Q, q, d_down, d_up):
    V = [0] + N
    model = gp.Model("BD")
    # x:
    x = {e: model.addVar(vtype=GRB.BINARY, name="x[{}]".format(e)) for e in E}
    # r:  r<=sum(yl+sum y(u-l)x
    r = model.addVar(vtype=GRB.CONTINUOUS, ub=sum(d_up.values() )/2 , name="r")
    u = model.addVars(V, vtype=GRB.CONTINUOUS, name="u")
    u[0] = 0  #
    model.update()

    model.setObjective(gp.quicksum(d_up[e] * x[e] for e in x) - r, GRB.MINIMIZE)
    # flow cut
    #
    model.addConstrs(gp.quicksum(x[i, j] for j in V if j != i) == 1 for i in N)
    model.addConstrs(gp.quicksum(x[i, j] for i in V if i != j) == 1 for j in N)

    #
    model.addConstrs(((u[i] + q[j]) <= (u[j] + Q * (1 - x[i, j]))) for i in V for j in N if i != j)
    model.addConstrs(((u[i] + q[j])
                      >= (u[j] - (Q - q[i] - q[j]) * (1 - x[i, j]))) for i in V for j in N if i != j)

    model.Params.outputFlag = False
    model.Params.threads = 1
    model.Params.MIPGap = 0.0

    model._x, model._r = x, r
    model._N = N
    model._Q = Q
    model._q = q
    model._d_down, model._d_up = d_down, d_up
    #
    model.Params.lazyConstraints = 1
    model.update()

    return model, x, r


def gen_cut(mod, where): # where 哪个阶段触发callback
    """Callback to add a cut for branch-and-cut framework for benders mod"""
    # Execute the function when an incumbent is found
    if where != GRB.Callback.MIPSOL:
        return

    x_sol = mod.cbGetSolution(mod._x)

    # check Benders cut
    ttb = mod.cbGet(GRB.Callback.RUNTIME)

    # Obtain the incumbent solution
    r_sol = mod.cbGetSolution(mod._r)
    d_up, d_down = mod._d_up, mod._d_down

    # Prepare worst-case scenario
    ####
    d_wst = get_wst_scenario(n=len(mod._N), sol=x_sol, d_down=d_down, d_up=d_up)

    y_val, y_sol = solve_cvrp_bigM(N=mod._N, E=d_down.keys(), d=d_wst, Q=mod._Q, q=mod._q)
    ####

    y_sol = {e for e in d_down if y_sol[e] > 0.5} # set对象


    for e in x_sol:
        if x_sol[e]>0.5:
            i = e[0]
            j = e[1]
            x_sol[i,j] = x_sol[j,i]

    if y_val + EPS < r_sol:
        # Add bd cuts

        mod.cbLazy(gp.quicksum(d_down[e] + (d_up[e] - d_down[e]) * x_sol[e] for e in y_sol) >= mod._r)
        return
    mod._ttb = ttb
    return


def solve_bc(n, Q, q, d_down, d_up, time_limit):
    N = [i for i in range(1, n + 1)]
    model, x, _ = set_bd_model(N, d_down.keys(), Q, q, d_down, d_up)  # return model, x, r

    model.Params.timeLimit = time_limit
    model.optimize(gen_cut)
    if model.SolCount <= 0:
        return None, None, None

    sol = [e for e in x if x[e].x > 0.5]

    x_e = dict()
    for i in range(n+1):
        for j in range(n+1):
            if i!=j:
                x_e[(i,j)] = 0

    for i in range(n+1):
        for j in range(n+1):
            if (i,j) in sol:
                x_e[(i,j)] = 1
                # x_e[(j, i)] = 1

    obj = (model.objVal)
    bound = (model.objBound + 1 - EPS)
    return (obj, bound, sol, x_e)

Solving the CVRP variant problem

コメント

お探しのものを見つけられませんでしたか？