三默网为您带来有关“(Python)Gurobi求解CVRP问题”的文章内容,供您阅读参考。

(Python)Gurobi求解CVRP问题

2023-01-16 16:13:27

代码:

import math
import networkx
import gurobipy as gp
from gurobipy import GRB


def vrp(V, c, m, q, Q):
    """
    求解VRP问题
    :param V: The list of nodes
    :param c: c[i,j] represents the distance from node i to node j
    :param m: The number of vehicles
    :param q: q[i] represents the demand of client i
    :param Q: The volume of the vehicle
    :return: model, vrp_callback
    """
    def vrp_callback(model, where):
        # Add constraints to remove the infeasible solution
        if where != GRB.Callback.MIPSOL:
            return
        edges = []
        for (i, j) in x:
            if model.cbGetSolution(x[i, j]) > 0.5:
                if i != V[0] and j != V[0]:
                    edges.append((i, j))
        G = networkx.Graph()
        G.add_edges_from(edges)
        Components = networkx.connected_components(G)
        for S in Components:
            S_card = len(S)
            q_sum = sum(q[i] for i in S)
            NS = int(math.ceil(float(q_sum) / Q))
            S_edges = [(i, j) for i in S for j in S if i < j and (i, j) in edges]
            if S_card >= 3 and (len(S_edges) >= S_card or NS > 1):
                model.cbLazy(gp.quicksum(x[i, j] for i in S for j in S if j > i) <= S_card - NS)
                print("adding cut for", S_edges)
        return

    model = gp.Model("vrp")
    x = {}
    for i in V:
        for j in V:
            if j > i and i == V[0]:  # depot
                x[i, j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)" % (i, j))
            elif j > i:
                x[i, j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)" % (i, j))
    model.update()

    model.addConstr(gp.quicksum(x[V[0], j] for j in V[1:]) == 2 * m, "DegreeDepot")
    for i in V[1:]:
        model.addConstr(gp.quicksum(x[j, i] for j in V if j < i) +
                        gp.quicksum(x[i, j] for j in V if j > i) == 2, "Degree(%s)" % i)

    model.setObjective(gp.quicksum(c[i][j] * x[i, j] for i in V for j in V if j > i), GRB.MINIMIZE)

    model.update()
    model.__data = x
    return model, vrp_callback


def input_data():
    V = range(17)  # The list of nodes
    # 每个节点的需求量,0表示起始节点没有需求
    q = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]  # The list of the demand
    m = 4  # The number of vehicles
    Q = 15  # The volume of the vehicle
    c = [
        [0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354, 468, 776, 662],
        [548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674, 1016, 868, 1210],
        [776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164, 1130, 788, 1552, 754],
        [696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822, 1164, 560, 1358],
        [582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708, 1050, 674, 1244],
        [274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628, 514, 1050, 708],
        [502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856, 514, 1278, 480],
        [194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320, 662, 742, 856],
        [308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662, 320, 1084, 514],
        [194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388, 274, 810, 468],
        [536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764, 730, 388, 1152, 354],
        [502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114, 308, 650, 274, 844],
        [388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194, 536, 388, 730],
        [354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0, 342, 422, 536],
        [468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536, 342, 0, 764, 194],
        [776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274, 388, 422, 764, 0, 798],
        [662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730, 536, 194, 798, 0]
        ]
    return V, c, q, m, Q


def main():
    V, c, q, m, Q = input_data()
    model, vrp_callback = vrp(V, c, m, q, Q)

    model.params.DualReductions = 0
    model.params.LazyConstraints = 1
    model.optimize(vrp_callback)
    x = model.__data

    edges = []
    for (i, j) in x:
        if x[i, j].X > .5:
            if i != V[0] and j != V[0]:
                edges.append((i, j))

    print("*" * 80)
    print("The optimal objective value is:", model.ObjVal)
    print("The optimal path is:")
    print(sorted(edges))
    # [(1, 4), (1, 7), (2, 6), (2, 8), (3, 4), (5, 6), (9, 14), (10, 16), (11, 12), (11, 15), (13, 15), (14, 16)]


if __name__ == '__main__':
    main()

输出:

Using license file C:\Users\86158\gurobi.lic
Changed value of parameter DualReductions to 0
   Prev: 1  Min: 0  Max: 1  Default: 1
Changed value of parameter LazyConstraints to 1
   Prev: 0  Min: 0  Max: 1  Default: 0
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 17 rows, 136 columns and 272 nonzeros
Model fingerprint: 0x9b8d7ee1
Variable types: 0 continuous, 136 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+02, 2e+03]
  Bounds range     [1e+00, 2e+00]
  RHS range        [2e+00, 8e+00]
adding cut for [(2, 12), (2, 15), (3, 10), (3, 15), (4, 7), (4, 12), (5, 10), (5, 13), (7, 9), (8, 9), (8, 14), (13, 14)]
Presolve time: 0.01s
Presolved: 17 rows, 136 columns, 272 nonzeros
Variable types: 0 continuous, 136 integer (120 binary)
adding cut for [(1, 3), (1, 4), (3, 4)]
adding cut for [(2, 6), (5, 8), (5, 9), (6, 8), (9, 14), (14, 16)]
adding cut for [(5, 16), (5, 15), (6, 16), (6, 15)]
adding cut for [(8, 13), (8, 14), (7, 13), (7, 14)]
adding cut for [(9, 11), (9, 12), (10, 11), (10, 12)]

Root relaxation: objective 5.388000e+03, 22 iterations, 0.00 seconds
adding cut for [(1, 4), (2, 6), (2, 10), (3, 4), (3, 15), (6, 8), (10, 16), (11, 12), (11, 15), (12, 13), (13, 14), (14, 16)]
adding cut for [(1, 4), (1, 7), (3, 4), (3, 15), (11, 12), (11, 15)]
adding cut for [(2, 6), (2, 10), (5, 6), (5, 8), (10, 16), (14, 16)]
adding cut for [(1, 4), (1, 5), (2, 6), (2, 10), (3, 4), (3, 15), (5, 6), (9, 14), (10, 16), (11, 12), (11, 15), (14, 16)]
adding cut for [(2, 6), (2, 10), (5, 6), (9, 14), (10, 16), (14, 16)]
adding cut for [(1, 4), (1, 7), (2, 6), (2, 10), (3, 4), (3, 5), (5, 6), (10, 16), (14, 16)]

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 5996.00000    0   16          - 5996.00000      -     -    0s
adding cut for [(1, 4), (3, 4), (3, 15), (11, 12), (11, 15), (12, 13)]
adding cut for [(2, 6), (2, 10), (5, 6), (5, 8)]
     0     0 6016.00000    0    9          - 6016.00000      -     -    0s
adding cut for [(1, 4), (3, 4), (3, 15), (11, 12), (11, 15)]
adding cut for [(2, 6), (2, 10), (9, 14), (10, 16), (14, 16)]
adding cut for [(1, 2), (1, 4), (2, 6), (3, 4), (5, 6), (5, 8)]
     0     0 6196.00000    0   17          - 6196.00000      -     -    0s
adding cut for [(9, 10), (10, 16), (11, 12), (11, 15), (13, 14), (13, 15), (14, 16)]
adding cut for [(10, 16), (11, 12), (11, 15), (13, 14), (13, 15), (14, 16)]
     0     0 6196.00000    0   17          - 6196.00000      -     -    0s
adding cut for [(5, 6), (5, 8), (6, 10), (10, 16), (14, 16)]
adding cut for [(9, 12), (11, 12), (11, 15), (13, 15)]
H    0     0                    6436.0000000 6196.00000  3.73%     -    0s
*    0     0               0    6208.0000000 6208.00000  0.00%     -    0s

Cutting planes:
  Lazy constraints: 20

Explored 1 nodes (62 simplex iterations) in 0.08 seconds
Thread count was 8 (of 8 available processors)

Solution count 2: 6208 6436 

Optimal solution found (tolerance 1.00e-04)
Best objective 6.208000000000e+03, best bound 6.208000000000e+03, gap 0.0000%

User-callback calls 79, time in user-callback 0.01 sec
********************************************************************************
The optimal objective value is: 6208.0
The optimal path is:
[(1, 4), (1, 7), (2, 6), (2, 8), (3, 4), (5, 6), (9, 14), (10, 16), (11, 12), (11, 15), (13, 15), (14, 16)]

Process finished with exit code 0

附录:技术细节

参数名称:LazyConstraints

通过回调添加惰性约束的程序必须将此参数设置为值1。

该参数告诉Gurobi避免某些与惰性约束不兼容的缩减和转换。

注意:如果通过设置 Lazy 属性(而不是通过回调)使用 lazy 约束,则不需要设置此参数。

LazyConstraints

Programs that use lazy constraints must set this parameter

Type: int
Default value: 0
Minimum value: 0
Maximum value: 1

Programs that add lazy constraints through a callback must set this parameter to value 1.

The parameter tells the Gurobi algorithms to avoid certain reductions and transformations that are incompatible with lazy constraints.

Note that if you use lazy constraints by setting the Lazy attribute (and not through a callback), there’s no need to set this parameter.

Note: Only affects mixed integer programming (MIP) models.

对应语句:

model.params.LazyConstraints = 1

参数名称:DualReductions

确定是否在预处理中执行双重缩减

Controls dual reductions

Type: int
Default value: 1
Minimum value: 0
Maximum value: 1

Determines whether dual reductions are performed in presolve.

You should disable these reductions if you received an optimization status of INF_OR_UNBD and would like a more definitive conclusion.

对应语句:

model.params.DualReductions = 0

Model.optimize()

optimize  ( callback=None )

Optimize the model.

The algorithm used for the optimization depends on the model type (simplex or barrier for a continuous model; branch-and-cut for a MIP model).

Upon successful completion, this method will populate the solution related attributes of this model.

Arguments:

callback: Callback function.

The callback function should take two arguments, model and where.

During the optimization, the (callback) function will be called periodically, with model set to the model being optimized, and where indicating where in the optimization the callback is called from.

对应语句:

model.optimize(vrp_callback)

Callbacks

Gurobi callback class.

A callback is a user function that is called periodically by the Gurobi optimizer in order to allow the user to query or modify the state of the optimization.
(译:callback是一个用户函数,这个函数被Gurobi优化器定期调用,为了允许用户查询或修改优化的状态)

More precisely, if you pass a function that takes two arguments (model and where) as the argument to Model.optimize, your function will be called during the optimization.
(译:更准确地说,如果向Model.optimize传入一个以modelwhere作为实参的函数,那么您的函数将在优化的过程中被调用)

Your callback function can then call Model.cbGet to query the optimizer for details on the state of the optimization.
(译:您的callback函数可以进一步调用Model.cbGet来询问优化器 优化状态的细节)

Gurobi callbacks can be used both to monitor the progress of the optimization and to modify
the behavior of the Gurobi optimizer.
(译:Gurobi的callbacks可以用来监测优化过程和修改Gurobi优化器的行为)

A simple user callback function might call Model.cbGet to produce a custom display, or perhaps to terminate optimization early (using Model.terminate).
(译:一个简单的用户callback函数可能会调用Model.cbGet来产生自定义显示,或可能使用Model.terminate来提前终结优化)

More sophisticated MIP callbacks might use Model.cbGetNodeRel or Model.cbGetSolution to retrieve values from the solution to the current node, and then use Model.cbCut or Model.cbLazy to add a constraint to cut off that solution, or Model.cbSetSolution to import a heuristic solution built from that solution.
(译:更复杂的MIP callbacks也许会使用Model.cbGetNodeRelModel.cbGetSolution来取回当前节点解的数值,接下来使用Model.cbCutModel.cbLazy来加入约束来切掉那个解,或使用Model.cbSetSolution来导入从那个解构造的启发式解)

For multi-objective problems, you might use Model.cbStopOneMultiObj to interrupt the optimization process of one of the optimization steps in a multi-objective MIP problem without stopping the hierarchical optimization process.
(译:我暂时不研究多目标优化,所以不翻译此句了。)

The Gurobi callback class provides a set of constants that are used within the user callback function.
(译:Gurobi的callback类提供了一系列在用户callback函数中使用的常数)

The first set of constants in this class list the options for the where argument to the user callback function.
(译:callback类的第一部分常数 罗列了 用户callback函数中实参where的选择)

The where argument indicates from where in the optimization process the user callback is being called.
(译:实参where表明,在优化过程中,用户callback函数从哪里被调用)

Options are listed in the Callback Codes section of this document.

The other set of constants in this class list the options for the what argument to Model.cbGet.
(译:另一部分参数罗列了Model.cbGet的实参what的选项)

The what argument is used by the user callback to indicate what piece of status information it
would like to retrieve.
(译:用户callback函数使用实参what来表明用户想重新获取哪些状态信息)

As with the where argument, you refer to a what constant through GRB.Callback.
(译:就像实参where一样,您通过GRB.Callbcak来查询what

If you would like to pass data to your callback function, you can do so through the Model object.
(译:如果您想向callback函数传参,那么你可以通过Model对象来实现)

For example, if your program includes the statement model._value = 1 before the optimization begins, then your callback function can query the value of model._value.
(译:例如,如果您的程序在优化开始前包含语句model._value = 1,那么您的callback函数可以查询model._value的取值)

Note that the name of the user data field must begin with an underscore.
(译:注意用户数据域必须由下划线开始)!!!
这一点特别重要,这也是为什么gurobi的callback代码总是出现一些下划线!!!

When solving a model using multiple threads, the user callback is only ever called from a single thread, so you don’t need to worry about the thread-safety of your callback.
(多线程相关,暂不翻译)

Note that changing parameters from within a callback is not supported, doing so may lead to undefined behavior.

Model.addVar()

addVar  (  lb=0.0, ub=float('inf'), obj=0.0, vtype=GRB.CONTINUOUS, name="", column=None  )

Add a decision variable to a model. 为一个模型增加一个决策变量.

Arguments
lb (optional): Lower bound for new variable.
ub (optional): Upper bound for new variable.
obj (optional): Objective coefficient for new variable.
vtype (optional): Variable type for new variable.
name (optional): Name for new variable. 注意不要带有"➡"或空格.
column (optional): Column object that indicates the set of constraints in which the new
variable participates, and the associated coefficients.

Return value
New variable object.

Example usage:

import gurobipy as gp
from gurobipy import GRB

mdl = gp.Model()

x = mdl.addVar()  # all default arguments
y = mdl.addVar(vtype=GRB.INTEGER, obj=1.0, name="y")  # arguments by name
z = mdl.addVar(0.0, 1.0, 1.0, GRB.BINARY, "z")  # arguments by position

print("x", x)
print("y", y)
print("z", z)

vtype到底有哪些?

(GRB.CONTINUOUS, GRB.BINARY, GRB.INTEGER, GEB.SEMICONT, or GRB.SEMIINT)