function Info = synth_GS_OP08a_ACC_T4_yal(Ai,Bi,b)
% function Info = synth_GS_OP08a_ACC_T4_yal(Ai,Bi,b)
% synthesize a gain-scheduling state feedback controller for a discrete-time time-varying polytopic system with
% bounded parameter variation via affine parameter-dependent Lyapunov function.
% Theorem 4 of (Oliveira & Peres, ACC 2008)
% input:  Ai      -> vertices of the dynamic matrix (Ai = [A1 A2 ... AN])
%         Bi      -> vertices of the input matrix (Bi = [B1 B2 ... BN])    
%         b       -> parameters variation bound (0 -> time-invariant, 1 -> arbitrarily fast)
% output: 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.K        -> number of scalar variables
%         Info.L        -> number of LMI rows
%
%example: 2 states and 2 vertices
% A1 = [-0.1  1.1;   
%       -0.8  -0.3];  
% A2 = [-0.6  0.6;
%       0.4   -0.3];
% B1 = [1;0];
% B2 = [0;1];
% Ai = [A1 A2];
% Bi = [B1 B2];
% Info = synth_GS_OP08a_ACC_T4_yal(Ai,Bi,1)
% 
% 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.cpusec_m = clock;

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

for i=1:vertices
    P{i} = sdpvar(order,order,'symmetric');
    G{i} = sdpvar(order,order,'full');
    Z{i} = sdpvar(inputs,order,'full');
end

for m=1:vertices
    for i=1:vertices
        Aii = Ai(:,(i-1)*order+1:i*order);
        Bii = Bi(:,(i-1)*inputs+1:i*inputs);


        T = [((1-b)*P{i} + b*P{m}) (Aii*G{i}+Bii*Z{i});
            (Aii*G{i}+Bii*Z{i})' G{i}+G{i}'-P{i}];
        LMIs = LMIs + set(T > 0);
        Info.L = Info.L + 2*order;
        for j=i+1:vertices
            Ajj = Ai(:,(j-1)*order+1:j*order);
            Bjj = Bi(:,(j-1)*inputs+1:j*inputs);


            T = [((1-b)*P{i} + b*P{m})+((1-b)*P{j} + b*P{m}) (Aii*G{j}+Ajj*G{i}+Bii*Z{j} +Bjj*Z{i});
                (Aii*G{j}+Ajj*G{i}+Bii*Z{j} +Bjj*Z{i})' G{i}+G{i}'+G{j}+G{j}'-P{i}-P{j}];
            LMIs = LMIs + set(T > 0);
            Info.L = Info.L + 2*order;
        end
    end
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.Gi = [];
    Info.Zi = [];
    Info.Pi = [];
    for i=vertices:-1:1
        Info.Pi = [Info.Pi double(P{i})];
        Info.Gi = [Info.Gi double(G{i})];
        Info.Zi = [Info.Zi double(Z{i})];
    end
    Info.feas = 1;
end
