Resources

Cycle Detection for Kidney Donor Exchange (netsle02)

/************************************************************************/
/*                                                                      */
/*          S A S   S A M P L E   L I B R A R Y                         */
/*                                                                      */
/*    NAME: netsle02                                                    */
/*   TITLE: Cycle Detection for Kidney Donor Exchange (netsle02)        */
/* PRODUCT: OR                                                          */
/*  SYSTEM: ALL                                                         */
/*    KEYS: OR                                                          */
/*   PROCS: OPTMODEL, PRINT                                             */
/*    DATA:                                                             */
/*                                                                      */
/* SUPPORT:                             UPDATE:                         */
/*     REF:                                                             */
/*    MISC: Example 2 from the network solver documentation.            */
/*                                                                      */
/************************************************************************/

/* create random graph on n nodes with arc probability p
   and uniform(0,1) weight */
%let n = 100;
%let p = 0.02;
data LinkSetIn;
   call streaminit(1);
   do from = 0 to &n - 1;
      do to = 0 to &n - 1;
         if from eq to then continue;
         else if rand('UNIFORM') < &p then do;
            weight = rand('UNIFORM');
            output;
         end;
      end;
   end;
run;

%let max_length = 10;
proc optmodel;
   /* declare index sets and parameters, and read data */
   set <num,num> ARCS;
   num weight {ARCS};
   read data LinkSetIn into ARCS=[from to] weight;
   set<num,num,num> ID_ORDER_NODE;

   /* generate all cycles with 2 <= length <= max_length */
   solve with NETWORK /
      loglevel        = moderate
      graph_direction = directed
      links           = (include=ARCS)
      cycle           = (maxcycles=all minlength=2 maxlength=&max_length)
      out             = (cycles=ID_ORDER_NODE)
   ;

   /* extract <cid,from,to> triples from <cid,order,node> triples */
   set <num,num,num> ID_FROM_TO init {};
   num last init ., from, to;
   for {<cid,order,node> in ID_ORDER_NODE} do;
      from = last;
      to   = node;
      last = to;
      if order ne 1 then ID_FROM_TO = ID_FROM_TO union {<cid,from,to>};
   end;

   /* solve set packing problem to find maximum weight node-disjoint union
      of short directed cycles */
   set CYCLES = setof {<c,i,j> in ID_FROM_TO} c;
   set ARCS_c {c in CYCLES} = setof {<(c),i,j> in ID_FROM_TO} <i,j>;
   set NODES_c {c in CYCLES} = union {<i,j> in ARCS_c[c]} {i,j};
   set NODES = union {c in CYCLES} NODES_c[c];
   num cycle_weight {c in CYCLES} = sum {<i,j> in ARCS_c[c]} weight[i,j];

   /* UseCycle[c] = 1 if cycle c is used, 0 otherwise */
   var UseCycle {CYCLES} binary;

   /* declare objective */
   max TotalWeight
      = sum {c in CYCLES} cycle_weight[c] * UseCycle[c];

   /* each node appears in at most one cycle */
   con node_packing {i in NODES}:
      sum {c in CYCLES: i in NODES_c[c]} UseCycle[c] <= 1;

   /* call solver */
   solve with milp;

   /* output optimal solution */
   create data Solution from [c]={c in CYCLES: UseCycle[c].sol > 0.5}
      cycle_weight;
quit;
%put &_OROPTMODEL_;

proc print data=Solution noobs label;
run;