[docs] add constraint program formulation of Sudoku

docs/Project.toml

@@ -5,6 +5,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+HiGHS_jll = "8fd58aa0-07eb-5a78-9b36-339c94fd15ea"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -23,14 +24,14 @@ CSV = "0.10"
 DataFrames = "1"
 Documenter = "0.27.9"
 GLPK = "=1.0.1"
-HTTP = "0.9"
+HTTP = "0.9, 1"
 HiGHS = "=1.1.3"
 Interpolations = "0.13"
 Ipopt = "=1.0.2"
 JSON = "0.21"
 JSONSchema = "1"
 Literate = "2.8"
-MathOptInterface = "=1.4.0"
+MathOptInterface = "=1.6.0"
 Plots = "1"
 SCS = "=1.1.2"
 StatsPlots = "0.14"

docs/src/tutorials/linear/sudoku.jl

@@ -38,7 +38,9 @@
 # a *feasibility* problem: we wish to find a feasible solution that satisfies
 # these rules. You can think of it as an optimization problem with an objective
 # of 0.
-#
+
+# ## Mixed-integer linear programming formulation
+
 # We can model this problem using 0-1 integer programming: a problem where all
 # the decision variables are binary. We'll use JuMP to create the model, and
 # then we can solve it with any integer programming solver.
@@ -154,3 +156,55 @@ sol
 # Which is the correct solution:
 
 # ![Solved Sudoku](../../assets/full_sudoku.png)
+
+# ## Constraint programming formulation
+
+# We can also model this problem using constraint programming and the
+# all-different constraint, which says that no two elements of a vector can take
+# the same value.
+
+# Because of the reformulation system in MathOptInterface, we can still solve
+# this problem using HiGHS.
+
+model = Model(HiGHS.Optimizer)
+set_silent(model)
+## HiGHS v1.2 has a bug in presolve which causes the problem to be classified as
+## infeasible.
+set_optimizer_attribute(model, "presolve", "off")
+
+# Instead of the binary variables, we directly define a 9x9 grid of integer
+# values between 1 and 9:
+
+@variable(model, 1 <= x[1:9, 1:9] <= 9, Int)
+
+# Then, we enforce that the values in each row must be all-different:
+
+@constraint(model, [i = 1:9], x[i, :] in MOI.AllDifferent(9))
+
+# That the values in each column must be all-different:
+
+@constraint(model, [j = 1:9], x[:, j] in MOI.AllDifferent(9))
+
+# And that the values in each 3x3 sub-grid must be all-different:
+
+for i in (0, 3, 6), j in (0, 3, 6)
+    @constraint(model, vec(x[i.+(1:3), j.+(1:3)]) in MOI.AllDifferent(9))
+end
+
+# Finally, as before we set the initial solution and optimize:
+
+for i in 1:9, j in 1:9
+    if init_sol[i, j] != 0
+        fix(x[i, j], init_sol[i, j]; force = true)
+    end
+end
+
+optimize!(model)
+
+# Display the solution
+
+csp_sol = round.(Int, value.(x))
+
+# Which is the same as we found before:
+
+sol == csp_sol