function Info = synth_stab_MOPB07_IET_T2_yal(Ai,B2i,deg,d)
% function Info = synth_stab_MOPB07_IET_T2_yal(Ai,B2i,deg,d)
% synthesize a homogeneous polynomially parameter-dependent gain that ensure closed-loop quadratic stability
% for continuous-time arbitrarily time-varying polytopic linear systems. Pólya's
% relaxations can be used. For completeness, it is allowed to consider homogeneous polynomial gains
% of arbitrary degree, but the degree one (deg=1) is necessary and sufficient for stabilization purposes.
% This conditions was published in (Montagner et al., IET-CTA, 2007).
%
% input: Ai,B2i         -> vertices of the polytopes (Ai = [A1 A2 ... AN], Bi = [B21 B22 ... B2N])
%        deg           -> degree of the polynomially parameter-dependent gain Z(alpha)
%        d             -> number of Pólya's relaxations
%output: Info.feas     -> 1 if such gain exits, 0 if unfeasible
%        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.W        -> solution variable W
%        Info.Ki       -> synthesized gains ([Z_1 Z_2 .. Z_N]*inv(W)]
%
% Date: 24/05/2007
% Author: ricfow@dt.fee.unicamp.br
%
%example: 2 states and 2 vertices
% A1 = [-0.9  0.2;
%       -0.5  -1.9 ];
% A2 = [-1.1  -0.1;
%       0.5    -2.8];
% B21 = [1;0];  B22 = [0;-1];
% Ai = [A1 A2]; B2i = [B21 B22];
% Info = synth_stab_MOPB07_IET_T2_yal(Ai,B2i,1,0)
% 
% Date: 24/05/2007
% Author: ricfow@dt.fee.unicamp.br

%determine the number of vertices of polytope
[order,vertices] = size(Ai);
vertices = vertices / order;
[order,inputs] = size(B2i);
inputs = inputs / vertices;

Info.cpusec_m = clock;

%new LMI system
LMIs = set([]);

Info.L = 0;
%generate the exponents of the monomials to create the variables Z_k
K_g = gen_coefs(vertices,deg);
number_VAR = size(K_g,1);%number of variables

%generate the variables Z_k
for i=1:number_VAR
    mon_string = string_monomial(K_g(i,:));
    eval(['Z',mon_string,' = sdpvar(inputs,order,''full'');']);
end
W = sdpvar(order,order,'symmetric');
LMIs = LMIs + set(W > 0);
Info.L = Info.L + order;

Kfr_d = gen_coefs(vertices,d);
size_Kfr_d = size(Kfr_d,1);

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

% A(alpha)W + WA(alpha) + BiZ(alpha) + Z(alpha)'Bi' < 0
for i=1:number_LMIs

    k = Kfr_g_d_1(i,:);

    Term = zeros(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);
                    B2 = B2i(:,(p-1)*inputs+1:p*inputs);
                    mon_string = string_monomial(k-kprime-e);
                    Z = eval(['Z',mon_string]);
                    beta = factorial(deg)*(k(p)-kprime(p))/Pi(k-kprime);

                    Term = Term + Cfr*(beta*(A*W+W'*A') + B2*Z + Z'*B2');
                end
            end
        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.W = double(W);
    Info.Ki = [];
    for i=1:number_VAR
        mon = string_monomial(K_g(i,:));
        Z = eval(['Z',mon]);
        K = double(Z)*inv(Info.W);
        Info.Ki = [Info.Ki K];
    end
    Info.feas = 1;
end
%___________________________________________________________________________________________________
function R = Pi(monomial)

R = 1;
for i=1:size(monomial,2)
    R = R*factorial(monomial(1,i));
end

%___________________________________________________________________________________________________
function [comb] = gen_coefs(vertices,grau)

nome = strcat('comb_N',int2str(vertices),'d',int2str(grau),'.mat');
if exist(nome,'file') == 2
    load(nome);
else
    for i=0:vertices
        var(i+1)=0;
    end
    comb = [];
    comb = nomial(grau, grau, vertices, 1, var,comb);
end

%___________________________________________________________________________________________________
function Coefs = nomial(n,n2,v,v2,var,Coefs)

if v2 < v
    for i=0:n2
        var(v2+1) = i;
        Coefs = nomial(n, n2-i, v, v2+1, var,Coefs);
    end
else
    var(v+1)=n2;
    prod = part(n,var,v);
    este = [];
    for j=1:v
        este = [este var(j+1)];
    end
    Coefs = [Coefs;este];
end

%___________________________________________________________________________________________________
function prod = part(n,var,v)

prod = 1;
for i=1:v
    prod = prod * factorial(var(i+1));
end
prod = factorial(n)/prod;

%___________________________________________________________________________________________________
function [string] = string_monomial(monomial)

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