【CVX】Chebyshev center of a polyhedron

Chebyshev center of a polyhedron

We consider the problem of finding the largest Euclidean ball that lies in a polyhedron described by linear inequalities
P = { x ∈ R n ∣ a i T x ≤ b i , i = 1 , 2 , … , m } \mathcal{P}=\{x\in\mathbb{R}^n\mid a_i^Tx\leq b_i, i=1,2, \dots, m\} P={x∈Rn∣aiT​x≤bi​,i=1,2,…,m}
The ball is represented as follows:
B = { x c + u ∣ ∥ u ∥ 2 ≤ r } \mathcal{B}=\{x_c+u\mid\lVert u\rVert_2\leq r\} B={xc​+u∣∥u∥2​≤r}
And, we wish to maximize $r$ subject to the constraint $\mathcal{B}\subseteq \mathcal{P}$.
We start by considering the simpler constraint that $\mathcal{B}$ lies in one halfspace $a_i^{T}x\le b$, i.e.,
$$\Vert u{\Vert }_2\le r⟹a_i^T(x_c+u)\le b_i\text{\hspace{1em}(1)}$$
Since
sup ⁡ { a i T u ∣ ∥ u ∥ 2 ≤ r } = r ∥ a i ∥ 2 \sup\{a_i^Tu\mid \lVert u\rVert_2\leq r\}=r\lVert a_i \rVert_2 sup{aiT​u∣∥u∥2​≤r}=r∥ai​∥2​
we can write $(1)$ as
$$a_i^T{x}_c+r\Vert a_i{\Vert }_2\le b_i$$
which is a linear inequality in $x_c$ and $r$. Hence the Chebshev center can be determined by solving the LP
$$\begin{gathered} \max r \\ s.t.a_i^T{x}_c+r\Vert a_i{\Vert }_2\le b_i,i=1,2,\dots ,m \end{gathered}$$

CVX code

%% compute and display the Chebyshev center of a 2D polyhedron
% Generate the input data
a1 = [2; 1];
a2 = [2; -1];
a3 = [-1; 2];
a4 = [-1; -2];
b = ones(4, 1);

% solve model
cvx_begin
    variables r(1) x_c(2);
    maximize r;
    a1'*x_c+r*norm(a1, 2)<=b(1);
    a2'*x_c+r*norm(a2, 2)<=b(2);
    a3'*x_c+r*norm(a3, 2)<=b(3);
    a4'*x_c+r*norm(a4, 2)<=b(4);
cvx_end

% plot
x = linspace(-2, 2);
theta = 0:pi/100:2*pi;
hold on;
% constraints
plot(x, -x*a1(1)./a1(2)+b(1)./a1(2), 'b-');
plot(x, -x*a2(1)./a2(2)+b(2)./a2(2), 'b-');
plot(x, -x*a3(1)./a3(2)+b(3)./a3(2), 'b-');
plot(x, -x*a4(1)./a4(2)+b(4)./a4(2), 'b-');

% Euclidean ball
plot(x_c(1)+r*cos(theta), x_c(2)+r*sin(theta), 'r');

% chebshev center
plot(x_c(1), x_c(2), 'k+');

xlabel('x_1');
ylabel('x_2');
title('Largest Euclidean ball in polyhedron');
axis([-1 1 -1 1]);
axis equal
hold off;

Chebyshev center with random constraints

%%
rng(729);
n = 10; 
m = 2*n;
A = randn(m, n);
b = A*rand(n, 1)+2*rand(m, 1);
norm_ai = sum(A.^2, 2).^(0.5);

fprintf(1, 'Computing Chebyshev center...');
cvx_solver mosek;
cvx_begin
    variables r x_c(n);
    dual variable y
    maximize r;
    y: A*x_c+r*norm_ai<=b;
cvx_end

fprintf(1, 'The Chebyshev center coordinates are: \n');
disp(x_c);
fprintf(1, 'The radius of the largest Euclidean ball is: \n');
disp(r);

Reference

Convex Optimization. Boyd