function Info = design_OP10_AUT_T4_yal(Ai,Bi,deg,d,b)
% function Info = design_OP10_AUT_T4_yal(Ai,Bi,deg,d,b)
%
% Design a robust stabilizing state-feedback gain for 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 homogeneous polynomial parameter-dependent matrix. The approach relies 
% on the change of variables (alpha and delta-alpha) to a higher dimension space called gamma-space. The test 
% is parametrized in terms of the degree of the homogeneous polynomial Lyapunov matrix, the number of Pólya's
% relaxations and the bound of the parameters variation. The condition is published in Theorem 4 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 dynamic matrix (Ai = [A1 A2 ... A3])
%          Bi      -> vertices of the input matrix (Bi = [B1 B2 ... B3])
%          deg     -> degree of the polynomially parameter dependent Lyapunov function
%          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
%          Info.K        -> robust state-feedback gain
%          Info.Acl      -> closed-loop dynamic matrix 
%
% Example: 2 states and 2 vertices
% A1=[-0.0216   -0.2351;
%     0.8987   -1.2888];
% A2=[0.6780   -0.3487;
%    -0.4137   -1.0098];
% Ai = [A1 A2];
% Bi = [1 2;1 1];
% deg=1;
% d = 0;
% Info = design_OP10_AUT_T4_yal(Ai,Bi,deg,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;
inputs = size(Bi,2)/vertices;

Info.cpusec_m = clock;

%new LMI system
LMIs = set([]);
Info.L = 0;
Gam = gen_base_d(vertices,b);
   
Gam_1 = Gam(:,1:vertices)';
Gam_2 = Gam(:,vertices+1:end)';
M=size(Gam_1,2);

polyP.mon = gen_coefs(vertices,deg);
number_VAR = size(polyP.mon,1);%number of variables

%generate the variables
for i=1:number_VAR
    polyP.coef{i} = sdpvar(order,order,'symmetric');    
end
polyP.N = vertices;
polyP.deg = deg;
C1 = alpha2gamma_HPPD(polyP,Gam_1);
C2 = alpha2gamma_HPPD(polyP,Gam_1+Gam_2);
for i=1:size(C1.mon,1)
    mon_string = string_monomial(C1.mon(i,:));
    eval(['P',mon_string,'  = C1.coef{i};']);
    eval(['Pb',mon_string,' = C2.coef{i};']);
end
G = sdpvar(order,order,'full');
Z = sdpvar(inputs,order,'full');

Ab = [];
Bb = [];
for j=1:M
    Anow = zeros(order);
    Bnow = zeros(order,inputs);
    for i=1:vertices
        Anow = Anow + Gam_1(i,j)*Ai(:,(i-1)*order+1:i*order);
        Bnow = Bnow + Gam_1(i,j)*Bi(:,(i-1)*inputs+1:i*inputs);
    end    
    Ab = [Ab Anow];
    Bb = [Bb Bnow];
end


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

for j=1:number_LMIs

    k = Kfr_g_d(j,:);

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

            Term_affine = zeros(2*order);
            for p=1:M
                if k(p)-kprime(p) > 0  
                    e = zeros(1,M);
                    e(p) = 1;
                    beta=factorial(deg-1)/prod(factorial(k-kprime-e));
                    A = Ab(:,(p-1)*order+1:p*order);
                    B = Bb(:,(p-1)*inputs+1:p*inputs);
                    T12 = beta*(A*G + B*Z);
                    T22 = beta*(G+G');
                    Term_affine = Term_affine + [zeros(order) T12;T12' T22];
                end
            end            
            mon_string = string_monomial(k-kprime);
            P = eval(['P',mon_string]);            
            Pb = eval(['Pb',mon_string]);               
            Term = Term + Cfr*(Term_affine + [Pb zeros(order);zeros(order) -P]);
        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','sedumi.eps',1e-14));
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
        Info.Pi = [Info.Pi double(polyP.coef{i})];
    end
    Info.feas = 1;
    Info.Ki = double(Z)*inv(double(G));
    Info.Acl = [];
    for i=1:vertices
        Info.Acl = [Info.Acl Ai(:,(i-1)*order+1:i*order)+Bi(:,(i-1)*inputs+1:i*inputs)*Info.Ki];
    end
end

%_____________________________________________________________________________________________________________
function Gam = gen_base_d(N,b)

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

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