Sudoku can be solved using a 3-index integer linear program (ILP) formulation using only binary variables, (google "ILP Sudoku" and read the first three hits).
But can it be solved with a two-index formulation using integer values [1,9]? (Yes.)
MOSS is based on finding a feasible solution to the following system of equations:
(a) Minimize SUM(x_ij)
With the following constraints:
(b) SUM(x_ij)==45 across j for all i
(c) SUM(x_ij)==45 across i for all j
(d) SUM(x_ij)==45 for x_i_j within each nonet
(e) ABS(x_ij - x_i_notj)>=1 for all i, j
(f) ABS(x_ij - x_noti_j)>=1 for all i, j
Where x is the Sudoku board, and, x_ij are cell-values within the board.
MOSS solves Sudoku entirely using an ILP solver (e.g. COIN-OR CBC, CPLEX, and Gurobi). Heuristics are restricted to use only within the solver, (i.e. only those heuristics integrated into COIN-OR CBC are used). Heuristics designed specifically for Sudoku are not used.
MOSS utilizes the absolute value constraint to enforce the all-different rule necesary for a Sudoku solution. For details of the absolute value constraint, see Gurobi Modeling, "Non Convex Case", page 9: https://www.gurobi.com/pdfs/user-events/2017-frankfurt/Modeling-2.pdf.
For details of the formulation, see Appendix A: Mathematical Formulation.
As of 12/20/2019 MOSS satisfies all constraints. Now witness the firepower of this fully armed and operational battle station!
gentle_board_9 = np.array([
... [0, 0, 0, 2, 6, 0, 7, 0, 1],
... [6, 8, 0, 0, 7, 0, 0, 9, 0],
... [1, 9, 0, 0, 0, 4, 5, 0, 0],
... [8, 2, 0, 1, 0, 0, 0, 4, 0],
... [0, 0, 4, 6, 0, 2, 9, 0, 0],
... [0, 5, 0, 0, 0, 3, 0, 2, 8],
... [0, 0, 9, 3, 0, 0, 0, 7, 4],
... [0, 4, 0, 0, 5, 0, 0, 3, 6],
... [7, 0, 3, 0, 1, 8, 0, 0, 0]
...])
gb9_solved = solve_board(gentle_board_9, max_solve_time=600000)
...
(1) Initializing optimization model...
(2) Creating objective variables...
(3) Setting constraints...
0
Success!
Board of size 9 solved in 0 seconds, using 134 simplex iterations.
print(gb9_solved['solution'])
[[4 3 5 2 6 9 7 8 1]
[6 8 2 5 7 1 4 9 3]
[1 9 7 8 3 4 5 6 2]
[8 2 6 1 9 5 3 4 7]
[3 7 4 6 8 2 9 1 5]
[9 5 1 7 4 3 6 2 8]
[5 1 9 3 2 6 8 7 4]
[2 4 8 9 5 7 1 3 6]
[7 6 3 4 1 8 2 5 9]]
Thanks to Seth W for the inspiration for this toy problem, the COIN CBC developers for making a great open-source ILP solver, and the google-ortools team for the mathematical modeling interface.
MOSS is dependent on 'ortools' and 'numpy'.
It is possible solve Sudoku with a two-index ILP formulation. However, the two-index solution is slower than commonly available heuristic solvers. It would be interesting to compare solve times using a commercial ILP solver such as Gurobi or CPLEX, versus the current open-source COIN-OR CBC solver.
As noted above, MOSS solves Sudoku by finding a feasible solution to the following system of equations: With the following constraints:
(b) SUM(x_ij)==45 across j for all i
(c) SUM(x_ij)==45 across i for all j
(d) SUM(x_ij)==45 for x_i_j within each nonet
(e) ABS(x_ij - x_i_notj)>=1 for all i, j
(f) ABS(x_ij - x_noti_j)>=1 for all i, j
Where x is the Sudoku board, and, x_ij are cell-values within the board.
While constraints [b, c, d] are relatively simple, constraints [e, f] are quite tricky to implement. Absolute value variables are non-linear, and therefore cannot be directly programmed into an ILP. The following approach is used to implement constraints [e, f]:
First, set variable 't' equal to the difference of each objective value pair. Set 't' bounds to [-10, 10]. For example:
(g) t_01_02 = x01 - x02
Second, initiate variables 'p' and 'n' with bounds [0, 10]. For example:
(h) 0 <= p_01_02 <= 10
(i) 0 <= n_01_02 <= 10
Third, initiate variable 'z' with bounds [1, 10]. For example:
(j) 1 <= z_01_02 <= 10
Fourth, initiate binary variable 'y'. For example:
(k) 0 <= y_01_02 <= 1
Fifth, set 't' equal to the difference of 'p' and 'n'. For example:
(l) t_01_02 = p_01_02 - n_01_02
Sixth, constraint that the difference of p and 10*y must be less than or equal to zero. For example:
(m) p_01_02 - 10*y_01-02 <= 0
Seventh, constraint that the sum of n and 10*y must be less than or equal to ten. For example:
(n) n_01_02 + 10*y_01_02 <= 10
Given steps 1-7, then 'z' will be the absolute value of the objective value pair. I.e.:
z_01_02 = abs(x01 - x02) iff constraints [g, h, i, j, k, l, m, n] are TRUE
Therefore, we apply [g - n] for each objective value pair and thereby constrain z to be greater than or equal to one. Thus, the all-different constraint is satisfied. (Assumes x_ij in [1, 9].)