Example: Mixed integer optimization (MILP)
This example demonstrates the extension of ConstraintTree
bounds structures to accommodate 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
. These include ConstraintTrees.EqualTo
and ConstraintTrees.Between
, 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)
Main.var"Main".IntegerFromTo(1, 6)
...and include this bound in constraints and variables:
dice_system = C.variables(keys = [:first_dice, :second_dice], bounds = IntegerFromTo(1, 6))
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 2 elements:
:first_dice => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".Integ…
:second_dice => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".Integ…
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_constraint (generic function with 2 methods)
JuMP does not support direct integrality constraints, so we need to make a small digression 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))
C.traverse(cs) do c
isnothing(c.bound) || jump_constraint(model, x, c.value, c.bound)
end
JuMP.set_silent(model)
JuMP.optimize!(model)
JuMP.value.(model[:x])
end
milp_optimized_vars (generic function with 1 method)
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),
)
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 3 elements:
:first_dice => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".I…
:points_distance => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), ConstraintTrees.Between(1.5, Inf))
:second_dice => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".I…
For solving, we use GLPK (it has MILP capabilities).
import GLPK
dices_thrown = C.substitute_values(
dice_system,
milp_optimized_vars(
dice_system,
dice_system.first_dice.value + dice_system.second_dice.value,
GLPK.Optimizer,
),
)
ConstraintTrees.Tree{Float64} with 3 elements:
:first_dice => 6.0
:points_distance => 2.0
:second_dice => 4.0
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))
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 3 elements:
:a => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".IntegerFromTo(1…
:b => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".IntegerFromTo(1…
:c => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), Main.var"Main".IntegerFromTo(1…
For simpliclty, we make a shortcut for "values" in all variables:
v = C.map(C.value, vars, C.Value)
ConstraintTrees.Tree{ConstraintTrees.Value} with 3 elements:
:a => LinearValue([1], [1.0])
:b => LinearValue([2], [1.0])
:c => LinearValue([3], [1.0])
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)
ConstraintTrees.Tree{ConstraintTrees.Constraint} with 5 elements:
:a_less_than_b => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), ConstraintTrees.Between(0.0, Inf))
:b_less_than_c => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#), ConstraintTrees.Between(0.0, Inf))
:circumference => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#))
:right_angled => ConstraintTrees.Constraint(ConstraintTrees.QuadraticValue(#= ... =#), ConstraintTrees.EqualTo(0.0))
:sides => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 3 elements =#)
We will need a solver that supports both quadratic and integer optimization:
import SCIP
triangle_sides =
C.substitute_values(
triangle_system,
milp_optimized_vars(
triangle_system,
-triangle_system.circumference.value,
SCIP.Optimizer,
),
).sides
ConstraintTrees.Tree{Float64} with 3 elements:
:a => 3.0
:b => 4.0
:c => 5.0
This page was generated using Literate.jl.