The OPTMODEL Procedure

Example 5.7 Sparse Modeling

(View the complete code for this example.)

This example demonstrates how to rewrite certain models for more efficient processing. Sometimes optimization models that run out of memory during problem generation can be rewritten to take advantage of sparsity to use memory more efficiently. This often occurs when a large array is modeled in a dense format but most of its entries are zeros. Usually, the array provides problem coefficients or it contains optimization variables.

The model for this example solves the facility location problem that is described in Example 9.3. This example is concerned with the resources that are required for PROC OPTMODEL problem generation and solver initialization. So the size of the problem has been increased, but the model has also been modified to make it easier to solve. In order to handle the larger problem size, the model eliminates a large number of the potential assignments of customers to facilities based on distance, making the problem sparse.

The following code generates a random instance of the facility location problem:

%let NumCustomers  = 1500;
%let NumSites      = 250;
%let SiteCapacity  = 50;
%let MaxDemand     = 10;
%let xmax          = 200;
%let ymax          = 100;
%let seed          = 938;

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

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

The following code uses a dense version of the facility location model. This model is equivalent to the model from Example 9.3 except for the added constraint distance_at_most_30. This constraint eliminates from consideration the assignment of customers to facilities over long distances by forcing the corresponding Assign variables to 0.

proc optmodel printlevel=2;
   profile on percent=0.1;
   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];

   /* do not assign customer to site more than 30 units away */
   con distance_at_most_30 {i in CUSTOMERS, j in SITES: dist[i,j] > 30}:
      Assign[i,j] = 0;

   /* solve the MILP */
   solve with milp/timetype=real;

quit;

If you inspect the log after running the preceding code, then you will see that the MILP presolver has pruned down the problem size considerably. If you also run the code with the SAS option FULLSTIMER enabled on a 64-bit system, then you will notice that about 1.1GB of memory is required for the PROC OPTMODEL step when you are running on a single CPU.

The solution and timing results for the dense model are shown in Output 5.7.1. The PROFILE ON statement requests further timing details, including evaluation time for declarations used by problem generation.

Output 5.7.1: Dense Model Results

Solution Summary
SolverMILP
AlgorithmBranch and Cut
Objective FunctionCostFixedCharge
Solution StatusOptimal within Relative Gap
Objective Value18001.140353
  
Relative Gap0.0000287924
Absolute Gap0.5182815753
Primal Infeasibility6.57252E-13
Bound Infeasibility1.042499E-13
Integer Infeasibility1.101341E-13
  
Best Bound18000.622072
Nodes12
Solutions Found5
Iterations33823
Presolve Time20.77
Solution Time268.52

Procedure Task Timing
TaskTime
(sec.)
% Time
Problem Generation25.928.61%
Solver Initialization5.511.83%
Code Generation0.160.05%
Presolve20.776.89%
Root Node Processing212.1670.44%
Branch And Cut18.656.19%
Synchronization2.290.76%
Idle10.643.53%
Other Tasks4.021.34%
Solver Postprocessing1.060.35%

Profile Information
ItemLineCol.Execution
Count
Net Time
(sec.)
Wait Time
(sec.)
% Total
Time
SOLVE362641278.576.7271.5%
Constraint distance_at_most_3036228 55.720.0014.3%
Constraint link36138 29.030.007.5%
Constraint assign_def36098 10.037.202.6%
Constraint capacity36178 5.230.171.3%
Number dist35948 5.030.001.3%
Min CostNoFixedCharge36038 3.230.000.8%
Var Assign36008 2.580.000.7%
Other profiled items   0.060.000.0%

Note:Total profiled time is 389.48 seconds.



The best approach for reducing the memory requirements is to eliminate the Assign variables that are always going to be 0. This is accomplished in the following sparse version of the code. Instead of indexing Assign over the crossproduct of CUSTOMERS and SITES, now the code defines a new set of pairs that satisfy the distance requirement, CUSTOMERS_SITES. This set replaces the constraint distance_at_most_30 in the dense model. The objective and constraints have been modified to use the new indexing scheme, with implicit set slicing (as described in the section More on Index Sets) for constraints assign_def and capacity.

proc optmodel printlevel=2;
   profile on percent=0.1;
   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;

   set CUSTOMERS_SITES = {i in CUSTOMERS, j in SITES: dist[i,j] <= 30};
   var Assign {CUSTOMERS_SITES} binary;
   var Build {SITES} binary;

   min CostNoFixedCharge
       = sum {<i,j> in CUSTOMERS_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 {<(i),j> in CUSTOMERS_SITES} Assign[i,j] = 1;

   /* if customer i assigned to site j, then facility must be built at j */
   con link {<i,j> in CUSTOMERS_SITES}:
      Assign[i,j] <= Build[j];

   /* each site can handle at most &SiteCapacity demand */
   con capacity {j in SITES}:
      sum {<i,(j)> in CUSTOMERS_SITES} demand[i] * Assign[i,j] <=
         &SiteCapacity * Build[j];

   /* solve the MILP */
   solve with milp/timetype=real;

quit;

The log from running the preceding code shows that the MILP presolver does not find anything to improve with this version of the model. On a 64-bit system, the FULLSTIMER option shows that memory requirements have been reduced to about 320MB when you are running on a single CPU, less than half the requirements of the previous model.

The solution and timing results for the sparse model are shown in Output 5.7.2. Note that the dense model (Output 5.7.1) and the sparse model (Output 5.7.2) are equivalent after presolver processing and generate the same result using similar amounts of solver time. On the other hand, problem generation time is significantly reduced as are other times including presolve time. Both models used the solver option TIMETYPE=REAL so that all times are reported in seconds of real time. You can see from the "Profile Information" table that the overhead associated with problem declarations has been significantly reduced.

Output 5.7.2: Sparse Model Results

Solution Summary
SolverMILP
AlgorithmBranch and Cut
Objective FunctionCostFixedCharge
Solution StatusOptimal within Relative Gap
Objective Value18001.140353
  
Relative Gap0.0000287924
Absolute Gap0.5182815753
Primal Infeasibility6.57252E-13
Bound Infeasibility1.042499E-13
Integer Infeasibility1.101341E-13
  
Best Bound18000.622072
Nodes12
Solutions Found7
Iterations33823
Presolve Time37.35
Solution Time282.70

Procedure Task Timing
TaskTime
(sec.)
% Time
Problem Generation10.123.45%
Solver Initialization0.630.22%
Code Generation0.020.01%
Presolve37.3512.72%
Root Node Processing209.2171.26%
Branch And Cut20.126.85%
Synchronization1.970.67%
Idle9.953.39%
Other Tasks4.111.40%
Solver Postprocessing0.110.04%

Profile Information
ItemLineCol.Execution
Count
Net Time
(sec.)
Wait Time
(sec.)
% Total
Time
SOLVE367541284.110.3995.0%
Number dist36468 4.810.001.6%
Set CUSTOMERS_SITES36528 3.980.001.3%
Constraint link36668 2.540.000.8%
Constraint assign_def36628 1.800.000.6%
Constraint capacity36708 0.980.050.3%
Min CostNoFixedCharge36568 0.370.000.1%
Other profiled items   0.340.000.1%

Note:Total profiled time is 298.94 seconds.



Last updated: October 07, 2022