function Info = h2_OdOP08_SCL_T11_star_yal(Ai,Bi,Ci,kappa,param)
% function Info = h2_OdOP08_SCL_T11_star_yal(Ai,Bi,Ci,kappa,param)
% Evaluate a H-2 guaranteed cost for a continuous-time time-invariant polytopic linear 
% system the notion of kappa-robust H-2 performance presented in Theorem 11 of (Oliveira, 
% de Oliveira & Peres, SCL 2008) with P_k(alpha) = P_k (independent of alpha). The LMIs 
% are programmed using YALMIP  and can be solved  by any LMI solver supported by 
% YALMIP (SeDuMi is the default).
%
% inputs:  (Ai,Bi,Ci) -> vertices of the system matrices (Ai = [A1 .. AN], Bi = [B1 .. BN], ..)
%           kappa     -> value of kappa  (k >= 1)%
%           param   -> optional parameters
%           param.with_P_def_pos  -> add the constraint P_k(alpha) > 0. If at least one point inside
%                                 the polytope (for instance, one of the vertices) is known 
%                                 to be stable, this constraint can be dropped.
%           param.h2              -> if this parameter is different of zero, the routine will
%                                 check the feasibility of the LMIs for this value of the h2
%                                 norm.
%output: Info.h2       -> H-2 guaranteed cost (0 if unfeasible)
%        Info.cpusec   -> cpu time to solve the LMIs (seconds)
%        Info.cpusec_m -> cpu time to mount the LMIs (seconds)
%        Info.Pk       -> Solution variable Pk
%        Info.K        -> number of scalar variables used in the optmization problem
%        Info.L        -> number of LMI rows used in the optmization problem
%        Info.delta    -> minimal primal residual returned by the LMI solver (SeDuMi is the default).
%
%example: 2 states and 2 vertices
% A1 = [-0.9  0.2;
%       -0.5  -1.9 ];
% A2 = [-1.1  -0.1;
%       0.5    -2.8];
% Ai = [A1 A2];
% Bi = [1 2;0 0]; Ci = [1 0 1 0];
% Info = h2_OdOP08_SCL_T11_star_yal(Ai,Bi,Ci,2)
%
% Date: 20/06/2007
% Author: ricfow@dt.fee.unicamp.br

if nargin == 5
    if isfield(param,'with_P_def_pos')
        with_P_def_pos = param.with_P_def_pos;
    else
        with_P_def_pos = 0;
    end
    if isfield(param,'h2')
        h2 = param.h2;
    else
        h2 = 0;
    end
else
    with_P_def_pos = 0;
    h2 = 0;
end

[order,vertices] = size(Ai);
vertices = vertices / order;
[order,inputs] = size(Bi);
inputs = inputs / vertices;
outputs = size(Ci,1);

Info.cpusec_m = clock;
LMIs = set([]);

Info.L = 0;

L_k = kron([eye(kappa) zeros(kappa,1)]',eye(order));
R_k = kron([zeros(kappa,1) eye(kappa)],eye(order));

if h2 == 0
    gamma = sdpvar(1,1);
    obj = gamma;
else
    gamma = h2^2;
    obj = [];
end

Pk = sdpvar(kappa*order,kappa*order,'symmetric');
if kappa == 1 & with_P_def_pos == 1
    LMIs = LMIs + set(Pk > 0);
    Info.L = Info.L + kappa*order;
end
Q2 = blkdiag(eye(outputs),L_k*Pk*R_k + R_k'*Pk*L_k');
%creating the variables
for i=1:vertices
    
    W{i} = sdpvar(inputs,inputs,'symmetric');
    Q1{i} = blkdiag(-W{i},Pk);        
    
    if kappa > 1 & with_P_def_pos == 1
        Y{i} = sdpvar(kappa*order,(kappa-1)*order,'full');
    end
    V{i} = sdpvar(inputs+kappa*order,kappa*order,'full');
    Z{i} = sdpvar(outputs+(kappa+1)*order,outputs+kappa*order,'full');
end

%programming the LMIs
for i=1:vertices
    A = Ai(:,order*(i-1)+1:i*order);
    B = Bi(:,inputs*(i-1)+1:i*inputs);
    C = Ci(:,order*(i-1)+1:i*order);

    %P_k > YC_k + C_k'Y'
    Ck_minus_one_i = kron([eye(kappa-1) zeros(kappa-1,1)],A) - kron([zeros(kappa-1,1) eye(kappa-1)],eye(order));
    if kappa > 1 & with_P_def_pos == 1
        LMIs = LMIs + set(Pk > Y{i}*Ck_minus_one_i + Ck_minus_one_i'*Y{i}');
        Info.L = Info.L + kappa*order;
    end    

    %trace condition    
    Fk_i = [B -eye(order) zeros(order,(kappa-1)*order);
             zeros((kappa-1)*order,inputs) Ck_minus_one_i];
    LMIs = LMIs + set(Q1{i} + V{i}*Fk_i + Fk_i'*V{i}' < 0);
    Info.L = Info.L + inputs + kappa*order;

    %grammiam condition    
    Gk_i = [eye(outputs) -C zeros(outputs,kappa*order);
             zeros(kappa*order,outputs) kron([eye(kappa) zeros(kappa,1)],A) - kron([zeros(kappa,1) eye(kappa)],eye(order))];
    LMIs = LMIs + set(Q2 + Z{i}*Gk_i + Gk_i'*Z{i}' < 0);
    Info.L = Info.L + outputs+(kappa+1)*order;
    
    LMIs = LMIs + set(trace(W{i}) < gamma);
    Info.L = Info.L + 1;

    for j=i+1:vertices
        Aj = Ai(:,order*(j-1)+1:j*order);
        Bj = Bi(:,inputs*(j-1)+1:j*inputs);
        Cj = Ci(:,order*(j-1)+1:j*order);

        Ck_minus_one_j = kron([eye(kappa-1) zeros(kappa-1,1)],Aj) - kron([zeros(kappa-1,1) eye(kappa-1)],eye(order));
        if kappa > 1 & with_P_def_pos == 1            
            LMIs = LMIs + set(2*Pk  > Y{i}*Ck_minus_one_j + Ck_minus_one_j'*Y{i}' + Y{j}*Ck_minus_one_i + Ck_minus_one_i'*Y{j}');
            Info.L = Info.L + kappa*order;
        end        

        %trace condition        
        Fk_j = [Bj -eye(order) zeros(order,(kappa-1)*order);
             zeros((kappa-1)*order,inputs) Ck_minus_one_j];
        
        LMIs = LMIs + set(Q1{i} + Q1{j} + V{i}*Fk_j + Fk_j'*V{i}' + V{j}*Fk_i + Fk_i'*V{j}' < 0);
        Info.L = Info.L + inputs + kappa *order;

        %grammiam condition
        Gk_j = [eye(outputs) -Cj zeros(outputs,kappa*order);
             zeros(kappa*order,outputs) kron([eye(kappa) zeros(kappa,1)],Aj) - kron([zeros(kappa,1) eye(kappa)],eye(order))];
        
        LMIs = LMIs + set(2*Q2 + Z{i}*Gk_j + Gk_j'*Z{i}' + Z{j}*Gk_i + Gk_i'*Z{j}' < 0);
        Info.L = Info.L + outputs+(kappa+1)*order;
    end
end

Info.cpusec_m = etime(clock,Info.cpusec_m);

Info.K = size(getvariables(LMIs),2);

sol = solvesdp(LMIs,obj,sdpsettings('verbose',0,'solver','sedumi'));
Info.cpusec = sol.solvertime;
p=min(checkset(LMIs));
Info.delta = p;
if p > 0
    p = 0;
end

Info.h2 = 0;
%capturing the solutions (if ones exist)
if abs(p) < 1e-7 %adopted precision for the minimum primal residual
    Info.Pk = double(Pk);    
    Info.h2 = sqrt(double(gamma));
end
