% Boyd & Vandenberghe, "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% This is an example of D-optimal, A-optimal, and E-optimal
% experiment designs.

% problem data
m = 10;
angles1 = linspace(3*pi/4,pi,m);
angles2 = linspace(0,-pi/2,m);

% sensor positions
V = [3.0*[cos(angles1); sin(angles1)], ...
     1.5*[cos(angles2); sin(angles2)]];
p = size(V,2);
n = 2;
noangles = 5000;

% D-optimal design
%
%      maximize    log det V*diag(lambda)*V'
%      subject to  sum(lambda)=1,  lambda >=0
%

% setup the problem and solve it
cvx_begin
  variable lambda(p)
  maximize ( det_rootn( V*diag(lambda)*V' ) )
  subject to
    sum(lambda) == 1;
    lambda >= 0;
cvx_end
lambdaD = lambda; % save the solution for confidence ellipsoids

% plot results
figure(1)
% draw ellipsoid v'*W*v <= 2
W = inv(V*diag(lambda)*V');
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(2)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--', 0,0,'+');
set(d, 'Color', [0 0.5 0]); set(d(2),'MarkerFaceColor',[0 0.5 0]);
hold on;

dot=plot(V(1,:),V(2,:),'o');
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

% print out nonzero lambda
disp('Nonzero lambda values for D design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;

%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
hold off, axis off
% print -deps Ddesign.eps

% A-optimal design
%
%      minimize    Trace (sum_i lambdai*vi*vi')^{-1}
%      subject to  lambda >= 0, 1'*lambda = 1
%

% SDP formulation
e = eye(2,2);
cvx_begin sdp
  variables lambda(p) u(n)
  minimize ( sum(u) )
  subject to
    for k = 1:n
      [ V*diag(lambda)*V'  e(:,k);
        e(k,:)             u(k)   ] >= 0;
    end
    sum(lambda) == 1;
    lambda >= 0;
cvx_end
lambdaA = lambda; % save the solution for confidence ellipsoids

% plot results
figure(2)
% draw ellipsoid v'*W*v <= mu
W = inv(V*diag(lambda)*V')^2;
mu = diag(V'*W*V);
mu = mean(mu(ind));
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--',0,0,'+');
set(d, 'Color', [0 0.5 0]);
set(d(2), 'MarkerFaceColor', [0 0.5 0]);
hold on

dot = plot(V(1,:),V(2,:),'o');
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

disp('Nonzero lambda values for A design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;
%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off, hold off
% print -deps Adesign.eps

% E-optimal design
%
%      minimize    w
%      subject to  sum_i lambda_i*vi*vi' >= w*I
%                  lambda >= 0,  1'*lambda = 1;
%

cvx_begin sdp
  variables t lambda(p)
  maximize ( t )
  subject to
    V*diag(lambda)*V' >= t*eye(n,n);
    sum(lambda) == 1;
    lambda >= 0;
cvx_end

lambdaE = lambda; % save the solution for confidence ellipsoids

figure(3)
% draw ellipsoid v'*W*v <= mu
mu = diag(V'*W*V);
mu = mean(mu(ind));
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--', 0, 0, '+');
set(d, 'Color', [0 0.5 0]);
set(d(2), 'MarkerFaceColor', [0 0.5 0]);
hold on

dot = plot(V(1,:),V(2,:),'o');
lambda = lambda(1:p);
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

disp('Nonzero lambda values for E design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;
%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off, hold off
% print -deps Edesign.eps


% confidence ellipsoids
eta = 6.2514; % chi2inv(.9,3) value (command available in stat toolbox)
% draw 90 percent confidence ellipsoid  for D design
W = V*diag(lambdaD)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);

figure(4)
plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
text(ellipsoid(1,1100),ellipsoid(2,1100),'D');
hold on

% draw 90 percent confidence ellipsoid  for A design
W = V*diag(lambdaA)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
text(ellipsoid(1,1),ellipsoid(2,1),'A');

% draw 90 percent confidence ellipsoid  for E design
W = V*diag(lambdaE)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
d=plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
set(d,'Color',[0 0.5 0]);
text(ellipsoid(1,4000),ellipsoid(2,4000),'E');

% draw 90 percent confidence ellipsoid  for uniform design
W_u = inv(V*V'/p);
R = chol(W_u);  % W = R'*R
ellipsoid_u = sqrt(eta)*(R\[cos(angles); sin(angles)]);
plot(ellipsoid_u(1,:), ellipsoid_u(2,:), '--');
text(ellipsoid_u(1),ellipsoid_u(2),'U');
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off
% print -deps confidence.eps
hold off
 
Calling Mosek 9.1.9: 34 variables, 10 equality constraints
------------------------------------------------------------

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:32:15)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: MACOSX/64-X86

MOSEK warning 710: #2 (nearly) zero elements are specified in sparse col '' (12) of matrix 'A'.
MOSEK warning 710: #2 (nearly) zero elements are specified in sparse col '' (22) of matrix 'A'.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 10              
  Cones                  : 1               
  Scalar variables       : 24              
  Matrix variables       : 1               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 10              
  Cones                  : 1               
  Scalar variables       : 24              
  Matrix variables       : 1               
  Integer variables      : 0               

Optimizer  - threads                : 8               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 10
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 23                conic                  : 3               
Optimizer  - Semi-definite variables: 1                 scalarized             : 10              
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 49                after factor           : 49              
Factor     - dense dim.             : 0                 flops                  : 1.07e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   2.8e-01  2.8e-01  1.5e-01  9.00e-01   -6.099419583e-01  -6.151349439e-01  2.8e-01  0.01  
2   1.0e-01  1.0e-01  6.2e-02  3.60e-01   -2.100567633e+00  -1.935844370e+00  1.0e-01  0.01  
3   3.7e-02  3.7e-02  1.3e-02  1.00e+00   -2.477214979e+00  -2.423814987e+00  3.7e-02  0.01  
4   1.2e-02  1.2e-02  2.5e-03  8.96e-01   -2.871811315e+00  -2.851360554e+00  1.2e-02  0.01  
5   7.4e-04  7.4e-04  4.0e-05  8.85e-01   -3.161111814e+00  -3.159716259e+00  7.4e-04  0.01  
6   2.1e-05  2.1e-05  2.0e-07  9.94e-01   -3.181508258e+00  -3.181464753e+00  2.1e-05  0.01  
7   3.4e-06  3.4e-06  1.3e-08  9.98e-01   -3.181896330e+00  -3.181889011e+00  3.4e-06  0.01  
8   3.8e-07  3.8e-07  4.9e-10  9.99e-01   -3.181969805e+00  -3.181968986e+00  3.8e-07  0.01  
9   3.9e-08  3.9e-08  1.6e-11  1.00e+00   -3.181979263e+00  -3.181979179e+00  3.9e-08  0.01  
10  3.5e-09  3.6e-09  4.2e-13  1.00e+00   -3.181980399e+00  -3.181980392e+00  3.5e-09  0.01  
Optimizer terminated. Time: 0.02    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -3.1819803992e+00   nrm: 7e+00    Viol.  con: 7e-09    var: 2e-09    barvar: 0e+00    cones: 0e+00  
  Dual.    obj: -3.1819803917e+00   nrm: 3e+00    Viol.  con: 0e+00    var: 1e-15    barvar: 8e-09    cones: 0e+00  
Optimizer summary
  Optimizer                 -                        time: 0.02    
    Interior-point          - iterations : 10        time: 0.01    
      Basis identification  -                        time: 0.00    
        Primal              - iterations : 0         time: 0.00    
        Dual                - iterations : 0         time: 0.00    
        Clean primal        - iterations : 0         time: 0.00    
        Clean dual          - iterations : 0         time: 0.00    
    Simplex                 -                        time: 0.00    
      Primal simplex        - iterations : 0         time: 0.00    
      Dual simplex          - iterations : 0         time: 0.00    
    Mixed integer           - relaxations: 0         time: 0.00    

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +3.18198
 
Nonzero lambda values for D design:
lambda(1) = 0.50003
lambda(10) = 0.49997
 
Calling Mosek 9.1.9: 32 variables, 11 equality constraints
------------------------------------------------------------

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:32:15)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: MACOSX/64-X86

MOSEK warning 710: #4 (nearly) zero elements are specified in sparse col '' (9) of matrix 'A'.
MOSEK warning 710: #4 (nearly) zero elements are specified in sparse col '' (19) of matrix 'A'.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 11              
  Cones                  : 0               
  Scalar variables       : 20              
  Matrix variables       : 2               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 11              
  Cones                  : 0               
  Scalar variables       : 20              
  Matrix variables       : 2               
  Integer variables      : 0               

Optimizer  - threads                : 8               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 11
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 20                conic                  : 0               
Optimizer  - Semi-definite variables: 2                 scalarized             : 12              
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 46                after factor           : 46              
Factor     - dense dim.             : 0                 flops                  : 1.54e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   2.0e+00  1.0e+00  3.0e+00  0.00e+00   2.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   3.6e-01  1.8e-01  3.0e-01  1.43e-01   1.338526336e+00   9.160379630e-01   1.8e-01  0.01  
2   1.1e-01  5.6e-02  4.3e-02  1.11e+00   1.046442677e+00   9.128418787e-01   5.6e-02  0.01  
3   4.2e-02  2.1e-02  8.9e-03  1.47e+00   9.808576717e-01   9.399794279e-01   2.1e-02  0.01  
4   9.1e-03  4.6e-03  9.6e-04  1.10e+00   8.851796681e-01   8.769996235e-01   4.6e-03  0.01  
5   6.7e-04  3.3e-04  2.1e-05  9.61e-01   8.520571236e-01   8.514922333e-01   3.3e-04  0.01  
6   1.6e-05  7.8e-06  7.5e-08  9.94e-01   8.495899872e-01   8.495766969e-01   7.8e-06  0.01  
7   6.8e-07  3.4e-07  6.9e-10  9.99e-01   8.495305055e-01   8.495299248e-01   3.4e-07  0.01  
8   7.8e-08  3.9e-08  2.6e-11  1.00e+00   8.495281947e-01   8.495281283e-01   3.9e-08  0.01  
9   9.4e-09  4.7e-09  1.1e-12  1.00e+00   8.495279546e-01   8.495279466e-01   4.7e-09  0.01  
Optimizer terminated. Time: 0.01    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 8.4952795464e-01    nrm: 5e+00    Viol.  con: 5e-09    var: 0e+00    barvar: 0e+00  
  Dual.    obj: 8.4952794664e-01    nrm: 1e+00    Viol.  con: 0e+00    var: 1e-16    barvar: 5e-09  
Optimizer summary
  Optimizer                 -                        time: 0.01    
    Interior-point          - iterations : 9         time: 0.01    
      Basis identification  -                        time: 0.00    
        Primal              - iterations : 0         time: 0.00    
        Dual                - iterations : 0         time: 0.00    
        Clean primal        - iterations : 0         time: 0.00    
        Clean dual          - iterations : 0         time: 0.00    
    Simplex                 -                        time: 0.00    
      Primal simplex        - iterations : 0         time: 0.00    
      Dual simplex          - iterations : 0         time: 0.00    
    Mixed integer           - relaxations: 0         time: 0.00    

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.849528
 
Nonzero lambda values for A design:
lambda(1) = 0.2966
lambda(10) = 0.37795
lambda(20) = 0.32544
 
Calling Mosek 9.1.9: 23 variables, 3 equality constraints
------------------------------------------------------------

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:32:15)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: MACOSX/64-X86

MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (3) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (12) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (22) of matrix 'A'.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 3               
  Cones                  : 1               
  Scalar variables       : 23              
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 3               
  Cones                  : 1               
  Scalar variables       : 23              
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer  - threads                : 8               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 3
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 23                conic                  : 3               
Optimizer  - Semi-definite variables: 0                 scalarized             : 0               
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 6                 after factor           : 6               
Factor     - dense dim.             : 0                 flops                  : 2.44e+02        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.0e+00  9.0e+00  2.0e+00  0.00e+00   1.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   2.1e-01  1.9e+00  3.4e-01  -2.77e-01  -1.792668816e+00  -1.817470354e+00  2.1e-01  0.01  
2   4.0e-02  3.6e-01  1.2e-02  1.42e+00   -1.687864078e+00  -1.759043337e+00  4.0e-02  0.01  
3   7.3e-03  6.6e-02  1.1e-03  1.17e+00   -1.797933288e+00  -1.808227493e+00  7.3e-03  0.01  
4   1.1e-04  9.6e-04  2.0e-06  1.05e+00   -1.799891394e+00  -1.800034024e+00  1.1e-04  0.01  
5   2.1e-07  1.9e-06  1.7e-10  1.00e+00   -1.799999788e+00  -1.800000066e+00  2.1e-07  0.01  
6   2.9e-13  2.6e-12  2.4e-16  1.00e+00   -1.800000000e+00  -1.800000000e+00  2.9e-13  0.01  
Optimizer terminated. Time: 0.01    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -1.8000000000e+00   nrm: 1e+00    Viol.  con: 3e-13    var: 1e-13    cones: 0e+00  
  Dual.    obj: -1.8000000000e+00   nrm: 2e+00    Viol.  con: 0e+00    var: 3e-13    cones: 0e+00  
Optimizer summary
  Optimizer                 -                        time: 0.01    
    Interior-point          - iterations : 6         time: 0.01    
      Basis identification  -                        time: 0.00    
        Primal              - iterations : 0         time: 0.00    
        Dual                - iterations : 0         time: 0.00    
        Clean primal        - iterations : 0         time: 0.00    
        Clean dual          - iterations : 0         time: 0.00    
    Simplex                 -                        time: 0.00    
      Primal simplex        - iterations : 0         time: 0.00    
      Dual simplex          - iterations : 0         time: 0.00    
    Mixed integer           - relaxations: 0         time: 0.00    

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.8
 
Nonzero lambda values for E design:
lambda(10) = 0.2
lambda(20) = 0.8