function Info = estab_OP07_IEEETAC_T4_c_yal(Ai,deg,d)
%function Info = estab_OP07_IEEETAC_T4_c_yal(Ai,deg,d)
% test the Hurwitz stability of the polytope Ai via a homogeneous
% polynomially parameter-dependent Lyapunov function of arbitray degree
% using Polya's relaxations. Parameter-dependent slack variables are also used. 
%Theorem 4 of (Oliveira & Peres, IEEE-TAC 2007)
%input:  Ai      -> vertices of the polytope (Ai = [A1 A2 ... A3])
%        deg     -> degree of the polynomially parameter dependent Lyapunov function
%        d       -> number of relaxations
%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.9  0.2;   
%       -0.5  -1.9 ];  
% A2 = [-1.1  -0.1;
%       0.5    -2.8];
% Ai = [A1 A2];
% Info = estab_OP07_IEEETAC_T4_c_yal(Ai,1,0); %affine parameter-dependent Lyapunov matrix
%
% Date: 13/03/2007
% 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]);
    eval(['Q',mon_string,'=[zeros(order) P;P zeros(order)];']);
    eval(['X',mon_string,'=sdpvar(2*order,order,''full'');']);
end

Kfr_g_d = gen_coefs(vertices,deg + d);
Kfr_d = gen_coefs(vertices,d);
number_LMIs = size(Kfr_g_d,1);
size_Kfr_d = size(Kfr_d,1);

% P(alpha) > 0
for j=1:number_LMIs

    k = Kfr_g_d(j,:);

    Term = zeros(order);
    for m=1:size_Kfr_d
        kprime = Kfr_d(m,:);
        if size(find(k-kprime>=0),2) == vertices
            Cfr  = factorial(d)/Pi(kprime);

            P = eval(['P',string_monomial(k-kprime)]);%get the variable index
            Term = Term + Cfr*P;
        end
    end
    LMIs = LMIs + set(Term > 0);
    Info.L = Info.L + order;
end

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

% A(alpha)'P(alplha) + P(alpha)A(alpha) < 0
for i=1:number_LMIs

    k = Kfr_g_d_1(i,:);

    Term = zeros(2*order);
    for j=1:size_Kfr_d
        kprime = Kfr_d(j,:);
        if size(find(k-kprime>=0),2) == vertices
            Cfr  = factorial(d)/Pi(kprime);
            for p=1:vertices
                if k(p) > kprime(p)
                    e = zeros(1,vertices);
                    e(p) = 1;

                    A = Ai(:,(p-1)*order+1:p*order);
                    Bfr = [A -eye(order)];
                    mon_string = string_monomial(k-kprime-e);
                    Q = eval(['Q',mon_string]);
                    X = eval(['X',mon_string]);

                    Term = Term + Cfr*[Q + X*Bfr + Bfr'*X'];
                end
            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:number_VAR
        mon_string = string_monomial(K_g(i,:));
        Info.Pi = [Info.Pi double(eval(['P',mon_string]))];
    end
    Info.feas = 1;
end

%___________________________________________________________________________________________________
function R = Pi(monomial)

R=prod(factorial(monomial));

%___________________________________________________________________________________________________
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

%___________________________________________________________________________________________________
function [string] = string_monomial(monomial)

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