# Minimal ellipses

This example comes from section 8.4.1 of the book Convex Optimization by Boyd and Vandenberghe (2004).

## Formulation

Given a set of $m$ ellipses of the form:

$$$E(A, b, c) = \{ x : x^\top A x + 2 b^\top x + c \leq 0 \},$$$

the minimal ellipse problem finds an ellipse with the minimum area that encloses the given ellipses.

It is convenient to parameterize the minimal enclosing ellipse as

$$$\{ x : || Px + q || \leq 1 \}.$$$

Then the optimal $P$ and $q$ are given by the convex semidefinite program;

\begin{aligned} \text{maximize } & \quad \log(\det(P)) \\ \text{subject to } & \quad \tau_i \geq 0, & i = 1, \ldots, m \\ & \quad\begin{bmatrix} P^2 - \tau_i A_i & P q - \tau_i b_i & 0 \\ (P q - \tau_i b_i)^\top & -1 - \tau_i c_i & (P q)^\top \\ 0 & (P q) & - P^2 \\ \end{bmatrix} \preceq 0 \text{ (PSD) } & i=1, \ldots, m \end{aligned}

with helper variables $\tau$.

## Required packages

This tutorial uses the following packages:

using JuMP
import LinearAlgebra
import Plots
import SCS

## Data

First, define the $m$ input ellipses (here $m = 6$), parameterized as $x^T A_i x + 2 b_i^T x + c \leq 0$:

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

We visualise the ellipses using the Plots package:

function plot_ellipse(plot, ellipse::Ellipse)
A, b, c = ellipse.A, ellipse.b, ellipse.c
θ = range(0, 2pi + 0.05; step = 0.05)
# Some linear algebra to convert θ into (x,y) coordinates.
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

## Build the model

Now let's build the model, using the change-of-variables P² = $P^2$ and P_q = $P q$. We'll recover the true value of P and q after the solve.

model = Model(SCS.Optimizer)
# We need to use a tighter tolerance for this example, otherwise the bounding
# ellipse won't actually be bounding...
set_optimizer_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 = [
#! format: off
(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²
#! format: on
]
@constraint(model, LinearAlgebra.Symmetric(X) <= 0, PSDCone())
end

We cannot directly represent the objective $\log(\det(P))$, so we introduce the conic reformulation:

@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$$$

Now, solve the program:

optimize!(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.04365e+00

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


## Results

After solving the model to optimality we can recover the solution in terms of $P$ and $q$:

P = sqrt(value.(P²))
q = P \ value.(P_q)
2-element Vector{Float64}:
-0.3964224834751066
-0.021382789710368295

Finally, overlaying the solution in the plot we see the minimal area enclosing ellipsoid:

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

Tip