function saida = iLyapm(A, Q, traco)
% function output = iLyap(A, Q)
%
% Resolve a inequação de Lyapunov por meio do seguinte problema de otimização:
%   min rho%   
%   s.a A'P+PA +Q < 0
%       P > 0
%       rho >= trace(P)
% As LMIs são programadas usando o YALMIP  e podem ser resolvidas
% por qualquer resolvedor de LMIs compatível com o YALMIP (SeDuM é o padrão).
%
% entradas : A -> matriz cuja estabilidade será testada.
%            Q -> matriz definida positiva.
%            traco -> valor da traço da matriz de Lyapunov.
%
% saídas : saida.factivel -> 1 se estável, 0 caso contrário.
%          saida.cpusec   -> tempo de cpu para resolver as LMIs (segundos)
%          saida.cpusec_m -> tempo de cpu para programar as LMIs (segundos)
%          saida.V        -> número de variáveis escalares.
%          saida.L        -> number de linhas de LMIs.
%          saida.P        -> variável soluçao P (se existir).
%          saida.traco    -> valor do traço da matriz P (em caso de factibilidade).
%
% Autor: ricfow@dt.fee.unicamp.br
% Data:  07/09/2010
%
% Exemplo:
%
% clear all;
% A = [0 1;-2 -3];
% Q = eye(2);
% traco = 1.5+1e-10;
% saida = iLyapm(A, Q, traco)

if nargin <= 2
    fprintf('Faltando parâmetro de entrada!\n');
    saida = [];
    return;
end
% Determina as dimensões do problema.
ordem = size(A,1);

% Tempo de programação.
saida.cpusec_m = clock;

% Inicializa a variável que armazenará as LMIs do problema.
LMIs = set([]);

% Cria as variáveis do problema.
P = sdpvar(ordem,ordem,'symmetric');
rho = sdpvar(1);

% Entra com as LMIs do problema.
LMIs = LMIs + set(P > 0) + set(-A'*P - P*A - Q > 0);

% Função objetivo.
LMIs = LMIs + set(traco >= trace(P));
obj = [];

% Finaliza o tempo de programação.
saida.cpusec_m = etime(clock,saida.cpusec_m);

% Determina o número de linhas de LMIs
saida.L = 0;
for i=1:size(LMIs,1)
    saida.L = saida.L + size(LMIs{i},1);
end

% Determina o número de variáveis escalares.
saida.V = size(getvariables(LMIs),2);

% Finalmente, o problema de otimização é resolvido.
solucao = solvesdp(LMIs,obj,sdpsettings('solver','sedumi'));
saida.cpusec = solucao.solvertime;

% Acessa o menor resíduo primal.
p=min(checkset(LMIs));
saida.delta = p;

saida.factivel = 0;
% Capturando a solução (se existir).
if p > 0
    saida.P     = double(P);    
    saida.traco = traco;
    saida.factivel = 1;
end