function Info = stab_OP10_AUT_T3_yal(Ai,L,d,b)
% function Info = stab_OP10_AUT_T3_yal(Ai,L,d,b)
%
% Test the robust stability of the a discrete-time time-varying polytopic system where the parameters have
% bounded rates of variation. The test is formulated in terms of the existence of a quadratic in the state
% Lyapunov function with a path-dependent matrix of arbitrary path length. The approach relies on the change 
% of variables (alpha(k), delta-alpha(k+r),r=1,...,L) to a higher dimension space called gamma-space. The test
% is parametrized in terms of path length of the path-dependent Lyapunov matrix, the number of Pólya's
% relaxations and the bound of the parameters variation. The condition is published in Theorem 3 of (Oliveira
% & Peres, Automatica, 2010). 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])
%          L       -> path length of the path-dependent Lyapunov matrix
%          d       -> number of Polya's relaxations
%          b       -> parameters variation bound (0 -> time-invariant, 1 -> arbitrarily fast)
%
% 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 P
%          Info.V        -> 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];
% L=1;
% d = 0;
% Info = stab_OP10_AUT_T3_yal(Ai,L,d,0.5)
% 
% Date: 26/06/2008
% 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;

Gam    = gen_base_d_path(vertices,b,L);
for i=1:L
    Gam_now{i} = Gam(:,(i-1)*vertices+1:i*vertices)';
end
Gam_plus{1} = Gam(:,1:vertices)' + Gam(:,vertices+1:2*vertices)';
for i=2:L
    Gam_plus{i} = Gam(:,(i-1+1)*vertices+1:(i+1)*vertices)';
end
M      = size(Gam,1);
poly.N = vertices;
poly.L = L;

poly.paths = all_paths(poly.N,poly.L);

for i=1:size(poly.paths,1)
    poly.coef{i} = sdpvar(order,order,'symmetric');    
end

poly_P     = alpha2gamma_path(poly,Gam_now);
poly_Pplus = alpha2gamma_path(poly,Gam_plus);
for i=1:size(poly_P.mon,1)
    mon_string = string_monomial(poly_P.mon(i,:));
    eval(['P',mon_string,' = poly_P.coef{i};']);
    eval(['Pb',mon_string,' = poly_Pplus.coef{i};']);
end


% construct the vertices of the polytopes A(a) and P(a)
% in the multidimensinal space (now A(gamma) and P(gamma)).
for j=1:M
    Ab{j} = zeros(order);
    for i=1:vertices
        Ab{j} = Ab{j} + Gam_now{1}(i,j)*Ai(:,(i-1)*order+1:i*order);
    end    
end

deg=poly.L;
Kfr_d = gen_coefs(M,d);
size_Kfr_d = size(Kfr_d,1);
Kfr_g_d_1 = gen_coefs(M,deg + d + 1);
number_LMIs = size(Kfr_g_d_1,1);

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) == M
            Cfr  = factorial(d)/prod(factorial(kprime));
            for p=1:M
                if k(p) > kprime(p)
                    e = zeros(1,M);
                    e(p) = 1;

                    A = Ab{p};
                    mon_string = string_monomial(k-kprime-e);
                    P = eval(['P',mon_string]);
                    Pb = eval(['Pb',mon_string]);

                    Term = Term + Cfr*([P A'*Pb;Pb*A Pb]);
                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.V = 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:size(poly_P.mon,1)
        mon_string = string_monomial(poly_P.mon(i,:));
        P=eval(['P',mon_string]);
        Info.Pi = [Info.Pi double(P)];        
    end
    Info.feas = 1;
end

%_____________________________________________________________________________________________________________
function Gam = gen_base_d_path(N,b,L)

for i = 1:N
    Gam((i-1)*N+1,i) = 1;
    putj    = 1:N;
    putj(i) = [];
    indj    = 2;
    for j = putj,
        Gam((i-1)*N+indj,[i N+i N+j]) = [ 1 -b b];
        Gam(N^2+(i-1)*(N-1)+indj-1,[i j N+i N+j]) = [ b 1-b -b b ];
        indj = indj+1;
    end
end

if nargin == 3 & L > 1
    Gam = compute_sucessive_deltas(N,b,Gam,L);
end

%_____________________________________________________________________________________________________________
function new_Gam = compute_sucessive_deltas(N,b,Gam,L)

new_Gam = Gam;

for ell=2:L
    next_Delta = new_Gam(:,1:N);
    for j=2:ell
        next_Delta = next_Delta + new_Gam(:,(j-1)*N+1:j*N);
    end
    next = [];
    for i=1:size(next_Delta,1)
        gnow = new_Gam(i,:);
        previous = find_previous(Gam,N,next_Delta(i,:));
        for j=1:size(previous,1)
            next = [next;gnow previous(j,:)];
        end
    end    
    new_Gam=next;
end


%________________________________________________________________________________________________________
function previous = find_previous(Gam,N,gnow)

previous = [];
for i=1:size(Gam,1)
    g=Gam(i,1:N);
    if sum(g==gnow) == size(gnow,2)
        previous = [previous;Gam(i,N+1:2*N)];
    end
end

%_________________________________________________________________________________________
function Paths = all_paths(N,L)

Paths = [];
for i=1:L    
    Paths = [Paths kron(ones(1,N^(i-1)),kron(1:N,ones(1,N^(L-i))))'];
end
%___________________________________________________________________________________________________
function [string] = string_monomial(monomial)

string = '';
for i=1:size(monomial,2)
    string = strcat(string,'_',int2str(monomial(1,i)));
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
