Merge pull request #14 from COBREXA/sew-milp

Add ability to constraint variables to binary or integer values.

Project.toml

@@ -12,8 +12,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Clarabel = "0.6"
 DataStructures = "0.18"
 DocStringExtensions = "0.8, 0.9"
-SBML = "1"
+GLPK = "1"
 JuMP = "1"
+SBML = "1"
+SCIP = "0.11"
 julia = "1.6"
 
 [extras]
@@ -22,7 +24,8 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
+SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Clarabel", "Downloads", "GLPK", "JuMP", "SBML", "Test"]
+test = ["Clarabel", "Downloads", "GLPK", "JuMP", "SBML", "SCIP", "Test"]

docs/Project.toml

@@ -7,6 +7,7 @@ GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
+SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 
 [compat]
 Documenter = "1"

docs/src/quadratic-optimization.jl → docs/src/2-quadratic-optimization.jl

@@ -95,7 +95,7 @@ s *=
 # translates the constraints into JuMP `Model`s to support the quadratic
 # constraints.
 import JuMP
-function optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer)
+function quad_optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer)
     model = JuMP.Model(optimizer)
     JuMP.@variable(model, x[1:C.var_count(cs)])
     JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x))
@@ -121,7 +121,7 @@ end
 # We can now load a suitable optimizer and solve the system by maximizing the
 # negative squared error:
 import Clarabel
-st = C.constraint_values(s, optimized_vars(s, -s.objective.value, Clarabel.Optimizer))
+st = C.constraint_values(s, quad_optimized_vars(s, -s.objective.value, Clarabel.Optimizer))
 
 # If the optimization worked well, we can nicely get out the position of the
 # closest point to the line that is in the elliptical area:

docs/src/3-mixed-integer-optimization.jl

@@ -0,0 +1,133 @@
+
+# # Example: Mixed integer optimization (MILP)
+#
+# This example demonstrates the extension of `ConstraintTree` bounds structures
+# to accomodate new kinds of problems. In particular, we create a new kind of
+# `Bound` that is restricting the value to be a full integer, and then solve a
+# geometric problem with that.
+
+# ## Creating a custom bound
+#
+# All bounds contained in constraints are subtypes of the abstract
+# [`ConstraintTrees.Bound`](@ref). These include
+# [`ConstraintTrees.EqualTo`](@ref) and [`ConstraintTrees.Between`](@ref), but
+# the types can be extended as necessary, given the final rewriting of the
+# constraint system to JuMP can handle the new bounds.
+#
+# Let's make a small "marker" bound for something that needs to be integer-ish,
+# between 2 integers:
+
+import ConstraintTrees as C
+
+mutable struct IntegerFromTo <: C.Bound
+    from::Int
+    to::Int
+end
+
+# We can now write e.g. a bound on the number on a thrown six-sided die as
+# follows:
+
+IntegerFromTo(1, 6)
+
+# ...and include this bound in constraints and variables:
+
+dice_system =
+    C.variables(keys = [:first_dice, :second_dice], bounds = [IntegerFromTo(1, 6)])
+
+# Now the main thing that is left is to be able to translate this bound to JuMP
+# for solving. We can slightly generalize our constraint-translation system
+# from the previous examples for this purpose, by separating out the functions
+# that create the constraints:
+
+import JuMP
+
+function jump_constraint(m, x, v::C.Value, b::C.EqualTo)
+    JuMP.@constraint(m, C.substitute(v, x) == b.equal_to)
+end
+
+function jump_constraint(m, x, v::C.Value, b::C.Between)
+    isinf(b.lower) || JuMP.@constraint(m, C.substitute(v, x) >= b.lower)
+    isinf(b.upper) || JuMP.@constraint(m, C.substitute(v, x) <= b.upper)
+end
+
+# JuMP does not support direct integrality constraints, so we need to make a
+# small disgression with a slack variable:
+function jump_constraint(m, x, v::C.Value, b::IntegerFromTo)
+    var = JuMP.@variable(m, integer = true)
+    JuMP.@constraint(m, var >= b.from)
+    JuMP.@constraint(m, var <= b.to)
+    JuMP.@constraint(m, C.substitute(v, x) == var)
+end
+
+function milp_optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer)
+    model = JuMP.Model(optimizer)
+    JuMP.@variable(model, x[1:C.var_count(cs)])
+    JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x))
+    function add_constraint(c::C.Constraint)
+        isnothing(c.bound) || jump_constraint(model, x, c.value, c.bound)
+    end
+    function add_constraint(c::C.ConstraintTree)
+        add_constraint.(values(c))
+    end
+    add_constraint(cs)
+    JuMP.set_silent(model)
+    JuMP.optimize!(model)
+    JuMP.value.(model[:x])
+end
+
+# Let's try to solve a tiny system with the dice first. What's the best value
+# we can throw if the dice are thrown at least 1.5 points apart?
+
+dice_system *=
+    :points_distance^C.Constraint(
+        dice_system.first_dice.value - dice_system.second_dice.value,
+        C.Between(1.5, Inf),
+    )
+
+# For solving, we use GLPK (it has MILP capabilities).
+import GLPK
+dices_thrown = C.constraint_values(
+    dice_system,
+    milp_optimized_vars(
+        dice_system,
+        dice_system.first_dice.value + dice_system.second_dice.value,
+        GLPK.Optimizer,
+    ),
+)
+
+@test isapprox(dices_thrown.first_dice, 6.0) #src
+@test isapprox(dices_thrown.second_dice, 4.0) #src
+
+# ## A more realistic example with geometry
+#
+# Let's find the size of the smallest right-angled triangle with integer side
+# sizes (aka a Pythagorean triple).
+
+vars = C.variables(keys = [:a, :b, :c], bounds = (IntegerFromTo(1, 100),))
+
+# For simpliclty, we make a shortcut for "values" in all variables:
+v = C.tree_map(vars, C.value, C.Value)
+
+# With that shortcut, the constraint tree constructs quite easily:
+triangle_system =
+    :sides^vars *
+    :circumference^C.Constraint(sum(values(v))) *
+    :a_less_than_b^C.Constraint(v.b - v.a, (0, Inf)) *
+    :b_less_than_c^C.Constraint(v.c - v.b, (0, Inf)) *
+    :right_angled^C.Constraint(C.squared(v.a) + C.squared(v.b) - C.squared(v.c), 0.0)
+
+# We will need a solver that supports both quadratic and integer optimization:
+import SCIP
+triangle_sides =
+    C.constraint_values(
+        triangle_system,
+        milp_optimized_vars(
+            triangle_system,
+            -triangle_system.circumference.value,
+            SCIP.Optimizer,
+        ),
+    ).sides
+
+@test isapprox(triangle_sides.a, 3.0) #src
+@test isapprox(triangle_sides.b, 4.0) #src
+@test isapprox(triangle_sides.c, 5.0) #src

test/runtests.jl

@@ -4,11 +4,15 @@ using Test
 
 @testset "ConstraintTrees tests" begin
     @testset "Metabolic modeling" begin
-        include("../docs/src/metabolic-modeling.jl")
+        include("../docs/src/1-metabolic-modeling.jl")
     end
 
     @testset "Quadratic optimization" begin
-        include("../docs/src/quadratic-optimization.jl")
+        include("../docs/src/2-quadratic-optimization.jl")
+    end
+
+    @testset "Mixed-integer optimization" begin
+        include("../docs/src/3-mixed-integer-optimization.jl")
     end
 
     @testset "Miscellaneous methods" begin