Facility Location (milpsl3)
/*****************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: milpsl3 */
/* TITLE: Facility Location (milpsl3) */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: OPTMODEL, GPLOT */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 3 from the Mixed Integer Linear Programming */
/* Solver chapter of Mathematical Programming. */
/* */
/*****************************************************************/
title 'Facility Location Problem';
%let NumCustomers = 50;
%let NumSites = 10;
%let SiteCapacity = 35;
%let MaxDemand = 10;
%let xmax = 200;
%let ymax = 100;
%let seed = 938;
/* generate random customer locations */
data cdata(drop=i);
length name $8;
do i = 1 to &NumCustomers;
name = compress('C'||put(i,best.));
x = ranuni(&seed) * &xmax;
y = ranuni(&seed) * &ymax;
demand = ranuni(&seed) * &MaxDemand;
output;
end;
run;
/* generate random site locations and fixed charge */
data sdata(drop=i);
length name $8;
do i = 1 to &NumSites;
name = compress('SITE'||put(i,best.));
x = ranuni(&seed) * &xmax;
y = ranuni(&seed) * &ymax;
fixed_charge = 30 * (abs(&xmax/2-x) + abs(&ymax/2-y));
output;
end;
run;
proc optmodel;
set <str> CUSTOMERS;
set <str> SITES init {};
/* x and y coordinates of CUSTOMERS and SITES */
num x {CUSTOMERS union SITES};
num y {CUSTOMERS union SITES};
num demand {CUSTOMERS};
num fixed_charge {SITES};
/* distance from customer i to site j */
num dist {i in CUSTOMERS, j in SITES}
= sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2);
read data cdata into CUSTOMERS=[name] x y demand;
read data sdata into SITES=[name] x y fixed_charge;
var Assign {CUSTOMERS, SITES} binary;
var Build {SITES} binary;
min CostNoFixedCharge
= sum {i in CUSTOMERS, j in SITES} dist[i,j] * Assign[i,j];
min CostFixedCharge
= CostNoFixedCharge + sum {j in SITES} fixed_charge[j] * Build[j];
/* each customer assigned to exactly one site */
con assign_def {i in CUSTOMERS}:
sum {j in SITES} Assign[i,j] = 1;
/* if customer i assigned to site j, then facility must be built at j */
con link {i in CUSTOMERS, j in SITES}:
Assign[i,j] <= Build[j];
/* each site can handle at most &SiteCapacity demand */
con capacity {j in SITES}:
sum {i in CUSTOMERS} demand[i] * Assign[i,j] <=
&SiteCapacity * Build[j];
/* solve the MILP with no fixed charges */
solve obj CostNoFixedCharge with milp / logfreq = 500;
/* clean up the solution */
for {i in CUSTOMERS, j in SITES} Assign[i,j] = round(Assign[i,j]);
for {j in SITES} Build[j] = round(Build[j]);
call symput('varcostNo',put(CostNoFixedCharge,6.1));
/* create a data set for use by GPLOT */
create data CostNoFixedCharge_Data from
[customer site]={i in CUSTOMERS, j in SITES: Assign[i,j] = 1}
xi=x[i] yi=y[i] xj=x[j] yj=y[j];
/* solve the MILP, with fixed charges with warm start */
solve obj CostFixedCharge with milp / primalin logfreq = 500;
/* clean up the solution */
for {i in CUSTOMERS, j in SITES} Assign[i,j] = round(Assign[i,j]);
for {j in SITES} Build[j] = round(Build[j]);
num varcost = sum {i in CUSTOMERS, j in SITES} dist[i,j] * Assign[i,j].sol;
num fixcost = sum {j in SITES} fixed_charge[j] * Build[j].sol;
call symput('varcost', put(varcost,6.1));
call symput('fixcost', put(fixcost,5.1));
call symput('totalcost', put(CostFixedCharge,6.1));
/* create a data set for use by GPLOT */
create data CostFixedCharge_Data from
[customer site]={i in CUSTOMERS, j in SITES: Assign[i,j] = 1}
xi=x[i] yi=y[i] xj=x[j] yj=y[j];
quit;
title1 h=1.5 "Facility Location Problem";
title2 "TotalCost = &varcostNo (Variable = &varcostNo, Fixed = 0)";
data csdata;
set cdata(rename=(y=cy)) sdata(rename=(y=sy));
run;
/* create Annotate data set to draw line between customer and assigned site */
%annomac;
data anno(drop=xi yi xj yj);
%SYSTEM(2, 2, 2);
set CostNoFixedCharge_Data(keep=xi yi xj yj);
%LINE(xi, yi, xj, yj, *, 1, 1);
run;
proc gplot data=csdata anno=anno;
axis1 label=none order=(0 to &xmax by 10);
axis2 label=none order=(0 to &ymax by 10);
symbol1 value=dot interpol=none
pointlabel=("#name" nodropcollisions height=1) cv=black;
symbol2 value=diamond interpol=none
pointlabel=("#name" nodropcollisions color=blue height=1) cv=blue;
plot cy*x sy*x / overlay haxis=axis1 vaxis=axis2;
run;
quit;
title1 "Facility Location Problem";
title2 "TotalCost = &totalcost (Variable = &varcost, Fixed = &fixcost)";
/* create Annotate data set to draw line between customer and assigned site */
data anno(drop=xi yi xj yj);
%SYSTEM(2, 2, 2);
set CostFixedCharge_Data(keep=xi yi xj yj);
%LINE(xi, yi, xj, yj, *, 1, 1);
run;
proc gplot data=csdata anno=anno;
axis1 label=none order=(0 to &xmax by 10);
axis2 label=none order=(0 to &ymax by 10);
symbol1 value=dot interpol=none
pointlabel=("#name" nodropcollisions height=1) cv=black;
symbol2 value=diamond interpol=none
pointlabel=("#name" nodropcollisions color=blue height=1) cv=blue;
plot cy*x sy*x / overlay haxis=axis1 vaxis=axis2;
run;
quit;
%let NumCustomers = 45;
%let NumSites = 8;
%let SiteCapacity = 35;
%let MaxDemand = 10;
%let xmax = 200;
%let ymax = 100;
%let seed = 2345;
/* generate random customer locations */
data cdata(drop=i);
length name $8;
do i = 1 to &NumCustomers;
name = compress('C'||put(i,best.));
x = ranuni(&seed) * &xmax;
y = ranuni(&seed) * &ymax;
demand = ranuni(&seed) * &MaxDemand;
output;
end;
run;
/* generate random site locations and fixed charge */
data sdata(drop=i);
length name $8;
do i = 1 to &NumSites;
name = compress('SITE'||put(i,best.));
x = ranuni(&seed) * &xmax;
y = ranuni(&seed) * &ymax;
fixed_charge = (abs(&xmax/2-x) + abs(&ymax/2-y)) / 2;
output;
end;
run;
proc optmodel;
set <str> CUSTOMERS;
set <str> SITES init {};
/* x and y coordinates of CUSTOMERS and SITES */
num x {CUSTOMERS union SITES};
num y {CUSTOMERS union SITES};
num demand {CUSTOMERS};
num fixed_charge {SITES};
/* distance from customer i to site j */
num dist {i in CUSTOMERS, j in SITES}
= sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2);
read data cdata into CUSTOMERS=[name] x y demand;
read data sdata into SITES=[name] x y fixed_charge;
var Assign {CUSTOMERS, SITES} binary;
var Build {SITES} binary;
min CostFixedCharge
= sum {i in CUSTOMERS, j in SITES} dist[i,j] * Assign[i,j]
+ sum {j in SITES} fixed_charge[j] * Build[j];
/* each customer assigned to exactly one site */
con assign_def {i in CUSTOMERS}:
sum {j in SITES} Assign[i,j] = 1;
/* if customer i assigned to site j, then facility must be built at j */
con link {i in CUSTOMERS, j in SITES}:
Assign[i,j] <= Build[j];
/* each site can handle at most &SiteCapacity demand */
con capacity {j in SITES}:
sum {i in CUSTOMERS} demand[i] * Assign[i,j] <= &SiteCapacity * Build[j];
/* assign priority to Build variables (y) */
for{j in SITES} Build[j].priority=10;
/* solve the MILP with fixed charges, using branching priorities */
solve obj CostFixedCharge with milp / logfreq=1000;
quit;