function Info = estab_OP05_LAA_T3_yal(Ai,d)
% function Info = estab_OP05_LAA_T3_yal(Ai,d)
% test the robust Hurwitz stability of the polytope of matrices using an affine parameter-
% dependent Lyapunov function and Polya's relaxations. The condition is given in Theorem 3
% (Oliveira & Peres, LAA, 2005). 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 ... A3])
%         d       -> number of Polya's relaxations
%
% outputs: Info.feas     -> feasible (1) or not (0)
%          Info.cpusec   -> cpu time to solve the LMIs (seconds)
%          Info.cpusec_m -> cpu time to mount the LMIs (seconds)
%          Info.Pi       -> Solution variables Pi
%          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 = estab_OP05_LAA_T3_yal(Ai,1)
% 
% Date: 09/08/2005
% Author: ricfow@dt.fee.unicamp.br

%determine the number of vertices of polytope
[order,vertices] = size(Ai);
vertices = vertices / order;

Info.cpusec_m = clock;

%new LMI system
LMIs = set([]);
Info.L = 0;

for i=1:vertices
    P{i} = sdpvar(order,order,'symmetric');
    LMIs = LMIs + set(P{i} > 0);
    Info.L = Info.L + order;
end

combinations = gen_coefs(vertices,d + 2);
number_COMB = size(combinations,1);%number of combinations of p_1 + p_2 + ... + p_N = n

for k=1:number_COMB
    comb_now = combinations(k,:);
    
    Term = zeros(order);
    for i=1:vertices
        Aii = Ai(:,(i-1)*order+1:i*order);
        comb_aux = comb_now;
        comb_aux(1,i) = comb_aux(1,i) - 2;
        if comb_aux(1,i) >= 0
            C = factorial(d)/prod(factorial(comb_aux));
            Term = Term + C*(Aii'*P{i} + P{i}*Aii);            
        end        
        for j=i+1:vertices
            Ajj = Ai(:,(j-1)*order+1:j*order);                
            comb_aux = comb_now;
            comb_aux(1,i) = comb_aux(1,i) - 1;
            comb_aux(1,j) = comb_aux(1,j) - 1;            
            if comb_aux(1,i) >= 0 & comb_aux(1,j) >= 0
                C = factorial(d)/prod(factorial(comb_aux));
                Term = Term + C*(Aii'*P{j} + P{j}*Aii + Ajj'*P{i} + P{i}*Ajj);
            end            
        end
    end    
    LMIs = LMIs + set(Term < 0);
    Info.L = Info.L + 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.Pi   = [];
    for i=1:vertices        
        Info.Pi = [Info.Pi double(P{i})];
    end
    Info.feas = 1;
end

%___________________________________________________________________________________________________
function [comb] = gen_coefs(vertices,degree)
%generates all the solutions of: a_1 + a_2 + ... + a_vertices = degree

name = strcat('comb_N',int2str(vertices),'d',int2str(degree),'.mat');
if exist(name,'file') == 2
    load(name);
else
    comb = [];
    if vertices == 1
        comb = degree;
    else
        for i = 0:degree
            temp = gen_coefs(vertices-1,degree-i);
            comb   = [ comb
                i*ones(size(temp,1),1) temp ];
        end
    end
end
