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.