The Decomposition Algorithm

Example 13.2 Generalized Assignment Problem

The generalized assignment problem (GAP) is that of finding a maximum profit assignment from $n$ tasks to $m$ machines such that each task is assigned to precisely one machine subject to capacity restrictions on the machines. With each possible assignment, associate a binary variable $x_{ij}$, which, if set to $1$, indicates that machine $i$ is assigned to task $j$. For ease of notation, define two index sets $M=\left\{ 1,\dots ,m\right\} $ and $N=\left\{ 1,\dots ,n\right\} $. A GAP can be formulated as a MILP as follows:

$\displaystyle  $
$\displaystyle \text {maximize}  $
$\displaystyle  \sum _{i \in M} \sum _{j \in N} p_{ij} x_{ij} $
$\displaystyle  $
$\displaystyle \text {subject to}  $
$\displaystyle  \sum _{i \in M} x_{ij}  $
$\displaystyle  = 1  $
$\displaystyle  $
$\displaystyle  j \in N  $
$\displaystyle  $
$\displaystyle  \text {(assignment)} $
$\displaystyle  $
$\displaystyle  $
$\displaystyle  \sum _{j \in N} w_{ij} x_{ij}  $
$\displaystyle  \leq b_ i  $
$\displaystyle  $
$\displaystyle  i \in M  $
$\displaystyle  $
$\displaystyle  \text {(knapsack)}  $
$\displaystyle  $
$\displaystyle  $
$\displaystyle  x_{ij}  $
$\displaystyle  \in \left\{ 0,1\right\}   $
$\displaystyle  $
$\displaystyle  i \in M, \  j \in N  $

In this formulation, constraints (assignment) ensure that each task is assigned to exactly one machine. Inequalities (knapsack) ensure that for each machine, the capacity restrictions are met.

Consider the following example taken from Koch et al. (2011) with $n=24$ tasks to be assigned to $m=8$ machines. The data set profit_data provides the profit for assigning a particular task to a particular machine:

%let NumTasks    = 24;
%let NumMachines = 8;

data profit_data;
   input p1-p&NumTasks;
   datalines;
25 23 20 16 19 22 20 16 15 22 15 21 20 23 20 22 19 25 25 24 21 17 23 17
16 19 22 22 19 23 17 24 15 24 18 19 20 24 25 25 19 24 18 21 16 25 15 20
20 18 23 23 23 17 19 16 24 24 17 23 19 22 23 25 23 18 19 24 20 17 23 23
16 16 15 23 15 15 25 22 17 20 19 16 17 17 20 17 17 18 16 18 15 25 22 17
17 23 21 20 24 22 25 17 22 20 16 22 21 23 24 15 22 25 18 19 19 17 22 23
24 21 23 17 21 19 19 17 18 24 15 15 17 18 15 24 19 21 23 24 17 20 16 21
18 21 22 23 22 15 18 15 21 22 15 23 21 25 25 23 20 16 25 17 15 15 18 16
19 24 18 17 21 18 24 25 18 23 21 15 24 23 18 18 23 23 16 20 20 19 25 21
;

The data set weight_data provides the amount of resources used by a particular task when assigned to a particular machine:

data weight_data;
   input w1-w&NumTasks;
   datalines;
 8 18 22  5 11 11 22 11 17 22 11 20 13 13  7 22 15 22 24  8  8 24 18  8
24 14 11 15 24  8 10 15 19 25  6 13 10 25 19 24 13 12  5 18 10 24  8  5
22 22 21 22 13 16 21  5 25 13 12  9 24  6 22 24 11 21 11 14 12 10 20  6
13  8 19 12 19 18 10 21  5  9 11  9 22  8 12 13  9 25 19 24 22  6 19 14
25 16 13  5 11  8  7  8 25 20 24 20 11  6 10 10  6 22 10 10 13 21  5 19
19 19  5 11 22 24 18 11  6 13 24 24 22  6 22  5 14  6 16 11  6  8 18 10
24 10  9 10  6 15  7 13 20  8  7  9 24  9 21  9 11 19 10  5 23 20  5 21
 6  9  9  5 12 10 16 15 19 18 20 18 16 21 11 12 22 16 21 25  7 14 16 10
;

Finally, the data set capacity_data provides the resource capacity for each machine:

data capacity_data;
   input b @@;
   datalines;
36 35 38 34 32 34 31 34
;

The following PROC OPTMODEL statements read in the data and define the necessary sets and parameters:

proc optmodel;
   /* declare index sets */
   set TASKS    = 1..&NumTasks;
   set MACHINES = 1..&NumMachines;

   /* declare parameters */ 
   num profit   {MACHINES, TASKS};
   num weight   {MACHINES, TASKS};
   num capacity {MACHINES};

   /* read data sets to populate data */
   read data profit_data   into [i=_n_] {j in TASKS} <profit[i,j]=col('p'||j)>;
   read data weight_data   into [i=_n_] {j in TASKS} <weight[i,j]=col('w'||j)>;
   read data capacity_data into [_n_] capacity=b;      

The following statements declare the optimization model:

   /* declare decision variables */
   var Assign {MACHINES, TASKS} binary;

   /* declare objective */
   max TotalProfit = 
      sum {i in MACHINES, j in TASKS} profit[i,j] * Assign[i,j];

   /* declare constraints */
   con AssignmentCon {j in TASKS}: 
      sum {i in MACHINES} Assign[i,j] = 1;

   con KnapsackCon {i in MACHINES}: 
      sum {j in TASKS} weight[i,j] * Assign[i,j] <= capacity[i];         

The following statements use two different decompositions to solve the problem. The first decomposition defines each assignment constraint as a block and uses the pure network simplex solver for the subproblem. The second decomposition defines each knapsack constraint as a block and uses the MILP solver for the subproblem.

   /* each assignment constraint defines a block */
   for{j in TASKS} 
      AssignmentCon[j].block = j;
      
   solve with milp / logfreq=1000 
      decomp        =() 
      decomp_subprob=(algorithm=nspure);

   /* each knapsack constraint defines a block */
   for{j in TASKS} 
      AssignmentCon[j].block = .;
   for{i in MACHINES} 
      KnapsackCon[i].block = i; 
      
   solve with milp / decomp=();
quit;

The solution summaries are displayed in Output 13.2.1.

Output 13.2.1: Solution Summaries

The OPTMODEL Procedure

Solution Summary
Solver MILP
Algorithm Decomposition
Objective Function TotalProfit
Solution Status Optimal within Relative Gap
Objective Value 563
Iterations 9652
Best Bound 563.05601358
Nodes 9079
   
Relative Gap 0.0000994814
Absolute Gap 0.0560135845
Primal Infeasibility 0
Bound Infeasibility 0
Integer Infeasibility 0

Solution Summary
Solver MILP
Algorithm Decomposition
Objective Function TotalProfit
Solution Status Optimal within Relative Gap
Objective Value 562.99999995
Iterations 23
Best Bound 563.0010001
Nodes 1
   
Relative Gap 1.7763306E-6
Absolute Gap 0.0010000759
Primal Infeasibility 2E-8
Bound Infeasibility 0
Integer Infeasibility 2E-8


The iteration log for both decompositions is shown in Output 13.2.2. This example is interesting because it shows the tradeoff between the strength of the relaxation and the difficulty of its resolution. In the first decomposition, the subproblems are totally unimodular and can be solved trivially. Consequently, each iteration of the decomposition algorithm is very fast. However, the bound obtained is as weak as the bound found in direct methods (the LP bound). The weaker bound leads to the need to enumerate more nodes overall. Alternatively, in the second decomposition, the subproblem is the knapsack problem, which is solved using MILP. In this case, the bound is much tighter and the problem solves in very few nodes. The tradeoff, of course, is that each iteration takes longer because solving the knapsack problem is not trivial. Another interesting aspect of this problem is that the subproblem coverage in the second decomposition is much smaller than that of the first decomposition. However, when dealing with MILP, it is not always the size of the coverage that determines the overall effectiveness of a particular choice of decomposition.

Output 13.2.2: Log

NOTE: There were 8 observations read from the data set WORK.PROFIT_DATA.              
NOTE: There were 8 observations read from the data set WORK.WEIGHT_DATA.              
NOTE: There were 8 observations read from the data set WORK.CAPACITY_DATA.            
NOTE: Problem generation will use 16 threads.                                         
NOTE: The problem has 192 variables (0 free, 0 fixed).                                
NOTE: The problem has 192 binary and 0 integer variables.                             
NOTE: The problem has 32 linear constraints (8 LE, 24 EQ, 0 GE, 0 range).             
NOTE: The problem has 384 linear constraint coefficients.                             
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).            
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 192 variables, 32 constraints, and 384 constraint     
      coefficients.                                                                   
NOTE: The MILP solver is called.                                                      
NOTE: The Decomposition algorithm is used.                                            
NOTE: The DECOMP method value USER is applied.                                        
NOTE: The subproblem solver chosen is an LP solver but at least one block has integer 
      variables.                                                                      
NOTE: The decomposition subproblems consist of 24 disjoint blocks.                    
NOTE: The decomposition subproblems cover 192 (100.00%) variables and 24 (75.00%)     
      constraints.                                                                    
NOTE: The deterministic parallel mode is enabled.                                     
NOTE: The Decomposition algorithm is using up to 16 threads.                          
      Iter         Best       Master         Best       LP       IP   CPU  Real       
                  Bound    Objective      Integer      Gap      Gap  Time  Time       
NOTE: Starting phase 1.                                                               
         1       0.0000       8.9248            . 8.92e+00        .     0     0       
         4       0.0000       0.0000            . 0.00e+00        .     0     0       
NOTE: Starting phase 2.                                                               
         5     589.9388     561.1588            .    4.88%        .     0     0       
         6     568.8833     568.5610            .    0.06%        .     0     0       
         7     568.6464     568.6464            .    0.00%        .     0     0       
         .     568.6464     568.6464     562.0000    0.00%    1.17%     0     0       
NOTE: Starting branch and bound.                                                      
         Node  Active   Sols         Best         Best      Gap     CPU    Real       
                                  Integer        Bound             Time    Time       
            0       1      1     562.0000     568.6464    1.17%       0       0       
         1000     838      1     562.0000     565.1733    0.56%       8       7       
         2000    1500      1     562.0000     564.5574    0.45%      16      14       
         3000    1930      1     562.0000     564.1714    0.38%      24      22       
         4000    2170      1     562.0000     563.9106    0.34%      33      30       
         5000    2174      1     562.0000     563.6909    0.30%      41      38       
         6000    1970      1     562.0000     563.5094    0.27%      51      45       
         7000    1586      1     562.0000     563.3436    0.24%      60      53       
         8000     992      1     562.0000     563.2000    0.21%      69      61       
         8447     635      2     563.0000     563.1429    0.03%      73      65       
         9000      82      2     563.0000     563.0657    0.01%      79      69       
         9078       4      2     563.0000     563.0560    0.01%      79      70       
NOTE: The Decomposition algorithm used 16 threads.                                    
NOTE: The Decomposition algorithm time is 70.34 seconds.                              
NOTE: Optimal within relative gap.                                                    
NOTE: Objective = 563.                                                                
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 192 variables, 32 constraints, and 384 constraint     
      coefficients.                                                                   
NOTE: The MILP solver is called.                                                      
NOTE: The Decomposition algorithm is used.                                            
NOTE: The DECOMP method value USER is applied.                                        
NOTE: The decomposition subproblems consist of 8 disjoint blocks.                     
NOTE: The decomposition subproblems cover 192 (100.00%) variables and 8 (25.00%)      
      constraints.                                                                    
NOTE: The deterministic parallel mode is enabled.                                     
NOTE: The Decomposition algorithm is using up to 16 threads.                          
      Iter         Best       Master         Best       LP       IP   CPU  Real       
                  Bound    Objective      Integer      Gap      Gap  Time  Time       
NOTE: Starting phase 1.                                                               
         1       0.0000      10.0000            . 1.00e+01        .     0     0       
         8       0.0000       0.0000            . 0.00e+00        .     0     0       
NOTE: Starting phase 2.                                                               
         9    1140.1732     496.8898            .   56.42%        .     0     0       
        10     810.8500     510.4200            .   37.05%        .     0     0       
        11     702.2908     521.9923            .   25.67%        .     0     0       
        12     641.1544     539.4899            .   15.86%        .     0     0       
        13     633.9641     543.8024            .   14.22%        .     0     0       
        14     632.1283     547.0870            .   13.45%        .     1     0       
        15     594.0000     550.5741            .    7.31%        .     1     0       
        16     588.1974     553.9880            .    5.82%        .     1     0       
        17     584.2143     555.2143            .    4.96%        .     1     0       
        19     571.0000     560.0000            .    1.93%        .     1     0       
        20     571.0000     562.0000            .    1.58%        .     1     0       
         .     571.0000     562.4000     555.0000    1.51%    2.80%     1     1       
        22     568.3333     563.3333     555.0000    0.88%    2.35%     1     1       
        23     564.0000     564.0000     555.0000    0.00%    1.60%     1     1       
         .     564.0000     564.0000     563.0000    0.00%    0.18%     1     1       
NOTE: The continuous bound was improved to 563.001 due to objective granularity.      
        23     563.0010     563.0010     563.0000    0.00%    0.00%     1     1       
         Node  Active   Sols         Best         Best      Gap     CPU    Real       
                                  Integer        Bound             Time    Time       
            0       0      2     563.0000     563.0010    0.00%       1       1       
NOTE: The Decomposition algorithm used 8 threads.                                     
NOTE: The Decomposition algorithm time is 1.37 seconds.                               
NOTE: Optimal within relative gap.                                                    
NOTE: Objective = 563.