function saida = normaH2(A, B, C, norma)
% function saida = normaH2(A, B, C, norma)
%
% Calcula a norma H-2 de um sistema linear a tempo contínuo por meio do problema de otimização:
%   min Tr(CPC')
%   s.a A'P+PA +BB' <= 0
%       P > 0
% 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 dinâmica do sistema
%            B -> matriz da entrada de ruído do sistema.
%            C -> matriz de saída do sistema.
%            norma -> limitante da norma a ser testado.
%
% saídas : saida.norma    -> valor da normal (infinito se instável)
%          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.W        -> variável soluçao W (se existir).
%
% Autor: ricfow@dt.fee.unicamp.br
% Data:  07/09/2010
%
% Exemplo:
% clear all;
% k1=1;
% k2=1;
% m1=1;
% m2=0.5;
% c0=2;
% A = [0 0 1 0;
% 0 0 0 1;
% -(k1+k2)/m1 k2/m1 -c0/m1 0;
% k2/m2 -k2/m2 0 -c0/m2];
% B = [0;0;1/m1;0];
% C = [0 1 0 0];
% saida = normaH2(A, B, C)

% Determina as dimensões do problema.
ordem    = size(A,1);
entradas = size(B,2);
saidas   = size(C,1);

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

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

% Cria as variáveis do problema.
W = sdpvar(ordem,ordem,'symmetric');
if nargin <= 3
    rho = sdpvar(1);
    obj = rho;
else
    rho = norma*norma;
    obj = [];
end

% Entra com as LMIs do problema.
LMIs = LMIs + set(W > 0) + set(A*W + W*A' + B*B' < 0);
LMIs = LMIs + set(rho >= trace(C*W*C'));

% 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.norma = inf;
% Capturando a solução (se existir).
if p > -1e-7
    saida.W     = double(W);    
    saida.norma = sqrt(double(rho));
end