Yes, Benders decomposition applies even when the LP subproblem has only one block. In fact, the original paper by Benders in 1962 treated only that case. But exploiting block structure can certainly boost performance.
For completeness, here is one implementation of Benders decomposition for the capacitated version of your problem. To avoid the complications of dynamically generated feasibility cuts, I made VALID_MAPS the fully dense cross product of ORDERS and LOCATIONS. That change and the introduction of the SufficientCapacity constraint together imply that the subproblem will always be feasible. With max_orders_per_location = 5 and an optimality gap of 5%, scenario 2 solves quickly on my laptop.
%let relobjgap = 0.05;
proc optmodel printlevel=0;
/* Benders decomposition parameters */
num infinity = constant('BIG');
num best_lower_bound init 0; /* some finite lower bound on the original objective */
num best_upper_bound init infinity;
num gap = (if best_lower_bound ne 0 and best_upper_bound < infinity
then (best_upper_bound - best_lower_bound) / abs(best_lower_bound)
else .);
/* declare MILP master problem over complicating variables */
*Define data elements;
put 'optmodel: Define data elements';
set <num> ORDERS;
set <num> LOCATIONS;
num distance {ORDERS, LOCATIONS} init &inf_dist.;
*Define variable 1;
var On {LOCATIONS} binary init 0;
*Read datasets 1;
put 'optmodel: Read datasets 1';
read data work.orders into ORDERS=[ORDER_ID];
read data work.points into LOCATIONS=[LOCATION_ID] On=on;
read data work.distances into [ORDER_ID LOCATION_ID] distance=&distance_var.;
*Define variable 2;
/* set VALID_MAPS = {p in ORDERS, w in LOCATIONS: distance[p,w] <= (&inf_dist.-10)};*/
set VALID_MAPS = ORDERS cross LOCATIONS; /* make dense to avoid feasibility cuts */
/* var Map {VALID_MAPS} binary init 0;*/
var Map {VALID_MAPS} >= 0; /* relax to get correct duals */
*Read datasets 2;
put 'optmodel: Read datasets 2';
read data work.distances into [ORDER_ID LOCATION_ID] Map=map;
*Calculations;
put 'optmodel: Calculations';
impvar total_LOCATIONs_on = sum{w in LOCATIONS} On[w];
impvar total_maps_on = sum{<p,w> in VALID_MAPS} Map[p,w];
impvar total_distance = sum{<p,w> in VALID_MAPS} Map[p,w]*distance[p,w];
impvar total_penalty = total_LOCATIONs_on * &LOCATION_distance_penalty.;
*Mechanical constraints;
put 'optmodel: Mechanical constraints';
con ORDERs_mapped_correctly {p in ORDERS}:
/* (sum{<(p),w> in VALID_MAPS} Map[p,w]) = 1;*/
sum{w in LOCATIONS} Map[p,w] = 1;
set SUPPORT init {};
* con LOCATIONs_mapped_correctly {w in LOCATIONS}:
(sum{<p,(w)> in VALID_MAPS} Map[p,w]) <= (w in SUPPORT) * &max_ORDERs_per_LOCATION.;
con LOCATIONs_mapped_correctly {<p,w> in VALID_MAPS}:
-Map[p,w] >= -(w in SUPPORT);
num capacity {w in LOCATIONS} = min(&max_ORDERs_per_LOCATION., card(ORDERS));
con CapacityCon {w in LOCATIONS}:
-sum{p in ORDERS} Map[p,w] >= -(w in SUPPORT) * capacity[w];
con SufficientCapacity:
sum {w in LOCATIONS} capacity[w] * On[w] >= card(ORDERS);
*Business constraints;
put 'optmodel: Business constraints';
con max_LOCATIONs_on:
total_LOCATIONs_on <= &final_max_LOCATIONs_on.;
*Minimize distance;
put 'optmodel: Set Goals';
min goal = total_distance + total_penalty;
*Initial solution check (seed);
print "VALUES BEFORE SOLVING" total_LOCATIONs_on total_maps_on goal total_distance total_penalty;
/* master objective */
var Eta >= 0;
num Eta_sol;
min MasterObjective = total_penalty + Eta;
/* optimality cuts */
num num_optcuts init 0;
set OPTCUTS = 1..num_optcuts;
num cut_constant {OPTCUTS};
num cut_coefficient {LOCATIONS, OPTCUTS} init 0;
con OptimalityCut {optcut in OPTCUTS}:
Eta >= cut_constant[optcut] + sum {w in LOCATIONS} cut_coefficient[w,optcut] * On[w];
/* omit feasibility cuts (assume SufficientCapacity implies feasibility of subproblem) */
/* additional valid cuts */
/* smallest positive distance from p */
num spd {p in ORDERS} = min {w in LOCATIONS diff {p}} distance[p,w];
/* if On[p] = 0 then Eta[p] >= spd[p] */
* con OnOrPositiveDistance {p in ORDERS}:
Eta[p] >= spd[p] * (1 - On[p]);
con OnOrPositiveDistance:
Eta >= sum {p in ORDERS} spd[p] * (1 - On[p]);
problem Master include
On Eta
MasterObjective
SufficientCapacity
max_LOCATIONs_on
OptimalityCut OnOrPositiveDistance;
/* declare subproblem over Map here (same as original but with On fixed) */
min SubproblemObjective = sum{<p,w> in VALID_MAPS} Map[p,w]*distance[p,w];
problem Subproblem include
Map
SubproblemObjective
ORDERs_mapped_correctly LOCATIONs_mapped_correctly CapacityCon;
num lower_bound;
num upper_bound = SubproblemObjective.sol + &LOCATION_distance_penalty. * card(SUPPORT);
set INCUMBENT init {};
num Map_sol {VALID_MAPS};
/* main loop */
do until (gap ne . and gap <= &relobjgap);
put 'Solving Master...';
use problem Master;
put num_optcuts=;
for {w in LOCATIONS} On[w] = (w in INCUMBENT);
Eta = max(Eta.lb,
sum {p in ORDERS} spd[p] * (1 - On[p]),
max {optcut in OPTCUTS}
(cut_constant[optcut] + sum {w in LOCATIONS} cut_coefficient[w,optcut] * On[w])
);
solve with MILP / primalin logfreq=10000 target=(best_lower_bound) relobjgap=&relobjgap;
if _solution_status_ = 'INFEASIBLE' then do;
put 'ERROR: master is infeasible';
stop;
end;
Eta_sol = Eta.sol;
lower_bound = _OROPTMODEL_NUM_['BEST_BOUND'];
best_lower_bound = max(best_lower_bound, lower_bound);
SUPPORT = {w in LOCATIONS: On[w].sol > 0.5};
put SUPPORT=;
put 'Solving Subproblem...';
use problem Subproblem;
solve with LP / algorithm=ip crossover=0 logfreq=0;
if _solution_status_ = 'INFEASIBLE' then do;
put 'ERROR: subproblem is infeasible';
stop;
end;
else do;
if best_upper_bound > upper_bound then do;
best_upper_bound = upper_bound;
INCUMBENT = SUPPORT;
for {<p,w> in VALID_MAPS} Map_sol[p,w] = Map[p,w].sol;
create data incumbent from [LOCATION_ID]=INCUMBENT;
end;
if lower_bound < upper_bound then do;
if Eta_sol < SubproblemObjective then do;
/* generate new optimality cut */
num_optcuts = num_optcuts + 1;
cut_constant[num_optcuts] = sum {p in ORDERS} ORDERs_mapped_correctly[p].dual;
for {w in LOCATIONS}
cut_coefficient[w,num_optcuts] =
- sum{p in ORDERS} LOCATIONs_mapped_correctly[p,w].dual
- capacity[w] * CapacityCon[w].dual;
end;
end;
end;
put "BOUNDS: " best_lower_bound= best_upper_bound= gap=;
end;
put INCUMBENT=;
/* recover optimal solution */
for {w in LOCATIONS} On[w] = (w in INCUMBENT);
for {<p,w> in VALID_MAPS} Map[p,w] = Map_sol[p,w];
quit;
