Tips and tricks

Originally Contributed by: Arpit Bhatia

Tip

A good source of tips is the Mosek Modeling Cookbook.

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.0
julia> @constraint(model, t >= -x)x + t ≥ 0.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.0
julia> @constraint(model, t >= y)t - y ≥ 0.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.0
julia> @constraint(model, t <= y)t - y ≤ 0.0

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

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.0, x[2] - x[3] ≤ 0.0, -x[1] - x[2] + x[3] ≤ 0.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.0, -x[2] + x[3] ≤ 0.0, -x[1] - x[2] + x[3] ≥ -1.0)

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.0

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.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

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])2-element Vector{VariableRef}: x[1] x[2]
julia> @variable(model, y, Bin)y
julia> M = 100100
julia> @constraint(model, x[1] <= 1 + M * y)x[1] - 100 y ≤ 1.0
julia> @constraint(model, x[2] <= 2 + M * (1 - y))x[2] + 100 y ≤ 102.0
Warning

If M is too small, the solution may be suboptimal. If M is too big, the solver may encounter numerical issues. Try to use domain knowledge to choose an M that is just right. Gurobi has a good documentation section on this topic.

Indicator constraints

Problem

Suppose we want to model that a certain linear inequality must be satisfied when some other event occurs, i.e., 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.

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

julia> model = Model();
julia> @variable(model, x[1:2])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.0}

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

julia> model = Model();
julia> @variable(model, x[1:2])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.0}

Trick 2

If the solver doesn't support indicator constraints, you an use the big-M trick.

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

julia> model = Model();
julia> @variable(model, x[1:2])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.0

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

julia> model = Model();
julia> @variable(model, x[1:2])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.0

Semi-continuous variables

Info

This section uses sets from MathOptInterface. By default, JuMP exports the MOI symbol as an alias for the MathOptInterface.jl package. We recommend making this more explicit in your code by adding the following lines:

import MathOptInterface
const MOI = MathOptInterface

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 MOI.Semicontinuous(1.0, 2.0))x

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 MOI.Semiinteger(5.0, 10.0))x

Special Ordered Sets of Type I

A Special Ordered Set of Type I 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 II

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 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.0, y - λ[1] - 0.25 λ[2] - 0.25 λ[4] - λ[5] - 2.25 λ[6] - 4 λ[7] = 0.0, λ[1] + λ[2] + λ[3] + λ[4] + λ[5] + λ[6] + λ[7] = 1.0, [λ[1], λ[2], λ[3], λ[4], λ[5], λ[6], λ[7]] ∈ MathOptInterface.SOS2{Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]))

Tip

This tutorial was generated using Literate.jl. View the source .jl file on GitHub.