I have spent some time on this and can solve scenario 2 (optimal objective value is 347,566) in a few minutes by exploiting the problem structure. Explicitly, I relaxed the Map variables to be >= 0 instead of binary and then applied classical Benders decomposition, with a MILP master problem and an LP subproblem. The motivation is that the problem is easy once the On variables are fixed. See also this SAS Global Forum paper for a simplified approach that does not use the .DUAL suffix or exploit the independent blocks in the subproblem.
%let relobjgap = 1e-4;
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)};
/* 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.;
set LOCATIONS_p {ORDERS} init {};
for {<p,w> in VALID_MAPS} LOCATIONS_p[p] = LOCATIONS_p[p] union {w};
*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_p[p]} Map[p,w] = 1;
con Cover {p in ORDERS}:
/* sum{<(p),w> in VALID_MAPS} On[w] >= 1;*/
sum{w in LOCATIONS_p[p]} On[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);
*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 {ORDERS} >= 0;
num Eta_sol {ORDERS};
min MasterObjective = total_penalty + sum {p in ORDERS} Eta[p];
/* optimality cuts */
num num_optcuts {ORDERS} init 0;
set OPTCUTS {p in ORDERS} = 1..num_optcuts[p];
num cut_constant {p in ORDERS, OPTCUTS[p]};
num cut_coefficient {p in ORDERS, w in LOCATIONS_p[p], OPTCUTS[p]} init 0;
con OptimalityCut {p in ORDERS, optcut in OPTCUTS[p]}:
Eta[p] >= cut_constant[p,optcut] + sum {w in LOCATIONS_p[p]} cut_coefficient[p,w,optcut] * On[w];
/* omit feasibility cuts (assume Cover implies feasibility of subproblem) */
problem Master include
On Eta
MasterObjective
Cover
max_LOCATIONs_on
OptimalityCut;
/* declare subproblem over Map here (same as original but with On fixed) */
impvar SubproblemObjectiveBlock {p in ORDERS} = sum{w in LOCATIONS_p[p]} Map[p,w]*distance[p,w];
min SubproblemObjective = sum {p in ORDERS} SubproblemObjectiveBlock[p];
problem Subproblem include
Map
SubproblemObjective
ORDERs_mapped_correctly LOCATIONs_mapped_correctly;
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[*];*/
put (sum {p in ORDERS} num_optcuts[p])=;
for {w in LOCATIONS} On[w] = (w in INCUMBENT);
for {p in ORDERS} Eta[p] = max(Eta[p].lb,
max {optcut in OPTCUTS[p]}
(cut_constant[p,optcut] + sum {w in LOCATIONS_p[p]} cut_coefficient[p,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;
for {p in ORDERS} Eta_sol[p] = Eta[p].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;
for {p in ORDERS: Eta_sol[p] < SubproblemObjectiveBlock[p]} do;
/* generate new optimality cut */
num_optcuts[p] = num_optcuts[p] + 1;
cut_constant[p,num_optcuts[p]] = ORDERs_mapped_correctly[p].dual;
for {w in LOCATIONS_p[p]}
cut_coefficient[p,w,num_optcuts[p]] = -LOCATIONs_mapped_correctly[p,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;
... View more
