Советы
Это руководство было изначально предоставлено Арпитом Бхатией (Arpit Bhatia).
Это руководство призвано служить упрощенным введением в коническое программирование с использованием JuMP.
В нем используются следующие пакеты.
using JuMP
import SCS
import LinearAlgebra|
В этом руководстве используются множества из MathOptInterface. JuMP по умолчанию экспортирует символ |
import MathOptInterface as MOI|
Хорошим источником дополнительной информации о функциях, которые можно моделировать с помощью конусов, является рецептурный справочник по моделированию Mosek. |
Базовая теория
Подмножество векторного пространства является конусом, если и положительные скаляры , а произведение .
Конус является выпуклым, если для любого ] и любых .
Задачи конического программирования — это задачи выпуклой оптимизации, в которых выпуклая функция минимизируется на пересечении аффинного подпространства и выпуклого конуса. Вот пример задач минимизации конической формы в прямом виде:
Соответствующая двойственная задача:
где каждое является замкнутым выпуклым конусом, а — его двойственным конусом.
Конус второго порядка
SecondOrderCone (или конус Лоренца) размерности представляет собой конус следующего вида:
Чаще всего он применяется для представления L2-нормы вектора :
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[1:3])
@variable(model, t)
@constraint(model, sum(x) == 1)
@constraint(model, [t; x] in SecondOrderCone())
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), value.(x)(0.5773503148525505, [0.3333333333638094, 0.33333333336381227, 0.33333333336381227])Повернутый конус второго порядка
Конус второго порядка, повернутый на угол в плоскости , известен как RotatedSecondOrderCone. Он имеет следующую форму:
При u = 0.5 он представляет сумму квадратов вектора :
data = [1.0, 2.0, 3.0, 4.0]
target = [0.45, 1.04, 1.51, 1.97]
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, θ)
@variable(model, t)
@expression(model, residuals, θ * data .- target)
@constraint(model, [t; 0.5; residuals] in RotatedSecondOrderCone())
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
value(θ), value(t)(0.4980000000000055, 0.004979422159620655)Экспоненциальный конус
MOI.ExponentialCone — это множество следующего вида:
Оно может использоваться для моделирования задач с log и exp.
Экспонента
Для моделирования используйте (x, 1, z):
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x == 1.5)
@variable(model, z)
@objective(model, Min, z)
@constraint(model, [x, 1, z] in MOI.ExponentialCone())
optimize!(model)
@assert is_solved_and_feasible(model)
value(z), exp(1.5)(4.481654619603649, 4.4816890703380645)Логарифм
Для моделирования используйте (x, 1, z):
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x)
@variable(model, z == 1.5)
@objective(model, Max, x)
@constraint(model, [x, 1, z] in MOI.ExponentialCone())
optimize!(model)
@assert is_solved_and_feasible(model)
value(x), log(1.5)(0.4055956586655662, 0.4054651081081644)Логарифм суммы экспонент
Для моделирования используйте:
N = 3
x0 = rand(N)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[i = 1:N] == x0[i])
@variable(model, t)
@objective(model, Min, t)
@variable(model, u[1:N])
@constraint(model, sum(u) <= 1)
@constraint(model, [i = 1:N], [x[i] - t, 1, u[i]] in MOI.ExponentialCone())
optimize!(model)
value(t), log(sum(exp.(x0)))(1.4726849472138828, 1.472772274699489)Энтропия
Задача максимизации энтропии состоит в максимизации функции энтропии с учетом линейных ограничений неравенства.
Эту задачу можно смоделировать с помощью экспоненциального конуса, используя следующее преобразование:
Таким образом, задача принимает следующий вид:
m, n = 10, 15
A, b = randn(m, n), rand(m, 1)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t[1:n])
@variable(model, x[1:n])
@objective(model, Max, sum(t))
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b)
@constraint(model, [i = 1:n], [t[i], x[i], 1] in MOI.ExponentialCone())
optimize!(model)
@assert is_solved_and_feasible(model)
objective_value(model)2.663206172203935У MOI.ExponentialCone есть двойственная форма MOI.DualExponentialCone, которая может быть более эффективной для некоторых формулировок.
Есть также MOI.RelativeEntropyCone для явного программирования функции относительной энтропии.
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, x[1:n])
@objective(model, Max, -t)
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b)
@constraint(model, [t; ones(n); x] in MOI.RelativeEntropyCone(2n + 1))
optimize!(model)
@assert is_solved_and_feasible(model)
objective_value(model)2.6631736431749626PowerCone
MOI.PowerCone — это трехмерное множество, параметризованное скалярным значением α. Оно имеет следующую форму:
Степенной конус допускает ряд переформулировок. Например, при можно смоделировать , используя степенной конус с . Следовательно, смоделировать при можно так:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, x >= 1.5)
@constraint(model, [t, 1, x] in MOI.PowerCone(1 / 3))
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), value(x)(3.374754835774137, 1.49995522631532)У MOI.PowerCone есть двойственная форма MOI.DualPowerCone, которая может быть более эффективной для некоторых формулировок.
P-норма
P-норму можно смоделировать с помощью MOI.PowerCone. Вывод см. в рецептурном справочнике по моделированию Mosek.
function p_norm(x::Vector, p)
N = length(x)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, r[1:N])
@variable(model, t)
@constraint(model, [i = 1:N], [r[i], t, x[i]] in MOI.PowerCone(1 / p))
@constraint(model, sum(r) == t)
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
return value(t)
end
x = rand(5);
LinearAlgebra.norm(x, 4), p_norm(x, 4)(0.9316922467512209, 0.931687599495156)Положительно полуопределенный конус
Множество положительно полуопределенных матриц (PSD) размерности образует конус в . Это множество математически записывается так:
Конус PSD представляется в JuMP с использованием множеств MOI PositiveSemidefiniteConeTriangle (для верхнего треугольника матрицы PSD) и PositiveSemidefiniteConeSquare (для полной матрицы PSD). Однако предпочтительнее использовать сокращенную запись PSDCone, как показано ниже.
Пример: наибольшее собственное значение симметричной матрицы
Предположим, что у есть собственные значения . Тогда у матрицы есть собственные значения . Обратите внимание, что является положительно полуопределенной матрицей только тогда, когда все эти собственные значения неотрицательны, и это верно в случае со значениями . Таким образом, смоделировать задачу нахождения наибольшего собственного значения симметричной матрицы можно следующим образом:
A = [3 2 4; 2 0 2; 4 2 3]
I = Matrix{Float64}(LinearAlgebra.I, 3, 3)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@objective(model, Min, t)
@constraint(model, t .* I - A in PSDCone())
optimize!(model)
@assert is_solved_and_feasible(model)
objective_value(model)8.000003377698668GeometricMeanCone
MOI.GeometricMeanCone — это конус следующего вида:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[1:4])
@variable(model, t)
@constraint(model, sum(x) == 1)
@constraint(model, [t; x] in MOI.GeometricMeanCone(5))
optimize!(model)
value(t), value.(x)(0.0, [0.25000000110304854, 0.2500000011030519, 0.2500000011030504, 0.25000000110305026])RootDetCone
MOI.RootDetConeSquare — это конус следующего вида:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, X[1:2, 1:2])
@objective(model, Max, t)
@constraint(model, [t; vec(X)] in MOI.RootDetConeSquare(2))
@constraint(model, X .== [2 1; 1 3])
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), sqrt(LinearAlgebra.det(value.(X)))(2.236116364743881, 2.236067957402722)В случае симметричности X вместо этого можно использовать MOI.RootDetConeTriangle. Это может быть эффективнее, так как решателю не нужно добавлять дополнительные ограничения, чтобы гарантировать симметричность X.
При построении функции используйте triangle_vec для получения верхнего треугольника матрицы по столбцам в виде вектора в требуемом для JuMP порядке.
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, X[1:2, 1:2], Symmetric)
@objective(model, Max, t)
@constraint(model, [t; triangle_vec(X)] in MOI.RootDetConeTriangle(2))
@constraint(model, X .== [2 1; 1 3])
optimize!(model)
value(t), sqrt(LinearAlgebra.det(value.(X)))(2.2361792479821463, 2.2360679863317)LogDetCone
MOI.LogDetConeSquare — это конус следующего вида:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, u)
@variable(model, X[1:2, 1:2])
@objective(model, Max, t)
@constraint(model, [t; u; vec(X)] in MOI.LogDetConeSquare(2))
@constraint(model, X .== [2 1; 1 3])
@constraint(model, u == 0.5)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))(1.497920445822091, 1.497866100640713)В случае симметричности X вместо этого можно использовать MOI.LogDetConeTriangle. Это может быть эффективнее, так как решателю не нужно добавлять дополнительные ограничения, чтобы гарантировать симметричность X.
При построении функции используйте triangle_vec для получения верхнего треугольника матрицы по столбцам в виде вектора в требуемом для JuMP порядке.
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, u)
@variable(model, X[1:2, 1:2], Symmetric)
@objective(model, Max, t)
@constraint(model, [t; u; triangle_vec(X)] in MOI.LogDetConeTriangle(2))
@constraint(model, X .== [2 1; 1 3])
@constraint(model, u == 0.5)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))(1.4979204481049089, 1.497866101098124)Другие конусы и функции
Сведения о других конусах, поддерживаемых JuMP, см. в руководстве по MathOptInterface.