function [Info] = synth_h2_D_MOPB07_IET_T4_yal(Ai,B1i,B2i,Ci,D1i,D2i,deg,d,param)
%function Info = synth_h2_D_MOPB07_IET_T4_yal(Ai,B1i,B2i,Ci,D1i,D2i,deg,d,param)
% synthesize a homogeneous polynomially parameter-dependent gain that ensure closed-loop quadratic stability
% and an H-2 guaranteed cost for continuous-time arbitrarily time-varying polytopic linear systems.
% Pólya' relaxations can be used. Dual approach. This conditions was published in (Montagner et al., IET-CTA, 2007).
%
% input: Ai,B1i,B2i,D1i,D2i,Ci -> vertices of the polytopes (Ai = [A1 A2 ... AN], B2i = [B21 B22 ... B2N],...)
%        deg           -> degree of the polynomially parameter-dependent gain Z(alpha)
%        d             -> number of Pólya's relaxations
%        param         -> optional parameters
%        param.h2      -> if this parameter is different of zero, the routine will
%                                 check the feasibility of the LMIs for this value of the H-2 norm.
%output: Info.h2       -> H-2 guaranteed cost (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_k*inv(W))
%        Info.delta    -> minimal primal residual returned by the LMI solver (SeDuMi is the default).
%
% 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];
% B1 = [1;0];  B2 = [0;-1];
% Ai = [A1 A2]; B1i = [B1 B2];
% B2i = B1i; Ci = [1 0 1 0;0 1 0 1];
% D1i = zeros(2,2); D2i = [0.3 -0.2;0.2 -0.2];
% Info = synth_h2_D_MOPB07_IET_T4_yal(Ai,B1i,B2i,Ci,D1i,D2i,1,0)
%
% Date: 24/05/2007
% Author: ricfow@dt.fee.unicamp.br

if nargin == 9
    if isfield(param,'h2')
        h2 = param.h2;
    else
        h2 = 0;
    end
else
    h2 = 0;
end

%determine the number of vertices of polytope
[order,vertices] = size(Ai);
vertices = vertices / order;
[order,inputs_c] = size(B2i);%control inputs
inputs_c = inputs_c / vertices;
[order,inputs_d] = size(B1i);%disturbance inputs
inputs_d = inputs_d / vertices;
outputs = size(Ci,1);

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_c,order,''full'');']);
    eval(['M',mon_string,'=sdpvar(outputs,outputs,''symmetric'');']);
end

W = sdpvar(order,order,'symmetric');
% LMIs = LMIs + set(W > 0);
% Info.L = Info.L + order;

if h2 == 0
    eta = sdpvar(1);
    obj = eta;
else
    eta = h2^2;
    obj = [];
end

if deg > 1
    Kfr_g_d = gen_coefs(vertices,deg + d);
    Kfr_d = gen_coefs(vertices,d);
else
    Kfr_g_d = gen_coefs(vertices,deg);
    Kfr_d = gen_coefs(vertices,0);
end

number_LMIs = size(Kfr_g_d,1);
size_Kfr_d = size(Kfr_d,1);

for j=1:number_LMIs

    k = Kfr_g_d(j,:);

    Term = 0;
    for m=1:size_Kfr_d
        kprime = Kfr_d(m,:);
        if size(find(k-kprime>=0),2) == vertices
            Cfr  = factorial(d)/Pi(kprime);

            M = eval(['M',string_monomial(k-kprime)]);
            beta = factorial(deg)/Pi(k-kprime);
            Term = Term + Cfr*(trace(M) - beta*eta);
        end
    end
    LMIs = LMIs + set(Term < 0);
    Info.L = Info.L + 1;
end

if deg <= 1
    Kfr_d = gen_coefs(vertices,d);
    size_Kfr_d = size(Kfr_d,1);
end

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

for i=1:number_LMIs

    k = Kfr_g_d_1(i,:);

    Term1 = zeros(outputs+order);
    Term2 = zeros(order+inputs_d);

    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(:,inputs_c*(p-1)+1:p*inputs_c);
                    B1 = B1i(:,inputs_d*(p-1)+1:p*inputs_d);
                    C = Ci(:,order *(p-1)+1:p*order);
                    D1 = D1i(:,inputs_d*(p-1)+1:p*inputs_d);
                    D2 = D2i(:,inputs_c*(p-1)+1:p*inputs_c);

                    mon_string = string_monomial(k-kprime-e);
                    Z = eval(['Z',mon_string]);
                    M = eval(['M',mon_string]);
                    beta = factorial(deg)*(k(p)-kprime(p))/Pi(k-kprime);

                    Term1 = Term1 + Cfr*[beta*W (beta*W*C'+Z'*D2');
                        (D2*Z+C*W*beta) M];

                    Term2 = Term2 + Cfr*[beta*(A*W+W*A')+B2*Z+Z'*B2' B1*beta;
                        B1'*beta -eye(inputs_d)*beta];

                end
            end
        end
    end
    LMIs = LMIs + set(Term1 > 0);
    LMIs = LMIs + set(Term2 < 0);
    Info.L = Info.L + 2*order + outputs + inputs_d;
end
Info.cpusec_m = etime(clock,Info.cpusec_m);

Info.K = size(getvariables(LMIs),2);

sol = solvesdp(LMIs,obj,sdpsettings('verbose',0,'solver','sedumi'));
Info.cpusec = sol.solvertime;
p=min(checkset(LMIs));
Info.delta = p;
if p > 0
    p = 0;
end

Info.h2 = 0;
%capturing the solutions (if ones exist)
if abs(p) < 1e-7 %adopted precision for the minimum primal residual
    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.h2 = sqrt(double(eta));
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
