function [Info] = estab_OP06_SCL_T3_yal(Ai,deg)
%function [Info] = estab_OP06_SCL_T3_yal(Ai,deg)
% test the Schur stability of the polytope Ai via a homogeneous
% polynomially parameter-dependent Lyapunov function of arbitray degree. 
% Theorem 3 of (Oliveira & Peres, SCL, 2006)
%input:  Ai      -> vertices of the polytope (Ai = [A1 A2 ... A3])
%        deg     -> degree of the polynomially parameter dependent Lyapunov function
%output: 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 P
%        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_OP06_SCL_T3_yal(Ai,1); %affine parameter-dependent Lyapunov matrix
%
% 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;

%generate the exponents of the monomials to create the variables P_k
K_g = gen_coefs(vertices,deg);
number_VAR = size(K_g,1);%number of variables

%generate the variables
for i=1:number_VAR
    mon_string = string_monomial(K_g(i,:));
    eval(['P',mon_string,' = sdpvar(order,order,''symmetric'');']);
    P = eval(['P',mon_string]);
    LMIs = LMIs + set(P > 0);
    Info.L = Info.L + order;
end

Kfr_g_1 = gen_coefs(vertices,deg + 1);
number_LMIs = size(Kfr_g_1,1);

for i=1:number_LMIs
    
    monomio_atual = Kfr_g_1(i,:);
    nonzero = find(monomio_atual);%nonzero exponents 
    n_nonzero = size(nonzero,2);%number of nonzeros exponents (determine the number of terms in the LMI)
    Term = zeros(2*order);    
    for j=1:n_nonzero        
        ind = nonzero(j);
        A = Ai(:,order*(ind-1)+1:ind*order);
        mono = monomio_atual;
        mono(ind) = mono(ind) - 1;
        
        mon_string = string_monomial(mono);
        P = eval(['P',mon_string]);
        Term = Term + [-P P*A;A'*P -P];
    end
    Info.L = Info.L + 2*order;
    LMIs = LMIs + set(Term < 0);
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:number_VAR
        mon_string = string_monomial(K_g(i,:));
        Info.Pi = [Info.Pi double(eval(['P',mon_string]))];
    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
    for i=0:vertices
        var(i+1)=0;
    end
    comb = [];
    comb = nomial(degree, degree, vertices, 1, var,comb);
end

%___________________________________________________________________________________________________
function Coefs = nomial(n,n2,v,v2,var,Coefs)

if v2 < v
    for i=0:n2
        var(v2+1) = i;
        Coefs = nomial(n, n2-i, v, v2+1, var,Coefs);
    end
else
    var(v+1)=n2;
    prod = part(n,var,v);
    este = [];
    for j=1:v
        este = [este var(j+1)];
    end
    Coefs = [Coefs;este];
end

%___________________________________________________________________________________________________
function prod = part(n,var,v)

prod = 1;
for i=1:v
    prod = prod * factorial(var(i+1));
end
prod = factorial(n)/prod;

%___________________________________________________________________________________________________
function [string] = string_monomial(monomial)

string = '';
for i=1:size(monomial,2)
    string = strcat(string,'_',int2str(monomial(1,i)));
end