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.

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 , there are a few options. In all cases, these reformulations only work if you are minimizing "down" into . They do not work if you are trying to maximize .

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 and :

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 , that is, , 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 , that is, , 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 , 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 .

Min

To model , 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 .

Modulo

To model , where is a constant modulus, we use the relationship , where is the number of times that can be divided by and is the remainder.

julia> n = 4
4

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

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

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

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

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 and , and we want at least one to hold.

Trick 1

Use an indicator constraint.

Example Either or .

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 true
y[1] + y[2] = 1

Trick 2

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

Example Either or .

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 = 100
100

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

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, that is, for a binary variable , we want to model the implication:

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

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

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

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 = 100
100

julia> @constraint(model, sum(x) <= 1 + M * (1 - z))
x[1] + x[2] + 100 z ≤ 101

Example if .

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 = 100
100

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 ] that also can assume the value zero, that is: .]

Example ]

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 ] and can also assume the value zero: \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> ŷ = x̂ .^ 2
7-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 ŷ.

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(ŷ[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]))