function Info = stab_MOPB09_IJC_C4_c_yal(Ai,deg,d,dot_alphas)
% function Info = stab_MOPB09_IJC_C4_c_yal(Ai,deg,d,dot_alphas)
% test the robust stability of a continuous-time time-varying polytopic linear system
% whose parameters have bounded rates variation. The test is formulated through the 
% existence of a Lyapunov function quadratic in the state and homogeneous polynomially 
% of arbitray degree in the parameters. The condition is given in Corollary 4 of (Montagner,
% Oliveira, Peres & Bliman, IJC, 2009). Slack variables dependending affinely on
% on \alpha and depending affinely on \eta are introduced in the problem 
% to improve convergence, at the price of a higher computational burden. Polya's 
% relaxations can also be used. The LMIs are programmed using YALMIP  and can be solved 
% by any LMI solver supported by  YALMIP (SeDuMi is the default). 
%
% input: Ai      -> vertices of the polytopes (Ai = [A1 A2 ... AN])
%        deg     -> degree of the polynomially parameter-dependent Lyapunov matrix
%        d       -> number of Pólya's relaxations
%        dot_alphas-> time-derivative bounds of the parameters given in the form:
%                        dot_alphas = [dot_alpha_1_min dot_alpha_1_max;
%                                      dot_alpha_2_min dot_alpha_2_max;
%                                           ..              ..
%                                      dot_alpha_N_min dot_alpha_N_max];
%
% output: 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.K        -> number of scalar variables
%        Info.L        -> number of LMIs rows
%        Info.Pk       -> solution variables Pk
%
% Example: 2 states, 2 vertices, 1 input
% A1 = [-0.9  0.2;
%       -0.5  -1.9 ];
% A2 = [-1.1  -0.1;
%       0.5    -2.8];
% Ai = [A1 A2];
% dot_alphas = [-5 5;-7 7];

% Info = stab_MOPB09_IJC_C4_c_yal(Ai,1,0,dot_alphas)
% 
% Date: 07/08/2007
% Author: ricfow@dt.fee.unicamp.br

if deg < 2
    deg = 2;
end
%determine the number of vertices of polytope
[order,vertices] = size(Ai);
vertices = vertices / order;

Info.cpusec_m = clock;

%new LMI system
[H] = gen_base(vertices,dot_alphas);

LMIs = set([]);
Info.L = 0;

%generate the exponents of alpha's (power 'deg') to create the variables P
coefs_VAR = gen_coefs(vertices,deg);
number_VAR = size(coefs_VAR,1);%number of variables

%generate the variables
for i=1:number_VAR
    mon_string = string_monomial(coefs_VAR(i,:));
    eval(['P',mon_string,' = sdpvar(order,order,''symmetric'');']);
end

coefs_VAR_affine = gen_coefs(vertices,1);
%generate the multipliers
for i=1:size(coefs_VAR_affine,1);
    mon_string = string_monomial(coefs_VAR_affine(i,:));
    for j=1:size(H,2)
        eval(['X',mon_string,'{j} = sdpvar(2*order,order,''full'');']);
    end
end

Kfr_g_d = gen_coefs(vertices,deg+d);
Kfr_d = gen_coefs(vertices,d);
size_Kfr_d = size(Kfr_d,1);
number_LMIs = size(Kfr_g_d,1);%number of LMIs to be generated
% positive definite Lyapunov matrix
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_minus_1 = gen_coefs(vertices,deg - 1);
size_Kfr_g_minus_1 = size(Kfr_g_minus_1,1);
Kfr_g_d_minus_2 = gen_coefs(vertices,deg  +d -2);
size_Kfr_g_d_minus_2 = size(Kfr_g_d_minus_2,1);

for ell=1:size(H,2)
    for i=1:number_LMIs

        k = Kfr_g_d(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);
                mon_string = string_monomial(k-kprime);
                P = eval(['P',mon_string]);
                Term_invariant = [zeros(order) P;
                    P zeros(order)];

                Term_varying = zeros(2*order);
                for m=1:vertices
                    em = zeros(1,vertices);
                    em(m) = 1;
                    Term_P = zeros(order);
                    for o=1:size_Kfr_g_minus_1
                        kbar = Kfr_g_minus_1(o,:);
                        if size(find(k-kprime>=kbar),2) == vertices
                            new_k = kbar+em;
                            mon_string = string_monomial(new_k);
                            P = eval(['P',mon_string]);
                            Term_P = Term_P + new_k(1,m)*P;
                        end
                    end
                    Term_varying = Term_varying + H(m,ell)*[Term_P zeros(order); zeros(order) zeros(order)];
                end
                Term = Term + Cfr*(Term_invariant + Term_varying);
            end
        end
        %just the quadratic part (X(alpha)B(alpha) + B(alpha)'X(alpha)')
        for j=1:size_Kfr_g_d_minus_2
            kprime = Kfr_g_d_minus_2(j,:);
            if size(find(k-kprime>=0),2) == vertices
                Cfr  = factorial(deg+d-2)/Pi(kprime);
                for p=1:vertices
                    e = zeros(1,vertices);
                    e(p) = 1;

                    if k(p) > kprime(p)
                        Bfr = [Ai(:,(p-1)*order+1:p*order) -eye(order)];
                        mon_string = string_monomial(k-kprime-e);
                        X = eval(['X',mon_string,'{ell}']);
                        Term = Term + Cfr*(X*Bfr + Bfr'*X');
                    end
                end
            end
        end
        LMIs = LMIs + set(Term < 0);
        Info.L = Info.L + 2*order;        
    end    
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.Pk   = [];
    for i=1:number_VAR
        mon_string = string_monomial(coefs_VAR(i,:));
        Info.Pk = [Info.Pk 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

%___________________________________________________________________________________________________
function [H] = gen_base(N,bounds)
%function [H] = gen_base(N,bounds)
% This function generates the representation base H  used by the diffential
% inclusion method
% inputs: N       -> number of vertices of the polytope
%         bounds  -> time-derivative bounds. ex.: [[-1 1];[-1 1]...[-1 1]];
%
% output: H -> representation base

H = solve_plane(N,bounds,0);

%___________________________________________________________________________________________________
function H = solve_plane(N,limites,plano)

if N > 3
    Hinter1 = solve_plane(N-1,limites(2:end,:),limites(1,1)+plano);
    Hinter2 = solve_plane(N-1,limites(2:end,:),limites(1,2)+plano);
    cima1 = [];
    cima2 = [];
    for i=1:size(Hinter1,2)
        cima1 = [cima1 limites(1,1)];
    end
    for i=1:size(Hinter2,2)
        cima2 = [cima2 limites(1,2)];
    end
    Hinter1 = [cima1;Hinter1];
    Hinter2 = [cima2;Hinter2];
    H = [Hinter1 Hinter2];
else
    H = [];
    for i=1:N
        for k=1:N
            if i == k
                continue;
            end
            sobra = [];
            pos = [];
            if N == 2
                if -limites(i,1) >= limites(k,1) & -limites(i,1) <= limites(k,2)
                    vec_final(1,i) = limites(i,1);
                    vec_final(1,k) = -limites(i,1);
                    if ja_contem(H,vec_final') == 0
                        H = [H vec_final'];
                    end
                end
                if -limites(i,2) >= limites(k,1) & -limites(i,2) <= limites(k,2)
                    vec_final(1,i) = limites(i,2);
                    vec_final(1,k) = -limites(i,2);
                    if ja_contem(H,vec_final') == 0
                        H = [H vec_final'];
                    end
                end
            else
                for o=1:N
                    if o ~= i & o ~= k
                        pos = o;
                        sobra = limites(o,:);
                        break;
                    end
                end
                for m=1:size(sobra',1);
                    soma1 = -sum(sobra(:,m)) - limites(i,1) - plano;
                    soma2 = -sum(sobra(:,m)) - limites(i,2) - plano;
                    if soma1 >= limites(k,1) & soma1 <= limites(k,2)
                        vec_final(1,i) = limites(i,1);
                        vec_final(1,k) = soma1;
                        vec_final(1,pos) = sobra(:,m);

                        if ja_contem(H,vec_final') == 0
                            H = [H vec_final'];
                        end
                    end
                    if soma2 >= limites(k,1) & soma2 <= limites(k,2)
                        vec_final(1,i) = limites(i,2);
                        vec_final(1,k) = soma2;
                        vec_final(1,pos) = sobra(:,m);
                        if ja_contem(H,vec_final') == 0
                            H = [H vec_final'];
                        end
                    end
                end
            end
        end
    end
end
%___________________________________________________________________________________________________
function R = ja_contem(H,h)

R = 0;
for i=1:size(H,2)
    if sum(H(:,i) == h) == size(h,1)
        R = 1;
        return;
    end
end