Минимальные эллипсы

Этот пример взят из раздела 8.4.1 книги Convex Optimization авторов Boyd and Vandenberghe (2004).

Формулировка

Если дано множество из эллипсов вида:

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

Минимальный охватывающий эллипс удобно параметризовать следующим образом:

Тогда оптимальные и определяются выпуклой полуопределенной программой:

со вспомогательными переменными .

Требуемые пакеты

В этом руководстве используются следующие пакеты.

using JuMP
import LinearAlgebra
import Plots
import SCS
import Test

Данные

Сначала определим входных эллипсов (в данном случае ), параметризованных как :

struct Ellipse
    A::Matrix{Float64}
    b::Vector{Float64}
    c::Float64
    function Ellipse(A::Matrix{Float64}, b::Vector{Float64}, c::Float64)
        @assert isreal(A) && LinearAlgebra.issymmetric(A)
        return new(A, b, c)
    end
end

ellipses = [
    Ellipse([1.2576 -0.3873; -0.3873 0.3467], [0.2722, 0.1969], 0.1831),
    Ellipse([1.4125 -2.1777; -2.1777 6.7775], [-1.228, -0.0521], 0.3295),
    Ellipse([1.7018 0.8141; 0.8141 1.7538], [-0.4049, 1.5713], 0.2077),
    Ellipse([0.9742 -0.7202; -0.7202 1.5444], [0.0265, 0.5623], 0.2362),
    Ellipse([0.6798 -0.1424; -0.1424 0.6871], [-0.4301, -1.0157], 0.3284),
    Ellipse([0.1796 -0.1423; -0.1423 2.6181], [-0.3286, 0.557], 0.4931),
];

Визуализируем эллипсы с помощью пакета Plots:

function plot_ellipse(plot, ellipse::Ellipse)
    A, b, c = ellipse.A, ellipse.b, ellipse.c
    θ = range(0, 2pi + 0.05; step = 0.05)
    # Немного линейной алгебры для преобразования θ в координаты (x,y).
    x_y = √A \ (√(b' * (A \ b) - c) .* hcat(cos.(θ), sin.(θ)) .- (√A \ b)')'
    Plots.plot!(plot, x_y[1, :], x_y[2, :]; label = nothing, c = :navy)
    return
end

plot = Plots.plot(; size = (600, 600))
for ellipse in ellipses
    plot_ellipse(plot, ellipse)
end
plot

Построение модели

Теперь построим модель, используя замену переменных P² = и P_q = . После решения мы восстановим истинные значения P и q.

model = Model(SCS.Optimizer)
# В этом примере нужно использовать более строгий допуск, иначе ограничивающий
# эллипс на самом деле не будет ограничивающим...
set_attribute(model, "eps_rel", 1e-6)
set_silent(model)
m, n = length(ellipses), size(first(ellipses).A, 1)
@variable(model, τ[1:m] >= 0)
@variable(model, P²[1:n, 1:n], PSD)
@variable(model, P_q[1:n])

for (i, ellipse) in enumerate(ellipses)
    A, b, c = ellipse.A, ellipse.b, ellipse.c
    X = [
        #! формат: выкл.
        (P² - τ[i] * A)   (P_q - τ[i] * b) zeros(n, n)
        (P_q - τ[i] * b)' (-1 - τ[i] * c)  P_q'
        zeros(n, n)       P_q              -P²
        #! формат: вкл.
    ]
    @constraint(model, LinearAlgebra.Symmetric(X) <= 0, PSDCone())
end

Целевую функцию представить напрямую невозможно, поэтому введем коническую переформулировку:

@variable(model, log_det_P)
@constraint(model, [log_det_P; 1; vec(P²)] in MOI.LogDetConeSquare(n))
@objective(model, Max, log_det_P)

log_det_P

Теперь решим программу:

optimize!(model)
Test.@test is_solved_and_feasible(model)
solution_summary(model)

* Solver : SCS

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "solved"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : -4.04358e+00
  Dual objective value : -4.04364e+00

* Work counters
  Solve time (sec)   : 2.70487e-01

Результаты

После решения модели можно восстановить решение в терминах и :

P = sqrt(value.(P²))
q = P \ value.(P_q)

2-element Vector{Float64}:
 -0.396421769654888
 -0.02139417906487509

Наконец, наложив решение на график, мы увидим охватывающий эллипсоид с минимальной площадью:

Plots.plot!(
    plot,
    [tuple(P \ [cos(θ) - q[1], sin(θ) - q[2]]...) for θ in 0:0.05:(2pi+0.05)];
    c = :crimson,
    label = nothing,
)