function Info = h2_MJLS_OVdVP09_IJC_T2_yal(A,E,C,G,mu,deg,d,param)
% function Info = h2_MJLS_OVdVP09_IJC_T2_yal(A,E,C,G,mu,deg,d,param)
% compute a H-2 guarantee cost of a discrete-time markovian jump linear
% system with polytopic uncertainties in the probability matrix.
% A homogeneous polynomially parameter-dependent Lyapunov matrix of arbitray degree
% and Polya's relaxations are used in the PRIMAL approach. The condition is given by Theorem 2 in
% (Oliveira, Vargas, do Val & Peres, IJC, 2009)
%
% inputs: A        -> dynamic matrix vertices: A{1} ... A{sigma}
%         E        -> disturbance matrix vertices: E{1} ... E{sigma}
%         C        -> output matrix vertices: C{1} ... C{sigma}
%         G        -> vertices of the probability matrix: G{1} ... G{v}
%         deg      -> degree of the polynomially parameter dependent Lyapunov matrix
%         d        -> number of Polya'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
%                                 h2 norm.
%
% outputs: 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.Pki      -> Solution variables Pk{i}, i=1,...,sigma
%          Info.K        -> number of scalar variables used in the optmization problem
%          Info.L        -> number of LMI rows used in the optmization problem
%          Info.delta    -> minimal primal residual returned by the LMI solver (SeDuMi is the default).
%
%
% Example: scalar system, two states of Markov, two vertices in the probability matrix
%
% A{1} = sqrt(0.9); A{2} = sqrt(1.1);
% G{1} = [0.9 0.1;0.3 0.7]; G{2} = [0.95 0.05;0.15 0.85];
% E{1} = 0.5; E{2} = 0.75; C{1} = 1; C{2} = 0.75;
% mu = [1 0];
% deg = 2;
% d = 0;
% Info = h2_MJLS_OVdVP09_IJC_T2_yal(A,E,C,G,mu,deg,d)
% 
% Date: 29/05/2007
% Author: ricfow@dt.fee.unicamp.br

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

sigma = size(A,2);      % number of markovian states
order = size(A{1},1);   % number of states of the system
vertices = size(G,2);   % number of vertices of the probability matrix
inputs = size(E{1},2);  % number of disturbance inputs
outputs = size(C{1},1); % number of outputs
Gamma = mount_Gamma(G,vertices,sigma,order);

Info.cpusec_m = clock;

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

%generate the exponents of the monomials to create the variables P_k
K_g = gen_coefs(vertices,deg);
number_VAR = size(K_g,1);%number of variables

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

for r=1:number_VAR
    mon_string = string_monomial(K_g(r,:));
    Y = [];
    for i=1:sigma
        eval(['P{i}.P',mon_string,' = sdpvar(order);']);
        eval(['Y = [Y;A{i}*P{i}.P',mon_string,'*A{i}''];']);
    end
    eval(['Yfr_',mon_string,' = Y;']);
end

Zfr = [];
for i=1:sigma
    Zfr = [Zfr;mu(i)*E{i}*E{i}'];
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);

%Lyapunov matrix positive definite: P_i(alpha) > 0, i=1,...,sigma
for i=1:sigma
    for j=1:number_LMIs

        k = Kfr_g_d(j,:);

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

                P_ = eval(['P{i}.P',string_monomial(k-kprime)]);
                Term_1 = Term_1 + Cfr*P_;
            end
        end
        LMIs = LMIs + set(Term_1 > 0);
        Info.L = Info.L + order;
    end
end

%Bounding the cost : rho > sum_j trace (Cj P_j(alpha) Cj')
for j=1:number_LMIs

    k = Kfr_g_d(j,:);

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

            mon_string = string_monomial(k-kprime);            
            
            sum_P = 0;
            for i=1:sigma
                P_ = eval(['P{i}.P',mon_string]);
                sum_P = sum_P + C{i}*(P_)*C{i}';
            end            

            Term_2 = Term_2 + Cfr*(trace(sum_P) - rho*(factorial(deg)/Pi(k-kprime)));
        end
    end
    LMIs = LMIs + set(Term_2 < 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);%number of LMIs to be generated

%Grammian condition: sum pij (Aj Pj(alpha) Aj' + mu EjEj') - Pi(alpha) < 0
for i=1:sigma
    for j=1:number_LMIs

        k = Kfr_g_d_1(j,:);

        Term_1 = zeros(order);
        for m=1:size_Kfr_d
            kprime = Kfr_d(m,:);
            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;

                        Gamma_i = Gamma{i,p};
                        Beta = factorial(deg)*(k(p)-kprime(p))/Pi(k-kprime);

                        mon_string = string_monomial(k-kprime-e);
                        P_ = eval(['P{i}.P',mon_string]);
                        Yf = eval(['Yfr_',mon_string]);

                        Term_1 = Term_1 + Cfr*(Gamma_i'*(Yf + Beta*Zfr) - P_);
                    end
                end
            end
        end
        LMIs = LMIs + set(Term_1 < 0);
        Info.L = Info.L + order;
    end
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
    for i=1:sigma
        Info.Pki{i}   = [];
        for j=1:number_VAR
            mon_string = string_monomial(K_g(j,:));
            Info.Pki{i} = [Info.Pki{i} double(eval(['P{i}.P',mon_string]))];
        end
    end
    Info.h2 = sqrt(double(rho));
end

%___________________________________________________________________________________________________
function R = Pi(monomial)

R=prod(factorial(monomial));

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

%___________________________________________________________________________________________________
function Gamma = mount_Gamma(G,v,sigma,n)

for i=1:sigma
    for r=1:v
        Vi = [];
        for j=1:sigma
            Vi = [Vi;eye(n)*G{r}(j,i)];
        end
        Gamma{i,r} = Vi;
    end
end
