Solving Sudoku with PROC OPTMODEL using MILP solver

From sasCommunity

/* Solving Sudoku puzzles with PROC OPTMODEL using MILP solver
/* Author: Rob Pratt, SAS/OR

/* Place a number into each box so that each row across, each 
column down, and each small 9-box square within the larger square
(there are 9 of these) will contain every number from 1 through 9. */

/* specify input matrix in dense format, with missing values for empty boxes */
data indata;
   input C1-C9;
   if (_n_=1 | _n_=4 | _n_=7) then put '-------------------------';
   put   '| ' C1 C2 C3 '| ' C4 C5 C6 '| ' C7 C8 C9 '|';
   if (_n_=9) then put '-------------------------'/;
   datalines;
1 . . . . 7 . 9 .
. 3 . . 2 . . . 8
. . 9 6 . . 5 . .
. . 5 3 . . 9 . .
. 1 . . 8 . . . 2
6 . . . . 4 . . .
3 . . . . . . 1 .
. 4 . . . . . . 7
. . 7 . . . 3 . .
;
run;

proc optmodel presolver=none;
   number n = 9;
   set DIGITS = 1..n;
   set ROWS = 1..n;
   set COLS = 1..n;
   number f {ROWS,COLS};

   /* read two-dimensional data */
   read data indata into [i = _n_] {j in COLS} <f[i,j] = col("c"||j)>;
   print f;

   /* decision variables: x[i,j,k] = 1 if cell (i,j) contains digit k, 0 otherwise */
   var x {ROWS,COLS,DIGITS} binary;

   /* dummy objective */
   max z = 0;

   /* digit constraints */
   con digit {i in ROWS, j in COLS}:
       sum {k in DIGITS} x[i,j,k] = 1;
   /* row constraints */
   con row_con {i in ROWS, k in DIGITS}:
       sum {j in COLS} x[i,j,k] = 1;
   /* column constraints */
   con col_con {j in COLS, k in DIGITS}:
       sum {i in ROWS} x[i,j,k] = 1;
   /* 9-box square constraints */
   con box_con {ri in 1..n by sqrt(n), rj in 1..n by sqrt(n), k in DIGITS}:
       sum {i in ri..ri+sqrt(n)-1, j in rj..rj+sqrt(n)-1} x[i,j,k] = 1;

   /* fix initial values */
   for {i in ROWS, j in COLS: f[i,j] ne .}
      fix x[i,j,f[i,j]] = 1;

   *expand;
   *solve;
   solve with milp / heuristics=aggressive cuts=none varsel=maxinfeas;
   put (symget("_OROPTMODEL_"));

   print _page_;
   print "Nonzero variables" {i in ROWS, j in COLS, k in DIGITS: round(x[i,j,k]) = 1} x[i,j,k];
   print _page_;
   print {i in ROWS, j in COLS} (sum {k in DIGITS} k * x[i,j,k]);

   create data outdata(drop=i) from [i]=ROWS {j in COLS} <col("C"||j) = (sum {k in DIGITS} k * x[i,j,k])>;

   /* to check for uniqueness, exclude this solution and solve again
   con exclude_con:
      sum {i in ROWS, j in COLS, k in DIGITS: round(x[i,j,k].sol) = 0} x[i,j,k]
    + sum {i in ROWS, j in COLS, k in DIGITS: round(x[i,j,k].sol) = 1} (1 - x[i,j,k])
   >= 1;
   print _page_;
   solve;
   put (symget("_OROPTMODEL_")); */
quit;

data _null_;
   set outdata;
   if (_n_=1 | _n_=4 | _n_=7) then put '-------------------------';
   put   '| ' C1 C2 C3 '| ' C4 C5 C6 '| ' C7 C8 C9 '|';
   if (_n_=9) then put '-------------------------'/;
run;