The OPTMILP Procedure

Example 6.3 Facility Location

This advanced example demonstrates how to warm start PROC OPTMILP by using the PRIMALIN= option. The model is constructed in PROC OPTMODEL and saved in an MPS-format data table for use in PROC OPTMILP. This problem can also be solved from within PROC OPTMODEL; see Chapter 14 for details.

Consider the classical facility location problem. Given a set L of customer locations and a set F of candidate facility sites, you must decide on which sites to build facilities and assign coverage of customer demand to these sites so as to minimize cost. All customer demand d Subscript i must be satisfied, and each facility has a demand capacity limit C. The total cost is the sum of the distances c Subscript i j between facility j and its assigned customer i, plus a fixed charge f Subscript j for building a facility at site j. Let y Subscript j Baseline equals 1 represent choosing site j to build a facility, and 0 otherwise. Also, let x Subscript i j Baseline equals 1 represent the assignment of customer i to facility j, and 0 otherwise. This model can be formulated as the following integer linear program:

StartLayout 1st Row 1st Column min 2nd Column sigma-summation Underscript i element-of upper L Endscripts sigma-summation Underscript j element-of upper F Endscripts c Subscript i j Baseline x Subscript i j plus sigma-summation Underscript j element-of upper F Endscripts f Subscript j Baseline y Subscript j 2nd Row 1st Column normal s period normal t period 2nd Column sigma-summation Underscript j element-of upper F Endscripts x Subscript i j 3rd Column equals 4th Column 1 5th Column for-all i element-of upper L 6th Column left-parenthesis normal a normal s normal s normal i normal g normal n normal bar normal d normal e normal f right-parenthesis 3rd Row 1st Column Blank 2nd Column x Subscript i j 3rd Column less-than-or-equal-to 4th Column y Subscript j 5th Column for-all i element-of upper L comma j element-of upper F 6th Column left-parenthesis normal l normal i normal n normal k right-parenthesis 4th Row 1st Column Blank 2nd Column sigma-summation Underscript i element-of upper L Endscripts d Subscript i Baseline x Subscript i j 3rd Column less-than-or-equal-to 4th Column upper C y Subscript j 5th Column for-all j element-of upper F 6th Column left-parenthesis normal c normal a normal p normal a normal c normal i normal t normal y right-parenthesis 5th Row 1st Column Blank 2nd Column x Subscript i j Baseline element-of StartSet 0 comma 1 EndSet 3rd Column Blank 4th Column Blank 5th Column for-all i element-of upper L comma j element-of upper F 6th Row 1st Column Blank 2nd Column y Subscript j Baseline element-of StartSet 0 comma 1 EndSet 3rd Column Blank 4th Column Blank 5th Column for-all j element-of upper F EndLayout

Constraint (assign_def) ensures that each customer is assigned to exactly one site. Constraint (link) forces a facility to be built if any customer has been assigned to that facility. Finally, constraint (capacity) enforces the capacity limit at each site.

Consider also a variation of this same problem where there is no cost for building a facility. This problem is typically easier to solve than the original problem. For this variant, let the objective be

StartLayout 1st Row 1st Column min 2nd Column sigma-summation Underscript i element-of upper L Endscripts sigma-summation Underscript j element-of upper F Endscripts c Subscript i j Baseline x Subscript i j EndLayout

First, construct a random instance of this problem by using the following DATA steps:

%let NumCustomers  = 30;
%let NumSites      = 6;
%let SiteCapacity  = 35;
%let MaxDemand     = 10;
%let xmax          = 200;
%let ymax          = 100;
%let seed          = 546;

/* generate random customer locations */
data mylib.cdata(drop=i);
   call streaminit(&seed);
   length name $8;
   do i = 1 to &NumCustomers;
      name = compress('C'||put(i,best.));
      x = rand('UNIFORM') * &xmax;
      y = rand('UNIFORM') * &ymax;
      demand = rand('UNIFORM') * &MaxDemand;
      output;
   end;
run;

/* generate random site locations and fixed charge */
data mylib.sdata(drop=i);
   call streaminit(&seed);
   length name $8;
   do i = 1 to &NumSites;
      name = compress('SITE'||put(i,best.));
      x = rand('UNIFORM') * &xmax;
      y = rand('UNIFORM') * &ymax;
      fixed_charge = 30 * (abs(&xmax/2-x) + abs(&ymax/2-y));
      output;
   end;
run;

The following PROC OPTMODEL statements generate the model and define both variants of the cost function:

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 mylib.cdata into CUSTOMERS=[name] x y demand;
   read data mylib.sdata into SITES=[name] x y fixed_charge;

   var Assign {CUSTOMERS, SITES} binary;
   var Build {SITES} binary;
   /* 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];

   min CostNoFixedCharge
      = sum {i in CUSTOMERS, j in SITES} dist[i,j] * Assign[i,j];

   save mps mylib.nofcdata;
   min CostFixedCharge
      = CostNoFixedCharge
        + sum {j in SITES} fixed_charge[j] * Build[j];

   save mps mylib.fcdata;

quit;

First solve the problem for the model with no fixed charge by using the following statements. The first PROC SQL call populates the macro variables varcostNo. This macro variable displays the objective value when the results are plotted. The second PROC SQL call generates a data table that is used to plot the results. The information printed in the log by PROC OPTMILP is displayed in Output 6.3.1.

proc optmilp data=mylib.nofcdata primalout=mylib.nofcout;
run;

proc sql noprint;
   select put(sum(_objcoef_ * _value_),6.1) into :varcostNo
   from mylib.nofcout;
quit;
proc sql;
   create table CostNoFixedCharge_Data as
   select
      scan(p._var_,2,'[],') as customer,
      scan(p._var_,3,'[],') as site,
      c.x as x1, c.y as y1, s.x as x2, s.y as y2
   from
      mylib.cdata as c,
      mylib.sdata as s,
      mylib.nofcout(where=(substr(_var_,1,6)='Assign' and
         round(_value_) = 1)) as p
   where calculated customer = c.name and calculated site = s.name;
quit;

Output 6.3.1: PROC OPTMILP Log for Facility Location with No Fixed Charges

NOTE: The problem nofcdata has 186 variables (186 binary, 0 integer, 0 free, 0  
      fixed).                                                                   
NOTE: The problem has 216 constraints (186 LE, 30 EQ, 0 GE, 0 range).           
NOTE: The problem has 726 constraint coefficients.                              
NOTE: The initial MILP heuristics are applied.                                  
NOTE: The MILP presolver value AUTOMATIC is applied.                            
NOTE: The MILP presolver removed 6 variables and 180 constraints.               
NOTE: The MILP presolver removed 366 constraint coefficients.                   
NOTE: The MILP presolver modified 0 constraint coefficients.                    
NOTE: The presolved problem has 180 variables, 36 constraints, and 360          
      constraint coefficients.                                                  
NOTE: The MILP solver is called.                                                
NOTE: The parallel Branch and Cut algorithm is used.                            
NOTE: The Branch and Cut algorithm is using up to 80 threads.                   
          Node   Active   Sols    BestInteger      BestBound      Gap    Time   
             0        1      1   1052.9414513    985.5035182    6.84%       0   
NOTE: The MILP presolver is applied again.                                      
             0        1      1   1052.9414513    985.5035182    6.84%       0   
             0        1      2   1045.3217618    985.5035182    6.07%       0   
             0        1      2   1045.3217618   1009.0043336    3.60%       0   
NOTE: The MILP presolver is applied again.                                      
             0        1      3   1043.4817129   1009.0043336    3.42%       0   
             0        1      3   1043.4817129   1009.0043336    3.42%       0   
             0        1      3   1043.4817129   1009.0043336    3.42%       0   
             0        1      5   1026.2068684   1009.0043336    1.70%       0   
             0        1      5   1026.2068684   1013.2272189    1.28%       0   
NOTE: The MILP presolver is applied again.                                      
             0        1      7   1026.0888959   1013.2272189    1.27%       0   
             0        1      7   1026.0888959   1013.2272189    1.27%       0   
             0        1      7   1026.0888959   1013.4298296    1.25%       0   
             0        1      7   1026.0888959   1015.1515946    1.08%       0   
             0        1      7   1026.0888959   1020.7316657    0.52%       0   
             0        1      8   1021.5193163   1021.5193163    0.00%       0   
NOTE: The MILP solver added 10 cuts with 66 cut coefficients at the root.       
NOTE: Optimal.                                                                  
NOTE: Objective = 1021.5193163.                                                 
NOTE: The Cloud Analytic Services server processed the request in 0.870211      
      seconds.                                                                  
NOTE: The data set MYLIB.NOFCOUT has 186 observations and 9 variables.          


Next, solve the fixed-charge model by using the following statements. Note that the solution to the model with no fixed charge is feasible for the fixed-charge model and should provide a good starting point for PROC OPTMILP. The PRIMALIN= option provides an incumbent solution ("warm start"). The two PROC SQL calls perform the same functions as in the case with no fixed charges. The results from this approach are shown in Output 6.3.2.

proc optmilp data=mylib.fcdata primalin=mylib.nofcout primalout=mylib.fcout;
run;

proc sql noprint;
   select put(sum(_objcoef_ * _value_), 6.1) into :varcost
   from mylib.fcout(where=(substr(_var_,1,6)='Assign'));
   select put(sum(_objcoef_ * _value_), 5.1) into :fixcost
   from mylib.fcout(where=(substr(_var_,1,5)='Build'));
   select put(sum(_objcoef_ * _value_), 6.1) into :totalcost
   from mylib.fcout;
quit;

proc sql;
   create table CostFixedCharge_Data as
   select
      scan(p._var_,2,'[],') as customer,
      scan(p._var_,3,'[],') as site,
      c.x as x1, c.y as y1, s.x as x2, s.y as y2
   from
      mylib.cdata as c,
      mylib.sdata as s,
      mylib.fcout(where=(substr(_var_,1,6)='Assign' and
         round(_value_) = 1)) as p
   where calculated customer = c.name and calculated site = s.name;
quit;

Output 6.3.2: PROC OPTMILP Log for Facility Location with Fixed Charges, Using Warm Start

NOTE: The problem fcdata has 186 variables (186 binary, 0 integer, 0 free, 0    
      fixed).                                                                   
NOTE: The problem has 216 constraints (186 LE, 30 EQ, 0 GE, 0 range).           
NOTE: The problem has 726 constraint coefficients.                              
NOTE: The initial MILP heuristics are applied.                                  
NOTE: The MILP presolver value AUTOMATIC is applied.                            
NOTE: The MILP presolver removed 0 variables and 0 constraints.                 
NOTE: The MILP presolver removed 0 constraint coefficients.                     
NOTE: The MILP presolver modified 0 constraint coefficients.                    
NOTE: The presolved problem has 186 variables, 216 constraints, and 726         
      constraint coefficients.                                                  
NOTE: The MILP solver is called.                                                
NOTE: The parallel Branch and Cut algorithm is used.                            
NOTE: The Branch and Cut algorithm is using up to 48 threads.                   
          Node   Active   Sols    BestInteger      BestBound      Gap    Time   
             0        1      4   9113.7292152   6814.6940297   33.74%       0   
             0        1      4   9113.7292152   9069.6292890    0.49%       0   
NOTE: The MILP presolver is applied again.                                      
             0        1      4   9113.7292152   9069.6292890    0.49%       0   
             0        1      4   9113.7292152   9069.6292890    0.49%       0   
             0        1      4   9113.7292152   9069.6292890    0.49%       0   
             0        1      7   9086.9946323   9069.6292890    0.19%       0   
             0        1      7   9086.9946323   9074.7236719    0.14%       0   
NOTE: The MILP presolver is applied again.                                      
             0        1      9   9086.8766598   9074.7236719    0.13%       0   
             0        1      9   9086.8766598   9074.7236719    0.13%       0   
             0        1     11   9082.3070802   9074.7236719    0.08%       0   
NOTE: The MILP solver added 28 cuts with 267 cut coefficients at the root.      
NOTE: Optimal within relative gap.                                              
NOTE: Objective = 9082.3070802.                                                 
NOTE: The Cloud Analytic Services server processed the request in 1.698669      
      seconds.                                                                  
NOTE: The data set MYLIB.FCOUT has 186 observations and 9 variables.            


The following two SAS programs produce a plot of the solutions for both variants of the model, using data tables produced by PROC SQL from the PRIMALOUT= data tables produced by PROC OPTMILP.

title1 "Facility Location Problem";
title2 "TotalCost = &varcostNo (Variable = &varcostNo, Fixed = 0)";
data mylib.csdata;
   set mylib.cdata(rename=(y=cy)) mylib.sdata(rename=(y=sy));
run;
/* create Annotate data set to draw line between customer and */
/* assigned site                                              */
data anno;
   retain function "line" drawspace "datavalue"
      linethickness 1 linecolor "black";
   set CostNoFixedCharge_Data(keep=x1 y1 x2 y2);
run;
proc sgplot data=mylib.csdata sganno=anno noautolegend;
   scatter x=x y=cy / datalabel=name datalabelattrs=(size=6pt)
      markerattrs=(symbol=circlefilled color=black size=6pt);
   scatter x=x y=sy / datalabel=name datalabelattrs=(size=6pt)
      markerattrs=(symbol=diamond color=blue size=6pt);
   xaxis display=(nolabel);
   yaxis display=(nolabel);
run;

The output from the first program appears in Output 6.3.3.

Output 6.3.3: Solution Plot for Facility Location with No Fixed Charges

Solution Plot for Facility Location with No Fixed Charges


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;
   retain function "line" drawspace "datavalue"
      linethickness 1 linecolor "black";
   set CostFixedCharge_Data(keep=x1 y1 x2 y2);
run;
proc sgplot data=mylib.csdata sganno=anno noautolegend;
   scatter x=x y=cy / datalabel=name datalabelattrs=(size=6pt)
      markerattrs=(symbol=circlefilled color=black size=6pt);
   scatter x=x y=sy / datalabel=name datalabelattrs=(size=6pt)
      markerattrs=(symbol=diamond color=blue size=6pt);
   xaxis display=(nolabel);
   yaxis display=(nolabel);
run;

The output from the second program appears in Output 6.3.4.

Output 6.3.4: Solution Plot for Facility Location with Fixed Charges

Solution Plot for Facility Location with Fixed Charges


The economic trade-off for the fixed-charge model forces you to build fewer sites and push more demand to each site.

Last updated: September 09, 2026