The Decomposition Algorithm

This formulation clearly has no decomposable structure. However, you can use a common modeling technique known as Lagrangian decomposition to bring the model into block-angular form. Lagrangian decomposition works by first partitioning the constraints into blocks. Then, each original variable is split into multiple copies of itself, one copy for each block in which the variable has a nonzero coefficient in the constraint matrix. Constraints are added to enforce the equality of each copy of the original variable. Then, the original constraints can be written in block-angular form by using the duplicate variables.

To apply Lagrangian decomposition to the resource allocation problem, define a set B of blocks and let define the set of start times for a given block b, such that . Given this partition of start times, let define the set of blocks in which task is scheduled to be running. Now, for each task , define duplicate variables for each . Let define the minimum block index for each class of variable that represents task i. The problem can now be modeled in block-angular form as follows:

In this formulation, the LinkDupVarsCon constraints ensure that the duplicate variables are equal to the original variables. Now, the five-task example has been transformed from a staircase structure to a block-angular structure:

To show how to apply Lagrangian decomposition in PROC OPTMODEL, consider the following data set TaskData from Caprara, Furini, and Malaguti (2010) which consists of tasks:

data TaskData;
   input profit weight start end;
   datalines;
100 74    1   12
 98 32    1    9
 73 27    1   22
 98 51    1   31
...
 23 40 2684 2689
 36 85 2685 2687
 65 44 2686 2689
 18 36 2687 2689
 88 57 2688 2689
;

The following PROC OPTMODEL statements read in the data and solve the original staircase formulation by calling the MILP solver directly:

%macro SetupData(task_data=, capacity=);
   set TASKS;
   num capacity=&capacity;
   num profit{TASKS}, weight{TASKS}, start{TASKS}, end{TASKS};

   read data &task_data into TASKS=[_n_] profit weight start end;
   /* the set of start times */

   set STARTS = setof{i in TASKS} start[i];
   /* the set of tasks i that are active at a given start time s */
   set TASKS_START{s in STARTS}
      = {i in TASKS: start[i] <= s < end[i]};
%mend SetupData;

%macro ResourceAllocation_Direct(task_data=, capacity=);
   proc optmodel;
      %SetupData(task_data=&task_data,capacity=&capacity);

      /* select task i to come online from period [start to end) */
      var x{TASKS} binary;

      /* maximize the total profit of running tasks */
      max TotalProfit = sum{i in TASKS} profit[i] * x[i];

      /* enforce that the shared resource capacity is not exceeded */
      con CapacityCon{s in STARTS}:
         sum{i in TASKS_START[s]} weight[i] * x[i] <= capacity;

      solve;
   quit;
%mend ResourceAllocation_Direct;

%ResourceAllocation_Direct(task_data=TaskData, capacity=100);

The problem summary and solution summary are displayed in Output 15.7.1.

The iteration log, which contains the problem statistics, the progress of the solution, and the optimal objective value, is shown in Output 15.7.2.

To transform this data into block-angular form, first sort the task data to help reduce the number of duplicate variables needed in the reformulation as follows:

proc sort data=TaskData;
   by start end;
run;

Then, create the partition of constraints into blocks of size block_size as follows:

%macro ResourceAllocation_Decomp(task_data=, capacity=, block_size=);
   proc optmodel;
      %SetupData(task_data=&task_data,capacity=&capacity);
      /* partition into blocks of size blocks_size */
      num block_size = &block_size;
      num num_blocks = ceil( card(TASKS) / block_size );
      set BLOCKS     = 1..num_blocks;

      /* the set of starts s for which task i is active */
      set STARTS_TASK{i in TASKS} = {s in STARTS: start[i] <= s < end[i]};

      /* partition the start times into blocks of size block_size */
      set STARTS_BLOCK{BLOCKS} init {};
      num block_id    init 1;
      num block_count init 0;
      for{s in STARTS} do;
         STARTS_BLOCK[block_id] = STARTS_BLOCK[block_id] union {s};
         block_count = block_count + 1;
         if(mod(block_count, block_size) = 0) then
             block_id = block_id + 1;
      end;

Then, the following PROC OPTMODEL statements define the block-angular formulation and solve the problem by using the decomposition algorithm, the PRESOLVER=BASIC option, and block_size=40. Because this reformulation is equivalent to the original staircase formulation, disabling some of the advanced presolver techniques ensures that the model maintains block-angularity.

      /* blocks in which task i is online */
      set BLOCKS_TASK{i in TASKS} =
         {b in BLOCKS: card(STARTS_BLOCK[b] inter STARTS_TASK[i]) > 0};

      /* minimum block id in which task i is online */
      num min_block{i in TASKS} = min{b in BLOCKS_TASK[i]} b;

      /* select task i to come online from period [start to end)
         in each block */
      var x{i in TASKS, b in BLOCKS_TASK[i]} binary;

      /* maximize the total profit of running tasks */
      max TotalProfit = sum{i in TASKS} profit[i] * x[i,min_block[i]];

      /* enforce that task selection is consistent across blocks */
      con LinkDupVarsCon{i in TASKS, b in BLOCKS_TASK[i] diff {min_block[i]}}:
         x[i,b] = x[i,min_block[i]];

      /* enforce that the shared resource capacity is not exceeded */
      con CapacityCon{b in BLOCKS, s in STARTS_BLOCK[b]}:
         sum{i in TASKS_START[s]} weight[i] * x[i,b] <= capacity;

      /* define blocks for decomposition algorithm */
      for{b in BLOCKS, s in STARTS_BLOCK[b]} CapacityCon[b,s].block = b;

      solve with milp / presolver=basic decomp;
   quit;
%mend ResourceAllocation_Decomp;

%ResourceAllocation_Decomp(task_data=TaskData, capacity=100, block_size=40);

The problem summary and solution summary are displayed in Output 15.7.3. Compared to the original formulation, the number of variables and constraints is increased by the number of duplicate variables.

The iteration log, which contains the problem statistics, the progress of the solution, and the optimal objective value, is shown in Output 15.7.4.

The decomposition algorithm solves the problem in fewer nodes due to the stronger bound obtained by the reformulation. However, it takes longer than the direct method to find a good feasible solution. The fact that the direct method seems to quickly find good feasible solutions but has weaker bounds motivates the use of a hybrid algorithm. In the macro %ResourceAllocation_Decomp, replace the statement,

      solve with milp / presolver=basic decomp;

with the following statements:

      solve with milp / relobjgap=0.1;
      solve with milp / presolver=basic primalin decomp;

These statements use the direct method with RELOBJGAP=0.1 to find a good starting solution and then use that result to seed the initial columns of the decomposition algorithm.

The solution summaries are displayed in Output 15.7.5.

The iteration log, which contains the problem statistics, the progress of the solution, and the optimal objective value, is shown in Output 15.7.6.

By using this hybrid method, you can take advantage of the direct method, which finds a good feasible solution quickly, and the strong bounds provided by the decomposition algorithm.

Alternatively, you can use the built-in hybrid method in one call to the solver by using the HYBRID=ON option in the DECOMP statement. This algorithm first processes the root node by using standard MILP techniques. It then proceeds with the decomposition algorithm to complete processing. In the macro %ResourceAllocation_Decomp, replace the statement

      solve with milp / presolver=basic decomp;

with the statement

      solve with milp / presolver=basic decomp=(hybrid=on);

The solution summary is displayed in Output 15.7.7.

The iteration log, which contains the problem statistics, the progress of the solution, and the optimal objective value, is shown in Output 15.7.8.

The reformulation of this resource allocation problem provides a nice example of the potential tradeoffs in modeling a problem for use with the decomposition algorithm. As seen in Example 15.2, the strength of the bound is an important factor in the overall performance of the algorithm, but it is not always correlated to the magnitude of the subproblem coverage. In this example, the block size determines the number of blocks. Moreover, it determines the number of linking variables that are needed in the reformulation. At one extreme, if the block size is set to be , then the number of blocks is 1, and the number of copies of original variables is 0. Using one block would be equivalent to the original staircase formulation and would not yield a model conducive to decomposition. As the number of blocks is increased, the number of linking variables increases (the size of the master problem), the strength of the decomposition bound decreases, and the difficulty of solving the subproblems decreases. In addition, as the number of blocks and their relative difficulty change, the efficient utilization of your machine’s parallel architecture can be affected.

The previous section used a block size of 40. The following statement calls the decomposition algorithm and uses a block size of 130:

%ResourceAllocation_Decomp(task_data=TaskData, capacity=100, block_size=130);

The solution summary is displayed in Output 15.7.9.

The iteration log, which contains the problem statistics, the progress of the solution, and the optimal objective value, is shown in Output 15.7.10.

This version of the model provides a stronger bound and requires a smaller branch-and-bound search tree to find an optimal solution.