The Decomposition Algorithm

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 one. If the donor and recipient are incompatible (because of blood types, tissue mismatch, and so on), the transplant cannot happen. Now suppose two donor-recipient pairs A and B are in this 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, shown graphically in Output 14.10.1.

More generally, an -way swap that involves donors and 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 exists if the donor from node is compatible with the recipient from node . Let define the set of nodes and define the set of arcs. The link weight, , 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 gives up a kidney if and only if the recipient from node 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 links. The kidney exchange problem is to find a maximum-weight node-disjoint union of short directed cycles.

Define an index set of candidate disjoint unions of short cycles (called matchings). Let be a binary variable, which, if set to , indicates that arc is in a matching . Let be a binary variable, which, if set to , indicates that node is covered by matching . In addition, let be a binary slack variable, which, if set to , indicates that node is not covered by any matching.

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

In this formulation, constraints (packing) ensure that each node is covered by at most one matching. Constraints (donate) and (receive) enforce the condition that if node is covered by matching , then the matching must use exactly one arc that leaves node (donate) and one arc that enters node (receive). Conversely, if node is not covered by matching , then no arcs that enter or leave node can be used by matching . Constraints (cardinality) enforce the condition that the number of arcs in matching must not exceed .

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 . Consider a decomposition by matching, where the packing constraints form the master problem and all other constraints form identical matching subproblems. As described in the section Special Case: Identical Blocks, 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 nodes with link probability 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;

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 node_packing {i in NODES}:
      sum {m in MATCHINGS} UseNode[i,m] + Slack[i] = 1;

   /* at most one donee 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 donee */
   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;

   /* decompose by matching (aggregate formulation) */
   for {i in NODES, m in MATCHINGS} donate[i,m].block = m;
   for {j in NODES, m in MATCHINGS} receive[j,m].block = m;
   for {m in MATCHINGS} cardinality[m].block = m;
   solve with milp / presolver=basic decomp=();

   /* 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 14.10.2.

The iteration log is displayed in Output 14.10.3.

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 8: The Network Solver). The optimal donor exchanges from the output data set Cycles are displayed in Output 14.10.4.