function Info = estab_OP05_LAA_T5_d_yal(Ai,d)
% function Info = estab_OP05_LAA_T5_d_yal(Ai,d)
% test the robust Schur stability of the polytope of matrices using an affine parameter-
% dependent Lyapunov function and Polya's relaxations. The condition is given in Theorem 5
% (Oliveira & Peres, LAA, 2005). Slack variables are used to improve the positivity test.
% 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.1  1.1;   
%       -0.8  -0.3];  
% A2 = [-0.6  0.6;
%       0.4   -0.3];
% Ai = [A1 A2];
% 
% Info = estab_OP05_LAA_T5_d_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;
    F{i} = sdpvar(order,order,'full');
    G{i} = sdpvar(order,order,'full');    
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(2*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));
            T11 = F{i}*Aii + Aii'*F{i}' - P{i};
            T12 = -F{i}+Aii'*G{i}';
            T22 = P{i}-G{i}-G{i}';
            Term = Term + C*[T11 T12;T12' T22];
        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));
                T11 = F{i}*Ajj + Ajj'*F{i}' + F{j}*Aii + Aii'*F{j}'-P{i}-P{j};
                T12 = -F{i}+Aii'*G{j}' -F{j}+Ajj'*G{i}';
                T22 = P{i}+P{j}-G{i}-G{i}' -G{j}-G{j}';
                Term = Term + C*[T11 T12;T12' T22];               
            end            
        end
    end    
    LMIs = LMIs + set(Term < 0);
    Info.L = Info.L + 2*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
