Pareto-optimal Cut

Prepare the problem for Benders decomposition.

from benderslib import ClassicalBenders, AnnotatedBenders, CallbackBase, ClassicalOC
from benderslib.solvers import Gurobi
from benderslib.utils import draw_curve

from gurobipy import Model, GRB


def make_original_problem():
    model = Model("Original")

    n_vars = 20
    y = model.addVars(n_vars, name="y", lb=1, ub=40, vtype=GRB.INTEGER)
    z = model.addVars(n_vars, name="z", lb=1, ub=40, vtype=GRB.CONTINUOUS)

    model.addConstr(y.sum() + z.sum() <= 50 * n_vars, "main_constr_yz")
    model.addConstrs((2 * y[i] <= 2 * (i + 1) for i in range(n_vars)), name="constr_y")
    model.addConstrs((2 * y[i] + z[i] >= i for i in range(n_vars)), name="constr_yz")
    model.addConstrs((3 * z[i] <= 15 for i in range(n_vars)), name="constr_z")

    model.setObjective(2 * y.sum() + 3 * z.sum(), sense=GRB.MINIMIZE)

    model.Params.OutputFlag = 0
    model.Params.LogToConsole = 0

    model.update()
    complicating_vars = [v.VarName for v in y.values()]
    return model, complicating_vars

Define the (approximate) Pareto-optimal cut callback.

class ParetoCut(CallbackBase):

    def __init__(self):
        # The core point will be iteratively updated.
        self.core_point = dict()

    def on_iteration_start(self, context):
        if context.state.n_iter <= 1:
            # No master solution is available at the beginning of the first iteration.
            return

        for var in context.current_comp_vals:

            # The approximate core point is updated as the average of
            # the current core point and the current master solution.
            # For the first iteration, the core point is simply set
            # to be the current master solution.

            if var in self.core_point:
                self.core_point[var] = (self.core_point[var] + context.current_comp_vals[var]) / 2
            else:
                self.core_point[var] = context.current_comp_vals[var]

        # Solve the Magnanti–Wong problem with the approximate core point.
        context.sub_problem.fix_vars(self.core_point)
        context.sub_problem.solve()

        # Generate an (approximate) Pareto-optimal cut.
        var_coefs = context.sub_problem.get_var_coefs(context.master_problem.complicating_vars)
        rhs = context.sub_problem.get_rhs()
        dual_values = context.sub_problem.get_dual_values()
        cut = ClassicalOC(context.master_problem.complicating_vars, var_coefs, dual_values, rhs)

        context.master_problem.add_cut(cut)

    def on_opti_cut_generated(self, context):
        if context.state.n_iter <= 1:
            # Use the master problem solution, instead of the core point,
            # to generate cuts in the first iteration, since the core point is not updated yet.
            return

        # Cut generation is taken over by the `on_iteration_start` callback,
        # so we clear the optimality cuts generated in the current iteration.
        context.current_opti_cuts = []

This is an implementation of the Algorithm 3 from the following paper. In the original Magnanti–Wong method, to obtain the Pareto-optimal cut, an additional optimization problem is required to be solved with the objective value of the original subproblem and a core point, which is usually difficult to be determined in practice. Papadakos provide a practical enhancement to the Magnanti–Wong method by using approximate core points, and the method does not require solution from the original subproblem, which is more efficient and easier to be implemented. However, it is worth noting that the cut generated can be not Pareto-optimal, relying on the choice of the core points.

See also

Papadakos, N. (2008). Practical enhancements to the Magnanti–Wong method. Operations Research Letters, 36(4), 444–449. https://doi.org/10.1016/j.orl.2008.01.005

Run with the callback.

model, complicating_vars = make_original_problem()
model_copy = model.copy()

BD = AnnotatedBenders(
    model,
    solver=Gurobi,
    complicating_vars=complicating_vars,
    benders=ClassicalBenders
)

callback = ParetoCut()
BD.register(callback)

BD.solve()
draw_curve(BD.result)

Hint

With stronger cuts, in the early stage of the algorithm, both the upper and lower bounds improve significantly.

Run without the callback.

BD = AnnotatedBenders(
    model_copy,
    solver=Gurobi,
    complicating_vars=complicating_vars,
    benders=ClassicalBenders
)

# BD.solve()
# draw_curve(BD.result)

See also

Tags: benders: classical, solver: gurobi, deterministic, callback, enhancement

Gallery generated by Sphinx-Gallery