The OPTMILP Procedure

Example 16.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 SAS data set for use in PROC OPTMILP. Note that this problem can also be solved from within PROC OPTMODEL; see Chapter 9 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 which sites to build facilities on and assign coverage of customer demand to these sites so as to minimize cost. All customer demand d_i must be satisfied, and each facility has a demand capacity limit c. The total cost is the sum of the distances c_{ij} between facility j and its assigned customer i, plus a fixed charge f_j for building a facility at site j. Let y_j = 1 represent choosing site j to build a facility, and 0 otherwise. Also, let x_{ij} = 1 represent the assignment of customer i to facility j, and 0 otherwise. This model can be formulated as the following integer linear program:

\min & \displaystyle\sum_{i \in l} \displaystyle\sum_{j \in f} c_{ij} x_{ij}    +...   ...\} & & & \forall i \in l, j \in f \    & y_{j} \in \{0,1\} & & & \forall j \in f

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.

Let us also consider 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

\min & \displaystyle\sum_{i \in l} \displaystyle\sum_{j \in f} c_{ij} x_{ij}

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

    %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;
 

In the following PROC OPTMODEL code, we generate the model and define both variants of the cost function:

    proc optmodel;                                                                                                                      
       set <str> CUSTOMERS;                                                                                                              
       set <str> SITES;                                                                                                                  
                                                                                                                                         
       /* 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;                                                                                                             
  
       /* 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 nofcdata;
  
       min CostFixedCharge                                                                                                               
           = CostNoFixedCharge 
             + sum {j in SITES} fixed_charge[j] * Build[j]; 
           
       save mps fcdata;
       
    quit;
 

We first solve the problem for the model with no fixed charge by using the following code. The first PROC SQL call populates the macro variables varcostNo. This macro variable is used to display the objective value when the results are plotted. The second PROC SQL call generates a data set which is used to plot the results. The information printed in the log by PROC OPTMILP is displayed in Output 16.3.1.

  
    proc optmilp data=nofcdata primalout=nofcout; 
    run; 
  
    proc sql noprint; 
       select put(sum(_objcoef_ * _value_),6.1) into :varcostNo 
       from 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 xi, c.y as yi, s.x as xj, s.y as yj 
       from 
          cdata as c, 
          sdata as s, 
          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 16.3.1: PROC OPTMILP Log for Facility Location with No Fixed Charges
 
NOTE: The problem nofcdata has 510 variables (510 binary, 0 integer, 0 free, 0
fixed).
NOTE: The problem has 560 constraints (510 LE, 50 EQ, 0 GE, 0 range).
NOTE: The problem has 2010 constraint coefficients.
NOTE: The OPTMILP presolver value AUTOMATIC is applied.
NOTE: The OPTMILP presolver removed 10 variables and 500 constraints.
NOTE: The OPTMILP presolver removed 1010 constraint coefficients.
NOTE: The OPTMILP presolver modified 0 constraint coefficients.
NOTE: The presolved problem has 500 variables, 60 constraints, and 1000
constraint coefficients.
NOTE: The MIXED INTEGER LINEAR solver is called.
Node Active Sols BestInteger BestBound Gap Time
0 1 0 . 961.2403449 . 0
0 1 2 966.4832160 966.4832160 0.00% 0
0 0 2 966.4832160 . 0.00% 0
NOTE: OPTMILP added 2 cuts with 94 cut coefficients at the root.
NOTE: Optimal.
NOTE: Objective = 966.483216.
NOTE: The data set WORK.NOFCOUT has 510 observations and 8 variables.
NOTE: PROCEDURE OPTMILP used (Total process time):
real time 0.10 seconds
cpu time 0.10 seconds
 
 
 



Next, we solve the fixed-charge model by using the following code. 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. We use the PRIMALIN= option to provide 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 16.3.2.

  
    proc optmilp data=fcdata primalin=nofcout; 
    run; 
  
    proc sql noprint; 
       select put(sum(_objcoef_ * _value_), 6.1) into :varcost 
       from fcout(where=(substr(_var_,1,6)='Assign')); 
       select put(sum(_objcoef_ * _value_), 5.1) into :fixcost 
       from fcout(where=(substr(_var_,1,5)='Build')); 
       select put(sum(_objcoef_ * _value_), 6.1) into :totalcost 
       from 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 xi, c.y as yi, s.x as xj, s.y as yj 
       from 
          cdata as c, 
          sdata as s, 
          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 16.3.2: PROC OPTMILP Log for Facility Location with Fixed Charges, Using Warm Start
 
NOTE: The problem fcdata has 510 variables (510 binary, 0 integer, 0 free, 0
fixed).
NOTE: The problem has 560 constraints (510 LE, 50 EQ, 0 GE, 0 range).
NOTE: The problem has 2010 constraint coefficients.
NOTE: The OPTMILP presolver value AUTOMATIC is applied.
NOTE: The OPTMILP presolver removed 0 variables and 0 constraints.
NOTE: The OPTMILP presolver removed 0 constraint coefficients.
NOTE: The OPTMILP presolver modified 0 constraint coefficients.
NOTE: The presolved problem has 510 variables, 560 constraints, and 2010
constraint coefficients.
NOTE: The MIXED INTEGER LINEAR solver is called.
Node Active Sols BestInteger BestBound Gap Time
0 1 1 18157.3952518 9946.2514269 82.56% 0
0 1 1 18157.3952518 9959.5657588 82.31% 0
0 1 3 10955.8156144 9970.3344501 9.88% 0
0 1 3 10955.8156144 9974.5475193 9.84% 0
0 1 3 10955.8156144 9974.9124195 9.83% 0
0 1 3 10955.8156144 9974.9216211 9.83% 0
0 1 3 10955.8156144 9975.7667099 9.82% 0
NOTE: OPTMILP added 16 cuts with 474 cut coefficients at the root.
21 18 4 10953.9474016 9975.7667099 9.81% 0
500 224 4 10953.9474016 10172.4076684 7.68% 3
1000 260 4 10953.9474016 10238.7585780 6.99% 7
1500 201 4 10953.9474016 10375.2261761 5.58% 10
2000 73 4 10953.9474016 10604.2672946 3.30% 13
2192 32 5 10950.0345546 10940.5483611 0.09% 13
2381 72 6 10949.9022616 10946.9424488 0.03% 14
2382 72 7 10949.9022584 10946.9424488 0.03% 14
2422 23 8 10948.4603436 10947.6159840 0.01% 14
NOTE: Optimal within relative gap.
NOTE: Objective = 10948.4603.
NOTE: The data set WORK.FCOUT has 510 observations and 8 variables.
NOTE: PROCEDURE OPTMILP used (Total process time):
real time 16.45 seconds
cpu time 16.45 seconds
 
 
 



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

Note: Execution of this code requires SAS/GRAPH software.

    title1 "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=0.7) cv=black;                                                                  
       symbol2 value=diamond interpol=none                                                                                                                                                      
          pointlabel=("#name" nodropcollisions color=blue height=0.7) cv=blue;                                                                   
       plot cy*x sy*x / overlay haxis=axis1 vaxis=axis2;                                                                                
    run;
    quit;
 

The output from the first program appears in Output 16.3.3.

Output 16.3.3: Solution Plot for Facility Location with No Fixed Charges
milpexfacloc4out.gif (16260 bytes)



    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=0.7) cv=black;                                                                    
       symbol2 value=diamond interpol=none                                                                                                                                                     
          pointlabel=("#name" nodropcollisions color=blue height=0.7) cv=blue;                                                                    
       plot cy*x sy*x / overlay haxis=axis1 vaxis=axis2;                                                                                
    run;
    quit;
 

The output from the second program appears in Output 16.3.4.

Output 16.3.4: Solution Plot for Facility Location with Fixed Charges
milpexfacloc5out.gif (17576 bytes)



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

Previous Page | Next Page | Top of Page