Tips and tricks

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial was originally contributed by Arpit Bhatia.

This tutorial collates some tips and tricks you can use when formulating mixed-integer programs. It uses the following packages:

julia> using JuMP

Absolute value

To model the absolute value function $t \ge |x|$, there are a few options. In all cases, these reformulations only work if you are minimizing $t$ "down" into $|x|$. They do not work if you are trying to maximize $|x|$.

Option 1

This option adds two linear inequality constraints:

julia> model = Model();
julia> @variable(model, x)x
julia> @variable(model, t)t
julia> @constraint(model, t >= x)-x + t ≥ 0
julia> @constraint(model, t >= -x)x + t ≥ 0

Option 2

This option uses two non-negative variables and forms expressions for $x$ and $t$:

julia> model = Model();
julia> @variable(model, z[1:2] >= 0)2-element Vector{VariableRef}:
 z[1]
 z[2]
julia> @expression(model, t, z[1] + z[2])z[1] + z[2]
julia> @expression(model, x, z[1] - z[2])z[1] - z[2]

Option 3

This option uses MOI.NormOneCone and lets JuMP choose the reformulation:

julia> model = Model();
julia> @variable(model, x)x
julia> @variable(model, t)t
julia> @constraint(model, [t; x] in MOI.NormOneCone(2))[t, x] ∈ MathOptInterface.NormOneCone(2)

L1-norm

To model $\min ||x||_1$, that is, $\min \sum\limits_i |x_i|$, use the MOI.NormOneCone:

julia> model = Model();
julia> @variable(model, x[1:3])3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
julia> @variable(model, t)t
julia> @constraint(model, [t; x] in MOI.NormOneCone(1 + length(x)))[t, x[1], x[2], x[3]] ∈ MathOptInterface.NormOneCone(4)
julia> @objective(model, Min, t)t

Infinity-norm

To model $\min ||x||_\infty$, that is, $\min \max\limits_i |x_i|$, use the MOI.NormInfinityCone:

julia> model = Model();
julia> @variable(model, x[1:3])3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
julia> @variable(model, t)t
julia> @constraint(model, [t; x] in MOI.NormInfinityCone(1 + length(x)))[t, x[1], x[2], x[3]] ∈ MathOptInterface.NormInfinityCone(4)
julia> @objective(model, Min, t)t

Max

To model $t \ge \max\{x, y\}$, do:

julia> model = Model();
julia> @variable(model, t)t
julia> @variable(model, x)x
julia> @variable(model, y)y
julia> @constraint(model, t >= x)t - x ≥ 0
julia> @constraint(model, t >= y)t - y ≥ 0

This reformulation does not work for $t \ge \min\{x, y\}$.

Min

To model $t \le \min\{x, y\}$, do:

julia> model = Model();
julia> @variable(model, t)t
julia> @variable(model, x)x
julia> @variable(model, y)y
julia> @constraint(model, t <= x)t - x ≤ 0
julia> @constraint(model, t <= y)t - y ≤ 0

This reformulation does not work for $t \le \max\{x, y\}$.

Modulo

To model $y = x \text{ mod } n$, where $n$ is a constant modulus, we use the relationship $x = n \cdot z + y$, where $z \in \mathbb{Z}_+$ is the number of times that $n$ can be divided by $x$ and $y$ is the remainder.

julia> n = 44
julia> model = Model();
julia> @variable(model, x >= 0, Int)x
julia> @variable(model, 0 <= y <= n - 1, Int)y
julia> @variable(model, z >= 0, Int)z
julia> @constraint(model, x == n * z + y)x - y - 4 z = 0

The modulo reformulation is often useful for subdividing a time increment into units of time like hours and days:

julia> model = Model();
julia> @variable(model, t >= 0, Int)t
julia> @variable(model, 0 <= hours <= 23, Int)hours
julia> @variable(model, days >= 0, Int)days
julia> @constraint(model, t == 24 * days + hours)t - hours - 24 days = 0

Boolean operators

Binary variables can be used to construct logical operators. Here are some example.

Or

\[x_3 = x_1 \lor x_2\]

julia> model = Model();
julia> @variable(model, x[1:3], Bin)3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
julia> @constraints(model, begin
           x[1] <= x[3]
           x[2] <= x[3]
           x[3] <= x[1] + x[2]
       end)(x[1] - x[3] ≤ 0, x[2] - x[3] ≤ 0, -x[1] - x[2] + x[3] ≤ 0)

And

\[x_3 = x_1 \land x_2\]

julia> model = Model();
julia> @variable(model, x[1:3], Bin)3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
julia> @constraints(model, begin
           x[3] <= x[1]
           x[3] <= x[2]
           x[3] >= x[1] + x[2] - 1
       end)(-x[1] + x[3] ≤ 0, -x[2] + x[3] ≤ 0, -x[1] - x[2] + x[3] ≥ -1)

Not

\[x_1 \neg x_2\]

julia> model = Model();
julia> @variable(model, x[1:2], Bin)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @constraint(model, x[1] == 1 - x[2])x[1] + x[2] = 1

Implies

\[x_1 \implies x_2\]

julia> model = Model();
julia> @variable(model, x[1:2], Bin)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @constraint(model, x[1] <= x[2])x[1] - x[2] ≤ 0

Disjunctions

Problem

Suppose that we have two constraints $a^\top x \leq b$ and $c^\top x \leq d$, and we want at least one to hold.

Trick 1

Use an indicator constraint.

Example Either $x_1 \leq 1$ or $x_2 \leq 2$.

julia> model = Model();
julia> @variable(model, x[1:2])2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @variable(model, y[1:2], Bin)2-element Vector{VariableRef}:
 y[1]
 y[2]
julia> @constraint(model, y[1] --> {x[1] <= 1})y[1] --> {x[1] ≤ 1}
julia> @constraint(model, y[2] --> {x[2] <= 2})y[2] --> {x[2] ≤ 2}
julia> @constraint(model, sum(y) == 1)  # Exactly one branch must be truey[1] + y[2] = 1

Trick 2

Introduce a "big-M" multiplied by a binary variable to relax one of the constraints.

Example Either $x_1 \leq 1$ or $x_2 \leq 2$.

julia> model = Model();
julia> @variable(model, x[1:2] <= 10)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @variable(model, y[1:2], Bin)2-element Vector{VariableRef}:
 y[1]
 y[2]
julia> M = 100100
julia> @constraint(model, x[1] <= 1 + M * y[1])x[1] - 100 y[1] ≤ 1
julia> @constraint(model, x[2] <= 2 + M * y[2])x[2] - 100 y[2] ≤ 2
julia> @constraint(model, sum(y) == 1)y[1] + y[2] = 1

Indicator constraints

Problem

Suppose we want to model that a certain linear inequality must be satisfied when some other event occurs, that is, for a binary variable $z$, we want to model the implication:

\[z = 1 \implies a^\top x \leq b\]

Trick 1

Some solvers have native support for indicator constraints. In addition, if the variables involved have finite domains, then JuMP can automatically reformulate an indicator into a mixed-integer program.

Example $x_1 + x_2 \leq 1$ if $z = 1$.

julia> model = Model();
julia> @variable(model, 0 <= x[1:2] <= 10)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @variable(model, z, Bin)z
julia> @constraint(model, z --> {sum(x) <= 1})z --> {x[1] + x[2] ≤ 1}

Example $x_1 + x_2 \leq 1$ if $z = 0$.

julia> model = Model();
julia> @variable(model, 0 <= x[1:2] <= 10)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @variable(model, z, Bin)z
julia> @constraint(model, !z --> {sum(x) <= 1})!z --> {x[1] + x[2] ≤ 1}

Trick 2

If the solver doesn't support indicator constraints and the variables do not have a finite domain, you can use the big-M trick.

Example $x_1 + x_2 \leq 1$ if $z = 1$.

julia> model = Model();
julia> @variable(model, x[1:2] <= 10)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @variable(model, z, Bin)z
julia> M = 100100
julia> @constraint(model, sum(x) <= 1 + M * (1 - z))x[1] + x[2] + 100 z ≤ 101

Example $x_1 + x_2 \leq 1$ if $z = 0$.

julia> model = Model();
julia> @variable(model, x[1:2] <= 10)2-element Vector{VariableRef}:
 x[1]
 x[2]
julia> @variable(model, z, Bin)z
julia> M = 100100
julia> @constraint(model, sum(x) <= 1 + M * z)x[1] + x[2] - 100 z ≤ 1

Semi-continuous variables

A semi-continuous variable is a continuous variable between bounds $[l,u]$ that also can assume the value zero, that is: $x \in \{0\} \cup [l,u].$

Example $x \in \{0\}\cup [1, 2]$

julia> model = Model();
julia> @variable(model, x in Semicontinuous(1.0, 2.0))x

You can also represent a semi-continuous variable using the reformulation:

julia> model = Model();
julia> @variable(model, x)x
julia> @variable(model, z, Bin)z
julia> @constraint(model, x <= 2 * z)x - 2 z ≤ 0
julia> @constraint(model, x >= 1 * z)x - z ≥ 0

When z = 0 the two constraints are equivalent to 0 <= x <= 0. When z = 1, the two constraints are equivalent to 1 <= x <= 2.

Semi-integer variables

A semi-integer variable is a variable which assumes integer values between bounds $[l,u]$ and can also assume the value zero: $x \in \{0\} \cup [l, u] \cap \mathbb{Z}.$

julia> model = Model();
julia> @variable(model, x in Semiinteger(5.0, 10.0))x

You can also represent a semi-integer variable using the reformulation:

julia> model = Model();
julia> @variable(model, x, Int)x
julia> @variable(model, z, Bin)z
julia> @constraint(model, x <= 10 * z)x - 10 z ≤ 0
julia> @constraint(model, x >= 5 * z)x - 5 z ≥ 0

When z = 0 the two constraints are equivalent to 0 <= x <= 0. When z = 1, the two constraints are equivalent to 5 <= x <= 10.

Special Ordered Sets of Type 1

A Special Ordered Set of Type 1 is a set of variables, at most one of which can take a non-zero value, all others being at 0.

They most frequently apply where a set of variables are actually binary variables. In other words, we have to choose at most one from a set of possibilities.

julia> model = Model();
julia> @variable(model, x[1:3], Bin)3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
julia> @constraint(model, x in SOS1())[x[1], x[2], x[3]] ∈ MathOptInterface.SOS1{Float64}([1.0, 2.0, 3.0])

You can optionally pass SOS1 a weight vector like

julia> @constraint(model, x in SOS1([0.2, 0.5, 0.3]))[x[1], x[2], x[3]] ∈ MathOptInterface.SOS1{Float64}([0.2, 0.5, 0.3])

If the decision variables are related and have a physical ordering, then the weight vector, although not used directly in the constraint, can help the solver make a better decision in the solution process.

Special Ordered Sets of Type 2

A Special Ordered Set of type 2 is a set of non-negative variables, of which at most two can be non-zero, and if two are non-zero these must be consecutive in their ordering.

julia> model = Model();
julia> @variable(model, x[1:3])3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
julia> @constraint(model, x in SOS2([3.0, 1.0, 2.0]))[x[1], x[2], x[3]] ∈ MathOptInterface.SOS2{Float64}([3.0, 1.0, 2.0])

The ordering provided by the weight vector is more important in this case as the variables need to be consecutive according to the ordering. For example, in the above constraint, the possible pairs are:

• Consecutive
  • (x[1] and x[3]) as they correspond to 3 and 2 resp. and thus can be non-zero
  • (x[2] and x[3]) as they correspond to 1 and 2 resp. and thus can be non-zero
• Non-consecutive
  • (x[1] and x[2]) as they correspond to 3 and 1 resp. and thus cannot be non-zero

Piecewise linear approximations

SOSII constraints are most often used to form piecewise linear approximations of a function.

Given a set of points for x:

julia> x̂ = -1:0.5:2-1.0:0.5:2.0

and a set of corresponding points for y:

julia> ŷ = x̂ .^ 27-element Vector{Float64}:
 1.0
 0.25
 0.0
 0.25
 1.0
 2.25
 4.0

the piecewise linear approximation is constructed by representing x and y as convex combinations of x̂ and ŷ.

julia> N = length(x̂)7
julia> model = Model();
julia> @variable(model, -1 <= x <= 2)x
julia> @variable(model, y)y
julia> @variable(model, 0 <= λ[1:N] <= 1)7-element Vector{VariableRef}:
 λ[1]
 λ[2]
 λ[3]
 λ[4]
 λ[5]
 λ[6]
 λ[7]
julia> @objective(model, Max, y)y
julia> @constraints(model, begin
           x == sum(x̂[i] * λ[i] for i in 1:N)
           y == sum(ŷ[i] * λ[i] for i in 1:N)
           sum(λ) == 1
           λ in SOS2()
       end)(x + λ[1] + 0.5 λ[2] - 0.5 λ[4] - λ[5] - 1.5 λ[6] - 2 λ[7] = 0, y - λ[1] - 0.25 λ[2] - 0.25 λ[4] - λ[5] - 2.25 λ[6] - 4 λ[7] = 0, λ[1] + λ[2] + λ[3] + λ[4] + λ[5] + λ[6] + λ[7] = 1, [λ[1], λ[2], λ[3], λ[4], λ[5], λ[6], λ[7]] ∈ MathOptInterface.SOS2{Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]))