function output = ssf_stab_quad_c_yal(A,B)
%function output = ssf_stab_quad_c_yal(A,B)
%
% Synthesize a robust static state feedback (ssf) stabilizing gain for a continuous-time polytopic time 
% varying system using quadratic stability relaxation. The LMI condition was published in (Geromel et al., 
% SICON, 1989). The LMIs are programmed using Yalmip and can solved using any SDP solver supported by
% Yalmip, as Sedumi.
%
%input:  A      -> dynamic matrix polytope (cell array)
%        B      -> input matrix polytope (cell array)
%
%output: output.feas     -> stabilizable (1) or not (0)
%        output.cpusec   -> cpu time in seconds to solve the LMIs
%        output.cpusec_m -> cpu time in seconds to mount the LMIs
%        output.L        -> number of scalar variables
%        output.P        -> Solution variable P
%        output.K        -> state feedback gain
%        output.L        -> number of LMIs rows
%
% Example:
%
% A{1} = [0.19 0.68;0.19 0.30];
% A{2} = [0.93 0.42;0.47 0.85];
% B{1} = [0.53;0.20];
% B{2} = [0.54;0.15];
% output = ssf_stab_quad_c_yal(A,B)
%
% Date: 26/09/2008
% Author: ricfow@dt.fee.unicamp.br

% determine the dimensions of the system
order    = size(A{1},1);
vertices = size(A,2);
inputs   = size(B{1},2);

output.cpusec_m = clock;

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

% generate the variables
W = sdpvar(order,order,'symmetric');
Z = sdpvar(inputs,order,'full');

% LMIs of quadratic stability ssf design relaxation
LMIs = LMIs + set(W > 0);
for i = 1:vertices
    LMIs = LMIs + set(A{i}*W + W*A{i}' + B{i}*Z + Z'*B{i}' < 0);
end

% determine the number of LMI rows
output.L = 0;
for i=1:size(LMIs,1)
    output.L = output.L + size(LMIs{i},1);
end
% determine the number of scalar variables
output.V = size(getvariables(LMIs),2);

% evaluate the elapsed time to mount the LMIs set
output.cpusec_m = etime(clock,output.cpusec_m);

% solve the LMIs
sol = solvesdp(LMIs,[],sdpsettings('verbose',0,'solver','sedumi'));

% evaluate the elapsed time to solve the LMIs set
output.cpusec_s = sol.solvertime;

% retrieving the minimal primal residual
p=min(checkset(LMIs));

output.feas = 0;
% capturing the solutions (if ones exist)
if p > 0
    output.feas = 1;
    output.W = double(W);
    output.K = double(Z)*inv(double(W));    
    for i=1:vertices
        output.Acl{i} = A{i}+B{i}*output.K;
    end
end
