Документация Engee

Советы

Это руководство было изначально предоставлено Арпитом Бхатией (Arpit Bhatia).

Хорошую подборку советов можно найти в рецептурном справочнике по моделированию Mosek.

В этом руководстве собраны некоторые советы и приемы, которые можно использовать при формулировании частично целочисленных программ. В нем используются следующие пакеты.

julia> using JuMP

Абсолютное значение

Есть несколько вариантов моделирования функции абсолютного значения . Все эти переформулировки работают только при минимизации «вниз» до . Они не работают при попытке максимизировать .

Вариант 1

В этом варианте добавляются два линейных ограничения неравенства:

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

Вариант 2

В этом варианте используются две неотрицательные переменные и формируются выражения для и :

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]

Вариант 3

В этом варианте используется MOI.NormOneCone, и JuMP выбирает переформулировку самостоятельно:

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-норма

Для моделирования , то есть , используйте 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

Норма бесконечности

Для моделирования , то есть , используйте 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

Максимум

Для моделирования сделайте следующее.

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

Эта переформулировка не работает для .

Минимум

Для моделирования сделайте следующее.

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

Эта переформулировка не работает для .

Деление по модулю

Для моделирования , где  — постоянный модуль, мы используем соотношение , где показывает, сколько раз можно разделить на , а  — остаток.

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

Переформулировка деления по модулю часто полезна для деления временного шага на единицы времени, такие как часы и дни:

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

Логические операторы

Двоичные переменные можно использовать для построения логических операторов. Вот ряд примеров.

Или

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)

И

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)

Не

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

Если…​, то…​

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

Дизъюнкции

Задача

Предположим, есть два ограничения и и должно соблюдаться хотя бы одно из них.

Прием 1

Используйте индикаторное ограничение.

Пример: либо , либо .

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)  # Одна и только одна ветвь должна иметь значение true
y[1] + y[2] = 1

Прием 2

Введите большое М (Big M), умноженное на двоичную переменную, чтобы ослабить одно из ограничений.

Пример: либо , либо .

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

Если M слишком мало, решение может быть неоптимальным. Если M слишком велико, решатель может столкнуться с численными проблемами. Попробуйте использовать знания предметной области, чтобы выбрать правильное значение M. В документации на сайте Gurobi есть хороший раздел по этой теме.

Индикаторные ограничения

Задача

Предположим, нужно смоделировать ситуацию, когда определенное линейное неравенство должно выполняться при наступлении какого-либо другого события, то есть для двоичной переменной нужно смоделировать импликацию:

Прием 1

В некоторых решателях есть встроенная поддержка индикаторных ограничений. Кроме того, если задействованные переменные имеют конечные области определения, JuMP может автоматически переформулировать индикатор в частично целочисленную программу.

Пример: , если .

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}

Пример: , если .

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}

Прием 2

Если решатель не поддерживает индикаторные ограничения, а переменные не имеют конечной области определения, можно использовать метод Big M.

Пример: , если .

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

Пример: , если .

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

Полунепрерывные переменные

Полунепрерывная переменная — это непрерывная переменная между границами ], которая также может принимать нулевое значение, то есть: .]

Пример: ]

julia> model = Model();


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

Полунепрерывную переменную можно также представить с помощью переформулировки:

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

При z = 0 эти два ограничения эквивалентны 0 <= x <= 0. При z = 1 эти два ограничения эквивалентны 1 <= x <= 2.

Частично целочисленные переменные

Частично целочисленная переменная — это переменная, которая принимает целые значения в пределах границ ], а также может принимать нулевое значение: \cap \mathbb{Z}.]

julia> model = Model();


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

Частично целочисленную переменную можно также представить с помощью переформулировки:

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

При z = 0 эти два ограничения эквивалентны 0 <= x <= 0. При z = 1 эти два ограничения эквивалентны 5 <= x <= 10.

Особые упорядоченные множества типа 1

Особое упорядоченное множество типа 1 — это множество переменных, максимум один элемент которого может принимать ненулевое значение, а остальные равны 0.

Чаще всего они применяются там, где множество переменных фактически содержит двоичные переменные. Иными словами, приходится выбирать максимум одну из множества возможностей.

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])

При необходимости можно передать для SOS1 весовой вектор, например:

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])

Если переменные решения связаны и упорядочены физически, весовой вектор хотя и не используется напрямую в ограничении, может помочь решателю принять более правильное решение.

Особые упорядоченные множества типа 2

Особое упорядоченное множество типа 2 — это множество неотрицательных переменных, из которых не более двух могут быть ненулевыми, причем они должны следовать по порядку.

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])

В этом случае более важен порядок, определяемый весовым вектором, так как переменные должны следовать по порядку. Например, в следующем ограничении возможны следующие пары:

  • Последовательные:

    • (x[1] и x[3]) соответствуют 3 и 2 и поэтому могут быть ненулевыми;

    • (x[2] и x[3]) соответствуют 1 и 2 и поэтому могут быть ненулевыми.

  • Непоследовательные:

    • (x[1] и x[2]) соответствуют 3 и 1 и поэтому не могут быть ненулевыми.

Кусочно-линейные аппроксимации

Ограничения SOSII чаще всего применяются для формирования кусочно-линейных аппроксимаций функции.

Для данного множества точек для x:

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

и множества соответствующих точек для y:

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

кусочно-линейная аппроксимация строится путем представления x и y как выпуклых комбинаций x̂ и ŷ.

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]))