function Info = synth_rob_OP08a_ACC_T3_yal(Ai,Bi,b)
% function Info = synth_rob_OP08a_ACC_T3_yal(Ai,Bi,b)
% synthesize a robust state feedback controller for a discrete-time time-varying polytopic system with
% bounded parameter variation via affine parameter-dependent Lyapunov function.
% Theorem 3 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        -> Robust controler (if exists)
%         Info.Acl      -> Closed loop dymamic matrix
%         Info.V        -> 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_rob_OP08a_ACC_T3_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;

G = sdpvar(order,order,'full');
Z = sdpvar(inputs,order,'full');
for i=1:vertices
    P{i} = sdpvar(order,order,'symmetric');
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 = [P{i} (Aii*G+Bii*Z);
            (Aii*G+Bii*Z)' G+G'-((1-b)*P{i} + b*P{m})];
        LMIs = LMIs + set(T > 0);
        Info.L = Info.L + 2*order;
    end
end

Info.cpusec_m = etime(clock,Info.cpusec_m);

Info.V = 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.K = double(Z)*inv(double(G));
    Info.Pi = [];
    Info.Acl = [];
    for i=vertices:-1:1
        Info.Pi = [Info.Pi double(P{i})];        
    end
    for i=1:vertices
        Info.Acl = [Info.Acl Ai(:,(i-1)*order+1:i*order)+Bi(:,(i-1)*inputs+1:i*inputs)*Info.K];
    end
    Info.feas = 1;
end
