Maximum volume inscribed ellipsoid in a polyhedron

% Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,...,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
%     maximize     log det B
%     subject to   || B a_i || + a_i^T d <= b,  for i = 1,...,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
  A(i,:) = null([px(i+1)-px(i) py(i+1)-py(i)])';
  b(i) = A(i,:)*.5*[px(i+1)+px(i); py(i+1)+py(i)];
  if A(i,:)*[pxint; pyint]-b(i)>0
    A(i,:) = -A(i,:);
    b(i) = -b(i);
  end
end

% formulate and solve the problem
cvx_begin
    variable B(n,n) symmetric
    variable d(n)
    maximize( det_rootn( B ) )
    subject to
       for i = 1:m
           norm( B*A(i,:)', 2 ) + A(i,:)*d <= b(i);
       end
cvx_end

% make the plots
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse_inner  = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer  = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), 'r--' );
plot( ellipse_outer(1,:), ellipse_outer(2,:), 'r--' );
axis square
axis off
hold off
 
Calling SDPT3: 28 variables, 9 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints =  9
 dim. of sdp    var  =  6,   num. of sdp  blk  =  2
 dim. of socp   var  = 15,   num. of socp blk  =  5
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      mean(obj)   cputime
-------------------------------------------------------------------
 0|0.000|0.000|9.1e+00|2.5e+00|1.6e+02| 7.205107e+00| 0:0:00| chol  1  1 
 1|1.000|0.871|2.8e-06|4.0e-01|2.4e+01| 5.831517e+00| 0:0:00| chol  1  1 
 2|1.000|1.000|4.6e-07|8.8e-03|2.7e+00| 1.671491e+00| 0:0:00| chol  1  1 
 3|0.830|1.000|8.5e-08|8.9e-04|4.8e-01| 8.822259e-01| 0:0:00| chol  1  1 
 4|1.000|0.901|2.8e-07|1.7e-04|5.4e-02| 9.487353e-01| 0:0:00| chol  1  1 
 5|0.967|0.974|1.0e-08|1.3e-05|1.6e-03| 9.522458e-01| 0:0:00| chol  1  1 
 6|0.969|0.976|6.5e-10|1.2e-06|4.4e-05| 9.523109e-01| 0:0:00| chol  1  1 
 7|0.941|0.978|8.0e-11|2.6e-08|2.1e-06| 9.523080e-01| 0:0:00| chol  1  1 
 8|1.000|1.000|5.4e-13|1.6e-11|2.3e-07| 9.523075e-01| 0:0:00| chol  1  1 
 9|0.994|1.000|7.5e-14|1.0e-12|6.7e-09| 9.523075e-01| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   =  9
 primal objective value =  9.52307512e-01
 dual   objective value =  9.52307506e-01
 gap := trace(XZ)       = 6.74e-09
 relative gap           = 2.32e-09
 actual relative gap    = 2.32e-09
 rel. primal infeas     = 7.55e-14
 rel. dual   infeas     = 1.00e-12
 norm(X), norm(y), norm(Z) = 1.9e+00, 2.6e+00, 4.6e+00
 norm(A), norm(b), norm(C) = 6.6e+00, 2.0e+00, 3.7e+00
 Total CPU time (secs)  = 0.2  
 CPU time per iteration = 0.0  
 termination code       =  0
 DIMACS: 7.5e-14  0.0e+00  1.1e-12  0.0e+00  2.3e-09  2.3e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.952308