function Info = stabLPV_OdOP_IETCTA_C1b_yal(Ai,kappa,H)
%function Info = stabLPV_OdOP_IETCTA_C1b_yal(Ai,kappa,H)
% test the robust stability of a time-varying polytopic system using the notion of kappa-robust
% stability adapted to the time-varying case. The conditions take into account informations about
% up the (kappa-1)-th time-derivative of the uncertain parameters. The slack variables (multipliers)
% are considered alpha-independent. This condition is formulated using the choice C1 given
% in the paper (Oliveira, de Oliveira & Peres, IETCTA, submitted) and by considering the relation
% dot{alpha} = H alpha, i.e. the time-varying parameters have an exponential dependece on time. This
% condition is indicated as bar{C}1 in the paper.
%
% inputs: Ai    -> vertices of the polytope (Ai = [A1 A2 ... AN])
%         kappa -> value of kappa
%         H     -> dynamic matrix describing the differential equation: dot{alpha} = H alpha
%                      
%
% outputs: Info.feas     -> 1 if robustly stable, 0 otherwise
%          Info.cpusec   -> cpu time to solve the LMIs (seconds)
%          Info.cpusec_m -> cpu time to mount the LMIs (seconds)
%          Info.Pk       -> Solution variables Pk
%          Info.K        -> number of scalar variables used in the optimization problem
%          Info.L        -> number of LMI rows used in the optmization problem
%
%
% Example: 2 states and 2 vertices
% beta = 0.22;
% A1 = [0 1; -(1 + 3*beta^2) -0.05];
% A2 = [0 1; -(1 - beta^2) -0.05];
% A3 = A2;
% Ai = [A1 A2 A3];
% omega = 2;
% q=omega/4;
% H = [q -3*q q; 3*q -q -q; -4*q 4*q 0];
% Info = stabLPV_OdOP_IETCTA_C1b_yal(Ai,3,H)
%
% Date: 08/05/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([]);
Info.L = 0;

if kappa > 1
    Y = sdpvar(kappa*order, (kappa-1)*order,'full');
end
Z = sdpvar((kappa+1)*order, kappa*order,'full');

R_k = [zeros(kappa,1) eye(kappa)];
L_k = [eye(kappa) zeros(kappa,1)];
for i=1:vertices
    Pk{i} = sdpvar(kappa*order,kappa*order,'symmetric');
    Qk{i} = kron(L_k',eye(order))*Pk{i}*kron(R_k,eye(order)) + kron(R_k',eye(order))*Pk{i}*kron(L_k,eye(order));
end

for i=1:vertices
    %(P_k)_i > Y(C_k-1)_i + (C_k-1)'Y'
    A = Ai(:,(i-1)*order+1:i*order);    
    Term = Pk{i};
    if kappa > 1
        Bkappa = generate_Bkappa(Ai,i,H,kappa-1,order,vertices);
        Term = Term - Y*Bkappa - Bkappa'*Y';
    end
    LMIs = LMIs + set(Term > 0);
    Info.L = Info.L + kappa*order;

    %(Q_k)i  + Z(C_k)i + (C_k)i'Z' < 0
    Bkappa = generate_Bkappa(Ai,i,H,kappa,order,vertices);
    dotP = zeros(kappa*order);
    for v=1:vertices
        dotP = dotP + Pk{v}*H(v,i);        
    end       
    dotP = [dotP zeros(kappa*order,order);zeros(order,kappa*order+order)];
    
    LMIs = LMIs + set(Qk{i} + dotP + Z*Bkappa + Bkappa'*Z' < 0);
    Info.L = Info.L + (kappa+1)*order;

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.feas = 1;
    Info.Pi = [];
    for i=1:vertices
        Info.Pi = [Info.Pi double(Pk{i})];
    end
end

%________________________________________________________________________________________________________
function Bfr = generate_Bkappa(Ai,vert,H,kappa,order,vertices)

Bfr = [Ai(:,(vert-1)*order+1:vert*order) -eye(order) zeros(order,(kappa-1)*order)];
for i=1:kappa-1
    B = [];
    for j=1:i
        varsigma = factorial(i)/(factorial(j-1)*factorial(i+1-j));
        xdotA = evaluate_xdotA(Ai,H,i-j+1,vert,order,vertices);
        B = [B xdotA*varsigma];
    end
    B = [B Ai(:,(vert-1)*order+1:vert*order) -eye(order) zeros(order,(kappa-i-1)*order)];
    Bfr = [Bfr; B];
end

%________________________________________________________________________________________________________
function A = evaluate_xdotA(Ai,H,kappa,i,order,vertices)

A = zeros(order);
Hp = H^kappa;
for k=1:vertices
    A = A + Hp(k,i)*Ai(:,(k-1)*order+1:k*order);
end