%clc
echo off
%*********************************************************
%
% Decay-rate estimation using non-convex SDP
%
%*********************************************************
%
% The problem we solve is estimatation of decay-rate of a 
% linear system x' = Ax. This can be formulated as a 
% generalized eigenvalue problem (GEVP)
%
% max alpha
% s.t A'P+PA < -2alphaP
%          P > I
%
% This time, we solve it as a BMI using PENBMI (hence you 
% need PENBMI to run this demo). Note, this is a quasi-convex
% problem, and PENBMI is actually guaranteed to find the
% global optima.
%pause
%clc

% Define the variables
A = [-1 2;-3 -4];
P = sdpvar(2,2);
alpha = sdpvar(1,1);

%pause

F0 = set(P>eye(2))+set(A'*P+P*A < -eye(2));
options = sdpsettings('verbose',0,'warning',0);
solvesdp(F0,trace(P),options);
P0 = double(P); setsdpvar(P,P0);
alpha0 = -max(eig(inv(P0)*(A'*P0+P0*A)))/2; setsdpvar(alpha,alpha0);

% Define the GEVP
F = set(P>eye(2))+set(A'*P+P*A < -2*alpha*P);
%F = F + set(-100 < recover(depends(F)) < 100);

% Maximize alpha (minimize -alpha)
   options=sdpsettings('Solver','penbmi','verbose',1,...
      'penbmi.pbm_eps',1e-4,'penbmi.p_eps',1e-6,'penbmi.alpha',1e-2,...
      'penbmi.p0',.01,'penbmi.precision_2',1e-4,'usex0',0);
solvesdp(F,-alpha,options);
double(alpha)
%pause