function Info = stab_OBP08_CDC_T2_yal(MSAi,deg)
% function Info = stab_OBP08_CDC_T2_yal(MSAi,deg)
% evaluates the Hurwitz robust stability of an uncertain linear system whose dynamic matrix
% is a homogeneous polynomial matrix with parameters belonging to the cartesian product of
% many simplexes, called multi-simplex. The test is performed through the existence of a 
% homogeneous polynomial matrix with parameters in the multi-simplex. The partial degrees 
% of the Lyapunov matrix, associated to each simplex inside the multi-simplex, can be 
% freely chosen by the user. The condition is given in Theorem 2 of (Oliveira, Bliman &
% Peres, CDC08, to appear). The LMIs are programmed using YALMIP 
% and can be solved by any LMI solver supported by YALMIP (SeDuMi is the default). 
%
% inputs: MSAi    -> homogeneous polynomial matrix in the multi-simplex representation to
%                    be tested. (see the routine affine2MS.m that converts an affine model 
%                    to a multi-simplex model)
%         deg     -> partial degrees [g1 g2 ... gm] associated to the polynomially 
%                    parameter-dependent Lyapunov matrix.
%
% outputs: Info.feas     -> stable (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 Pk
%          Info.K        -> number of scalar variables
%          Info.L        -> number of LMI rows
%
% Example: affine model A(theta)=A0 + A_1theta_1, |theta_1|< 1,   with
% A0 = [-4 2 -2;5 -6 1; -2 2 -7];
% A1 = [-5 -3  -13;-5 0 0;10 13 16];
% Atheta = [A0 A1];
% bounds = [-1 1];
% MSAi = affine2MS_LTI(Atheta,bounds);
% partial_degs = [4 4];
% Info = stab_OBP08_CDC_T2_yal(MSAi,partial_degs)
%
% Date: 18/02/2008
% Author: ricfow@dt.fee.unicamp.br

order = MSAi.nxm(1);
N = MSAi.m;

Info.cpusec_m = clock;

%new LMI system
LMIs = set([]);
%LMI rows counter
Info.L = 0;

%generate sets Kcal_i(g_i)
for i=1:N
    mon{i} = gen_coefs(MSAi.N(1,i),deg(1,i))';
end
K_g = all_possible_comb(N,mon);
number_VAR = size(K_g,2);
%generate the variables associated to the Lyapunov matrix
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

K_1 = MSAi.mon;
size_K_1 = size(K_1,2);

for i=1:N
    mon{i} = gen_coefs(MSAi.N(1,i),deg(1,i)+MSAi.deg(1,i))';
end
%generate the set Kcal(g+r)
Kfr_g_plus_r = all_possible_comb(N,mon);
number_LMIs = size(Kfr_g_plus_r,2);%number of variables

% A_MS(alpha)'P_MS(alpha) + P_MS(alpha)A_MS(alpha) < 0
for i=1:number_LMIs

    k = Kfr_g_plus_r(:,i);

    Term = zeros(order);
    for p=1:size(K_1,2)
        k_1 = K_1(:,p);
        if sum((k-k_1) >= 0) == size(k,1) %no negative component
        

            A = MSAi.coef{p}.coef_val;
            mon_string = string_monomial((k-k_1)');
            P = eval(['P',mon_string]);

            Term = Term + A'*P + P*A;            
        end
    end    
    LMIs = LMIs + set(Term < 0);
    Info.L = Info.L + 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

% Various auxiliary routines 
%
%________________________________________________________________________________________________________
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 Resp = all_possible_comb(n,H)

var = [];
Resp = call_recur(n,H,1,var,[],1);

%________________________________________________________________________________________________________
function Resp = call_recur(n,H,na,var,Resp,pos_ini)

if na <= n
    tam = size(H{na},1);
    for i=1:size(H{na},2)
        var(pos_ini:pos_ini+tam-1,1) = H{na}(:,i);
        Resp = call_recur(n,H,na+1,var,Resp,pos_ini+tam);
    end
else
    Resp = [Resp var];
end

%________________________________________________________________________________________________________
function string = string_monomial(monomial)

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