function Info = stab_OdOP08_SCL_T8_yal(Ai,kappa)
% function Info = stab_OdOP08_SCL_T8_yal(Ai,kappa)
% Test the Hurwitz robust stability of a continuous-time time-invariant polytopic linear 
% system via the notion of kappa-stability presented in Theorem 8 of (Oliveira, de Oliveira 
% & Peres, SCL 2008).  The LMIs are programmed using YALMIP and can be solved by any LMI 
% solver supported by YALMIP (SeDuMi is the default).
%
% inputs:  Ai      -> vertices of the polytope (Ai = [A1 A2 ... AN])
%          kappa   -> value of kappa  (k >= 1)
%
% outputs: Info.feas   -> stable (1) or not (0)
%          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
%          Info.L      -> number of LMI rows
%
% 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];
% Info = stab_OdOP08_SCL_T8_yal(Ai,2)
%
% Date: 20/06/2007
% Author: ricfow@dt.fee.unicamp.br

%dimensions of the system
[order,vertices] = size(Ai);
vertices = vertices / order;

Info.cpusec_m = clock;
%Initializes the description of a new LMI system.
LMIs = set([]);
%LMI counter
Info.L = 0;

Rk = [zeros(kappa,1) eye(kappa)];
Lk = [eye(kappa) zeros(kappa,1)];

Pk = sdpvar(kappa*order,kappa*order,'symmetric');
if kappa == 1
    LMIs = LMIs + set(Pk > 0);
    Info.L = Info.L + kappa*order;
end
Qk = kron(Lk',eye(order))*Pk*kron(Rk,eye(order)) + kron(Rk',eye(order))*Pk*kron(Lk,eye(order));

%creating the variables
for i=1:vertices
    if kappa > 1
        Y{i} = sdpvar(kappa*order,(kappa-1)*order,'full');
    end
    Z{i} = sdpvar((kappa+1)*order,kappa*order,'full');
end

%programming the LMIs
for i=1:vertices
    Aii = Ai(:,(i-1)*order+1:i*order);
    
    if kappa > 1        
        %(P_k)_i > Y(C_k-1)_i + (C_k-1)_i'Y'
        Ck1_i = kron([eye(kappa-1) zeros(kappa-1,1)],Aii) - kron([zeros(kappa-1,1) eye(kappa-1)],eye(order));
        LMIs = LMIs + set(Pk > Y{i}*Ck1_i + Ck1_i'*Y{i}');        
        Info.L = Info.L + kappa*order;
    end

    %(Q_k)i  + Z(C_k)i + (C_k)i'Z' < 0
    Ck2_i = [kron(eye(kappa),Aii) zeros(order*kappa,order)] - [zeros(order*kappa,order) eye(kappa*order)];
    LMIs = LMIs + set(Qk + Z{i}*Ck2_i + Ck2_i'*Z{i}' < 0);    
    Info.L = Info.L + (kappa+1)*order;    

    for j=i+1:vertices
        Ajj = Ai(:,(j-1)*order+1:j*order);
        if kappa > 1
            Ck1_j = kron([eye(kappa-1) zeros(kappa-1,1)],Ajj) - kron([zeros(kappa-1,1) eye(kappa-1)],eye(order));
            LMIs = LMIs + set(2*Pk > Y{i}*Ck1_j + Ck1_j'*Y{i}' + Y{j}*Ck1_i + Ck1_i'*Y{j}');
            Info.L = Info.L + kappa*order;
        end
        Ck2_j = [kron(eye(kappa),Ajj) zeros(order*kappa,order)] - [zeros(order*kappa,order) eye(kappa*order)];        
        LMIs = LMIs + set(2*Qk + Z{i}*Ck2_j + Ck2_j'*Z{i}' + Z{j}*Ck2_i + Ck2_i'*Z{j}' < 0);
        Info.L = Info.L + (kappa+1)*order;        
    end
end

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

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

sol = solvesdp(LMIs,[],sdpsettings('verbose',0,'solver','sedumi'));
Info.cpusec = sol.solvertime;
[p,d]=checkset(LMIs);

Info.feas = 0;
%capturing the solutions (if ones exist)
if sum(p<0) == 0    
    Info.Pk   = double(Pk);    
    Info.feas = 1;
end
