The Decomposition Algorithm

Example 15.11 Kidney Donor Exchange and METHOD=SET

This example looks at an application of integer programming to help create a kidney donor exchange. Suppose someone needs a kidney transplant and a family member is willing to donate a kidney. If the donor and recipient are incompatible (because of conflicting blood types, tissue mismatch, and so on), the transplant cannot proceed. Now suppose two donor-recipient pairs, A and B, are in this same situation, but donor A is compatible with recipient B and donor B is compatible with recipient A. Then two transplants can take place in a two-way swap, which is shown graphically in Figure 15.8.

Figure 15.8: Kidney Donor Exchange Two-Way Swap

Kidney Donor Exchange Two-Way Swap


More generally, an n-way swap that involves n donors and n recipients can be performed (Willingham 2009). To model this problem, define a directed graph as follows. Each node is an incompatible donor-recipient pair. Link $(i,j)$ exists if the donor from node i is compatible with the recipient from node j. Let N define the set of nodes and A define the set of arcs. The link weight, $w_{ij}$, is a measure of the quality of the match. By introducing dummy links whose weight is 0, you can also include altruistic donors who have no recipients or recipients who have no donors. The idea is to find a maximum-weight node-disjoint union of directed cycles. You want the union to be node-disjoint so that no kidney is donated more than once, and you want cycles so that the donor from node i gives up a kidney if and only if the recipient from node i receives a kidney.

Without any other constraints, the problem could be solved as a linear assignment problem. But doing so would allow arbitrarily long cycles in the solution. Because of practical considerations (such as travel) and to mitigate risk, each cycle must have no more than L links. The kidney exchange problem is to find a maximum-weight node-disjoint union of short directed cycles.

Define an index set $M = \{ 1,\dots ,|N|/2\} $ of candidate disjoint unions of short cycles (called matchings). Let $x_{ijm}$ be a binary variable, which, if set to 1, indicates that arc $(i,j)$ is in a matching m. Let $y_{im}$ be a binary variable that, if set to 1, indicates that node i is covered by matching m. In addition, let $s_ i$ be a binary slack variable that, if set to 1, indicates that node i is not covered by any matching.

The kidney donor exchange can be formulated as a MILP as follows:

\begin{align*} & \text {maximize} & \sum _{(i,j) \in A} \sum _{m \in M} w_{ij} x_{ijm}\\ & \text {subject to} & \sum _{m \in M} y_{im} + s_ i & = 1 & & i \in N & & \text {(Packing)}\\ & & \sum _{(i,j) \in A} x_{ijm} & = y_{im} & & i \in N, \ m \in M & & \text {(Donate)} \\ & & \sum _{(i,j) \in A} x_{ijm} & = y_{jm} & & j \in N, \ m \in M & & \text {(Receive)} \\ & & \sum _{(i,j) \in A} x_{ijm} & \leq L & & m \in M & & \text {(Cardinality)} \\ & & x_{ijm} & \in \left\{ 0,1\right\} & & (i,j) \in A, \ m \in M \\ & & y_{im} & \in \left\{ 0,1\right\} & & i \in N, \ m \in M \\ & & s_{i} & \in \left\{ 0,1\right\} & & i \in N \end{align*}

In this formulation, the Packing constraints ensure that each node is covered by at most one matching. The Donate and Receive constraints enforce the condition that if node i is covered by matching m, then the matching m must use exactly one arc that leaves node i (Donate) and one arc that enters node i (Receive). Conversely, if node i is not covered by matching m, then no arcs that enter or leave node i can be used by matching m. The Cardinality constraints enforce the condition that the number of arcs in matching m must not exceed L.

In this formulation, the matching identifier is arbitrary. Because it is not necessary to cover each incompatible donor-recipient pair (node), the Packing constraints can be modeled by using set partitioning constraints and the slack variable s. Consider a decomposition by matching, in which the Packing constraints form the master problem and all other constraints form identical matching subproblems. As described in the section Special Case: Identical Blocks and Ryan-Foster Branching, this is a situation in which an aggregate formulation and Ryan-Foster branching can greatly improve performance by reducing symmetry.

The following DATA step sets up the problem, first creating a random graph on n nodes with link probability p and Uniform(0,1) weight:

/* create random graph on n nodes with arc probability p
   and uniform(0,1) weight */
%let n = 100;
%let p = 0.02;
data ArcData;
   do i = 0 to &n - 1;
      do j = 0 to &n - 1;
         if i eq j then continue;
         else if ranuni(1) < &p then do;
            weight = ranuni(2);
            output;
         end;
      end;
   end;
run;

In this case, you can specify METHOD=SET and let the decomposition algorithm automatically detect the set partitioning master constraints (Packing) and each independent matching subproblem. The following PROC OPTMODEL statements read in the data, declare the optimization model, and use the decomposition algorithm to solve it:

%let max_length = 10;
proc optmodel;
   set <num,num> ARCS;
   num weight {ARCS};
   read data ArcData into ARCS=[i j] weight;
   print weight;
   set NODES = union {<i,j> in ARCS} {i,j};
   set MATCHINGS = 1..card(NODES)/2;

   /* UseNode[i,m] = 1 if node i is used in matching m, 0 otherwise */
   var UseNode {NODES, MATCHINGS} binary;

   /* UseArc[i,j,m] = 1 if arc (i,j) is used in matching m, 0 otherwise */
   var UseArc {ARCS, MATCHINGS} binary;

   /* maximize total weight of arcs used */
   max TotalWeight
      = sum {<i,j> in ARCS, m in MATCHINGS} weight[i,j] * UseArc[i,j,m];

   /* each node appears in at most one matching */
   /* rewrite as set partitioning (so decomp uses identical blocks)
      sum{} x <= 1 => sum{} x + s = 1, s >= 0 with no associated cost */
   var Slack {NODES} binary;
   con Packing {i in NODES}:
      sum {m in MATCHINGS} UseNode[i,m] + Slack[i] = 1;

   /* at most one recipient for each donor */
   con Donate {i in NODES, m in MATCHINGS}:
      sum {<(i),j> in ARCS} UseArc[i,j,m] = UseNode[i,m];

   /* at most one donor for each recipient */
   con Receive {j in NODES, m in MATCHINGS}:
      sum {<i,(j)> in ARCS} UseArc[i,j,m] = UseNode[j,m];

   /* exclude long matchings */
   con Cardinality {m in MATCHINGS}:
      sum {<i,j> in ARCS} UseArc[i,j,m] <= &max_length;

   /* automatically decompose using METHOD=SET */
   solve with milp / presolver=basic decomp=(method=set);

   /* save solution to a data set */
   create data Solution from
      [m i j]={m in MATCHINGS, <i,j> in ARCS: UseArc[i,j,m].sol > 0.5}
      weight[i,j];
quit;

In this case, the PRESOLVER=BASIC option ensures that the model maintains its specified symmetry, enabling the algorithm to use the aggregate formulation and Ryan-Foster branching. The solution summary is displayed in Output 15.11.1.

Output 15.11.1: Solution Summary

The OPTMODEL Procedure

Solution Summary
Solver MILP
Algorithm Decomposition
Objective Function TotalWeight
Solution Status Optimal within Relative Gap
Objective Value 26.020287142
   
Relative Gap 2.1560594E-6
Absolute Gap 0.0000561014
Primal Infeasibility 2.220446E-15
Bound Infeasibility 1.110223E-15
Integer Infeasibility 1.110223E-15
   
Best Bound 26.020343243
Nodes 19
Iterations 124
Presolve Time 0.02
Solution Time 20.48



The iteration log is displayed in Output 15.11.2.

Output 15.11.2: Log

NOTE: There were 194 observations read from the data set WORK.ARCDATA.                          
NOTE: Problem generation will use 4 threads.                                                    
NOTE: The problem has 14065 variables (0 free, 0 fixed).                                        
NOTE: The problem has 14065 binary and 0 integer variables.                                     
NOTE: The problem has 9457 linear constraints (48 LE, 9409 EQ, 0 GE, 0 range).                  
NOTE: The problem has 42001 linear constraint coefficients.                                     
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).                      
NOTE: The MILP presolver value BASIC is applied.                                                
NOTE: The MILP presolver removed 4786 variables and 3298 constraints.                           
NOTE: The MILP presolver removed 14290 constraint coefficients.                                 
NOTE: The MILP presolver modified 0 constraint coefficients.                                    
NOTE: The presolved problem has 9279 variables, 6159 constraints, and 27711 constraint          
      coefficients.                                                                             
NOTE: The MILP solver is called.                                                                
NOTE: The Decomposition algorithm is used.                                                      
NOTE: The Decomposition algorithm is executing in single-machine mode.                          
NOTE: The DECOMP method value SET is applied.                                                   
NOTE: All blocks are identical and the master model is set partitioning.                        
NOTE: The Decomposition algorithm is using an aggregate formulation and Ryan-Foster branching.  
NOTE: The number of block threads has been reduced to 1 threads.                                
NOTE: The problem has a decomposable structure with 48 blocks. The largest block covers 2.062%  
      of the constraints in the problem.                                                        
NOTE: The decomposition subproblems cover 9216 (99.32%) variables and 6096 (98.98%) constraints.
NOTE: The deterministic parallel mode is enabled.                                               
NOTE: The Decomposition algorithm is using up to 4 threads.                                     
      Iter         Best       Master         Best       LP       IP  CPU Real                   
                  Bound    Objective      Integer      Gap      Gap Time Time                   
         .     390.3714       9.2503       9.2503   97.63%   97.63%    0    0                   
         1     388.4993       9.2503       9.2503   97.62%   97.62%    0    0                   
         2     382.6496      10.9240      10.9240   97.15%   97.15%    0    0                   
         3     359.6462      10.9240      10.9240   96.96%   96.96%    0    0                   
         5     359.6462      14.1405      14.1405   96.07%   96.07%    0    0                   
         6     355.2267      21.3961      21.3961   93.98%   93.98%    0    0                   
         7     348.2985      21.3961      21.3961   93.86%   93.86%    0    0                   
         8     330.2734      21.3961      21.3961   93.52%   93.52%    0    0                   
         9     229.6871      22.1269      21.3961   90.37%   90.68%    1    1                   
         .     229.6871      22.1269      21.3961   90.37%   90.68%    1    1                   
        10     229.6871      22.1269      21.3961   90.37%   90.68%    1    1                   
        12     201.2476      23.8995      21.3961   88.12%   89.37%    1    2                   
        16     176.4324      24.5325      21.3961   86.10%   87.87%    2    2                   
        17     159.1914      24.5842      21.3961   84.56%   86.56%    2    2                   
        18     151.9138      24.8491      21.3961   83.64%   85.92%    2    2                   
        19     135.4435      25.0568      21.3961   81.50%   84.20%    2    2                   
         .     135.4435      25.0568      21.9864   81.50%   83.77%    2    2                   
        20     135.4435      25.0568      21.9864   81.50%   83.77%    2    2                   
        22     118.8012      25.2761      22.1565   78.72%   81.35%    2    2                   
        25     109.0404      25.9874      22.1565   76.17%   79.68%    2    2                   
        26     102.9869      26.0423      22.1565   74.71%   78.49%    2    3                   
        27      99.9984      26.0586      22.1565   73.94%   77.84%    2    3                   
        28      65.7400      26.2351      22.1565   60.09%   66.30%    3    3                   
        30      47.2447      26.4827      22.1565   43.95%   53.10%    3    3                   
        36      41.4607      26.7010      22.1565   35.60%   46.56%    4    4                   
        39      27.6471      26.7804      22.1565    3.13%   19.86%    4    4                   
         .      27.6471      26.7804      24.4269    3.13%   11.65%    5    4                   
        40      26.7804      26.7804      24.4269    0.00%    8.79%    5    5                   
         .      26.7804      26.7804      24.4269    0.00%    8.79%    5    5                   
NOTE: Starting branch and bound.                                                                
         Node  Active   Sols         Best         Best      Gap    CPU   Real                   
                                  Integer        Bound            Time   Time                   
            0       1      9      24.4269      26.7804    8.79%      5      5                   
            4       6     10      25.9955      26.5360    2.04%     10      9                   
            8       4     11      26.0203      26.2556    0.90%     16     14                   
           10       4     11      26.0203      26.2478    0.87%     17     15                   
           18       0     11      26.0203      26.0203    0.00%     23     20                   
NOTE: The Decomposition algorithm used 4 threads.                                               
NOTE: The Decomposition algorithm time is 20.40 seconds.                                        
NOTE: Optimal within relative gap.                                                              
NOTE: Objective = 26.020287142.                                                                 
NOTE: The data set WORK.SOLUTION has 42 observations and 4 variables.                           



The solution is a set of arcs that define a union of short directed cycles (matchings). The following call to PROC OPTNET extracts the corresponding cycles from the list of arcs and outputs them to the data set Cycles:

proc optnet
   direction  = directed
   data_links = Solution;
   data_links_var
      from    = i
      to      = j;
   cycle
      mode    = all_cycles
      out     = Cycles;
run;

For more information about PROC OPTNET, see SAS/OR User's Guide: Network Optimization Algorithms. Alternatively, you can extract the cycles by using the SOLVE WITH NETWORK statement in PROC OPTMODEL (see Chapter 9: The Network Solver). The optimal donor exchanges from the output data set Cycles are displayed in Figure 15.9.

Figure 15.9: Optimal Donor Exchanges

order node
1 5
2 19
3 56
4 12
5 33
6 70
7 63
8 43
9 15
10 5

order node
1 13
2 74
3 65
4 41
5 59
6 50
7 49
8 98
9 13

order node
1 16
2 91
3 17
4 57
5 87
6 72
7 64
8 22
9 88
10 16

order node
1 8
2 32
3 79
4 71
5 69
6 26
7 9
8 18
9 95
10 35
11 8

order node
1 52
2 77
3 94
4 81
5 52

order node
1 24
2 92
3 24