Note
Go to the end to download the full example code.
Linear Recourse (Gurobi)¶
Import necessary packages.
import json
import os
import sys
import time
from itertools import product
from benderslib import LShaped, CallbackBase, BendersContext, CST
from benderslib.solvers import Gurobi
from benderslib import LShapedOCGen
from gurobipy import LinExpr
try:
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
except NameError:
sys.path.insert(0, os.path.abspath("."))
from _utils import SMPSReader, first_stage_model, second_stage_model, deterministic_equivalent_model, draw, \
collect_data, limit_memory
Define a callback for the in-out stabilization.
See also
Fischetti, M., Ljubić, I., & Sinnl, M. (2016). Benders decomposition without separability: A computational study for capacitated facility location problems. European Journal of Operational Research, 253(3), 557–569. https://doi.org/10.1016/j.ejor.2016.03.002
Fischetti, M., Ljubić, I., & Sinnl, M. (2017). Redesigning benders decomposition for large-scale facility location. Management Science, 63(7), 2146–2162. https://doi.org/10.1287/mnsc.2016.2461
Warning
This callback is only compatible with single-cut Benders methods.
class InOut(CallbackBase):
def __init__(self, lambda_, alpha, n, m):
self.lambda_ = lambda_
self.alpha = alpha
self.n = n
self.m = m
self.core = None
self.master_linear = None
self.cut_generator = None
self.lb_not_improved_iter_num = 0
def on_sub_build(self, context: BendersContext):
time_start = time.perf_counter()
if self.master_linear is None:
# Initialize the linear relaxation of the master problem.
self.master_linear = context.master_problem.model.relax()
if self.core is None:
# Initialize the core point
self.master_linear.optimize()
self.core = {v: self.master_linear.getVarByName(v).X for v in context.master_problem.complicating_vars}
if self.cut_generator is None:
# Initialize the cut generator.
self.cut_generator = LShapedOCGen(context.master_problem, context.sub_problem, context.benders.params)
cuts = []
constrs = []
current_obj = -float('Inf')
for i in range(self.m):
lambda_ = self.lambda_
if self.lb_not_improved_iter_num >= self.n:
self.lb_not_improved_iter_num = 0
lambda_ = 1
if self.lb_not_improved_iter_num >= self.n * 2:
break
# Update points
point = dict()
for var_name in self.core:
x = self.master_linear.getVarByName(var_name).X
self.core[var_name] = self.alpha * x + (1 - self.alpha) * self.core[var_name]
point[var_name] = lambda_ * x + (1 - lambda_) * self.core[var_name]
# Generate cuts
context.sub_problem.fix_vars(point)
context.sub_problem.prl_solve()
cut = self.cut_generator.generate()[0]
# Add cuts
expr = LinExpr(cut.coefs, [self.master_linear.getVarByName(var_name) for var_name in cut.vars])
if cut.sense == CST.LE:
cons = self.master_linear.addConstr(expr <= cut.rhs)
elif cut.sense == CST.GE:
cons = self.master_linear.addConstr(expr >= cut.rhs)
else:
assert cut.sense == CST.EQ
cons = self.master_linear.addConstr(expr == cut.rhs)
constrs.append(cons)
cuts.append(cut)
# Check lower bound improvement
self.master_linear.optimize()
if self.master_linear.ObjVal > current_obj + 1e-4:
current_obj = self.master_linear.ObjVal
self.lb_not_improved_iter_num = 0
else:
self.lb_not_improved_iter_num += 1
# Detect constraint with positive slack
cut_added_num = 0
for cons, cut in zip(constrs, cuts):
if cons.Slack < 0:
context.master_problem.add_cut(cut)
cut_added_num += 1
# Add cuts to the master problem
end_time = time.perf_counter()
print(f"Generated <{cut_added_num}> cuts in <{(end_time - time_start):.2f}> seconds by InOut callback.")
Solve the instances using different methods and save the results.
@limit_memory(limit_gb=14.5)
def solve(smps_files, instance_name, sample_num, time_limit, solve_methods, seed=1024):
SMPS = SMPSReader(*smps_files, sample_num=sample_num, seed=seed)
SMPS.parse()
ins_file = f"./_ins/{instance_name}_{sample_num}.json"
SMPS.to_json(ins_file)
with open(ins_file, 'r') as f:
data = json.load(f)
# Solve using deterministic equivalent
if "de" in solve_methods:
model = deterministic_equivalent_model(data, enforce_integer=True)
model.setParam('TimeLimit', time_limit)
model.optimize()
model.write(f"./_sol/{instance_name}_de_{sample_num}.json")
# model.write("_temp.mps")
# from pyscipopt import Model
# scip = Model()
# scip.readProblem(filename="_temp.mps")
# scip.optimize()
# scip.writeStatisticsJson(f"./_sol/{instance_name}_de_{sample_num}.json")
# Solve using Benders decomposition
if "bd" in solve_methods:
master_model, complicating_vars = first_stage_model(data, enforce_integer=True)
sub_models, probs = second_stage_model(data)
BD = LShaped.from_models(
master_model=master_model,
master_solver=Gurobi,
sub_model=sub_models,
sub_solver=Gurobi,
complicating_vars=complicating_vars,
prob=probs,
)
BD.register(InOut(lambda_=0.2, alpha=0.3, n=5, m=30))
BD.params.parallel_sub = True
BD.params.use_bnc = True
BD.params.time_limit = time_limit
BD.params.theta_lb = 0
BD.solve()
BD.save(f"./_sol/{instance_name}_bd_{sample_num}.json")
Set 1 Instances
cargo, phone: https://www4.uwsp.edu/math/afelt/slptestset/download.html
lands, storm, gbd: https://pages.cs.wisc.edu/~swright/stochastic/sampling/
We selected these instances based on these criteria:
They are two-stage stochastic programming instances.
They have pure continuous and linear recourse, such that duality based cuts can be used.
Their SMPS
.stofiles are defined in theINDEPmode, allowing resampling scenarios.
Other stochastic programming instance collections:
SPS Resources: https://www.stoprog.org/resources
POSTS: https://users.iems.northwestern.edu/~jrbirge/html/dholmes/post.html
List of Optimization Problem Libraries: https://github.com/ekhoda/optimization_problem_libraries
MSPLib-Library: https://github.com/bonnkleiford/MSPLib-Library
RANDOMRHS 2013: https://users.wpi.edu/~atrapp/randomrhs_2013.htm
def run(solve_methods=None, draw_result=False, dry_run=True):
if solve_methods is None:
solve_methods = ["de", "bd"]
_dir = './set1'
smps_files = {
# Source: https://www4.uwsp.edu/math/afelt/slptestset/download.html
"cargo": (_dir + "/cargo/4node.cor.base", _dir + "/cargo/4node.tim", _dir + "/cargo/4node.sto.32768"),
"phone": (_dir + "/phone/phone.cor", _dir + "/phone/phone.tim", _dir + "/phone/phone.sto"),
# Source: https://pages.cs.wisc.edu/~swright/stochastic/sampling/
"storm": (_dir + "/storm/storm.cor", _dir + "/storm/storm.tim", _dir + "/storm/storm.sto"),
"lands": (_dir + "/lands/lands.cor", _dir + "/lands/lands.tim", _dir + "/lands/lands.sto"),
"gbd": (_dir + "/gbd/gbd.cor", _dir + "/gbd/gbd.tim", _dir + "/gbd/gbd.sto"),
# Note: *cargo* and *storm* were originated from the same problem, but the data is different.
}
sample_nums = [
64,
128,
256,
512,
1024,
]
ins_names = []
de_files = []
bd_files = []
ins_classes = []
s_nums = []
for (ins_class, smps_files), sample_num in product(smps_files.items(), sample_nums):
ins_name = f"{ins_class}"
de_file = f"./_sol/{ins_name}_de_{sample_num}.json"
bd_file = f"./_sol/{ins_name}_bd_{sample_num}.json"
ins_names.append(ins_name)
de_files.append(de_file)
bd_files.append(bd_file)
ins_classes.append(ins_class)
s_nums.append(sample_num)
if not dry_run:
solve(smps_files, ins_name, sample_num, time_limit=3600, solve_methods=solve_methods)
data_points = collect_data(
ins_names=ins_names,
de_files=de_files,
bd_files=bd_files,
ins_classes=ins_classes,
sample_nums=s_nums,
)
if draw_result:
draw(data_points)
return data_points
if __name__ == "__main__":
...
# run(solve_methods=[], draw_result=True)
See also
Tutorial of the Logic-based Benders Decomposition: L-shaped Method
This example uses the following class:
LShaped