function primal = ddstab_tvds(A,Ad,hmax,d,mysolver)
% DDSTAB_TVDS   delay-dependent stability for time-varying delay systems
%
% primal=ddstab_tvds(A,Ad,hmax,d,mysolver) provides a stability test for 
% continuous time delay systems of the form:
%
%     dx/dt = Ax(t) + Adx(t-h)   
%
% where A and Ad are nxn real constant matrices. The delay h is time 
% varying and the following constraints are assumed:
%
%    0 < h(t) < hmax    and    |dh/dt| < d
%
% The stability test is based on LMI conditions. The toolbox yalmip is
% required (tested with v3). The solver to be used (e.g. lmilab or sedumi)
% must be specified in the input argument 'mysolver' (string). The function
% outputs the primal of the LMI problem. The TVDS is stated stable if
% min(primal)>0.
%
% This code corresponds to Theorem 3 from the paper
% Y. Ariba and F. Gouaisbaut "Input-output framework for robust stability 
% of time-varying delay systems", IEEE Conference on Decision and Control 
% (CDC), 2009.
%
%
% Exemple: let define
%
%    A=[-2 0;0 -0.9]   and   Ad = [-1 0 ; -1 -1]
%
% then
%     >> ddstab_tvds(A,Ad,2.682,0.5,'sedumi');
%
%        feasible problem, the time-delay system is stable
%
% and
%    >> ddstab_tvds(A,Ad,2.683,0.5,'sedumi');
%
%        infeasible problem, stability cannot be stated
%
%
% We have proved that the tvds
%
%           [dx1/dt] =  [ -2     0  ] [x1(t)]  +   [ -1   0 ] [x1(t-h)]
%           [dx2/dt]    [  0   -0.9 ] [x2(t)]      [ -1  -1 ] [x2(t-h)]
%
% is stable forall h(t) such that
%
%      0 < h(t) < 2.682    and     |dh/dt| < 0.5
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% init
n=size(A,1);
dmax=d;
dmin=-d;

% filter
k=sqrt(8/(2-d));
b=sqrt(50);
c=sqrt(12.5);
a=sqrt(6.5+2*b);
delta=0.01;

phi=tf([k*hmax^2 k*c*hmax 0],[hmax^2 a*hmax b])+delta;
[Af Bf Cf Df]=tf2ss([k*hmax^2 k*c*hmax 0],[hmax^2 a*hmax b]);


% augmented system
E=[eye(2*n); [zeros(n) eye(n)]];
A_=[A zeros(n);zeros(n) A;eye(n) zeros(n)];
Ad_min=[Ad*(1-dmin) zeros(n);zeros(n) Ad;zeros(n) zeros(n)];
Ad_max=[Ad*(1-dmax) zeros(n);zeros(n) Ad;zeros(n) zeros(n)];


% building the interconnection
I=eye(n);
O=zeros(n,n);
I2=eye(2*n);
O2=zeros(2*n,2*n);
I4=eye(4*n);
O4=zeros(4*n,4*n);
Af2=blkdiag(Af,Af,Af,Af);
Bf2=blkdiag(Bf,Bf,Bf,Bf);
Cf2=blkdiag(Cf,Cf,Cf,Cf);
Df2=blkdiag(Df,Df,Df,Df);
EE=[E            zeros(3*n,4*n)  zeros(3*n,2*n)  zeros(3*n,2*n);
  zeros(4*n,2*n)      I4         zeros(4*n,2*n)  zeros(4*n,2*n);
    O2           zeros(2*n,4*n)      I2          zeros(2*n,2*n);
    O2           zeros(2*n,4*n)      O2                I2;
    O2           zeros(2*n,4*n)      O2                O2];

AA_min=[A_     zeros(3*n,4*n)    Ad_min               zeros(3*n,2*n);
        Bf2      Af2           zeros(4*n,2*n)         zeros(4*n,2*n);
        I2     zeros(2*n,4*n)      O2                     O2;
        Df2      Cf2           zeros(2*n,2*n)         zeros(2*n,2*n);
        I2     zeros(2*n,4*n)     -I2                    -I2];

AA_max=[A_     zeros(3*n,4*n)    Ad_max               zeros(3*n,2*n);
        Bf2      Af2           zeros(4*n,2*n)         zeros(4*n,2*n);
        I2     zeros(2*n,4*n)      O2                     O2;
        Df2      Cf2           zeros(2*n,2*n)         zeros(2*n,2*n);
        I2     zeros(2*n,4*n)     -I2                    -I2];

S_max=[EE -AA_max];
S_min=[EE -AA_min];


% separator
P=sdpvar(6*n,6*n);
Q0=sdpvar(2*n,2*n);
Q1=sdpvar(2*n,2*n);
X=sdpvar(20*n,13*n,'full');
Theta11=blkdiag(O2,O2,O2,-Q0,-Q1);
Theta12=blkdiag(-P,O2,O2);
Theta21=Theta12';
Theta22_min=blkdiag(O2,O2,O2,Q0*(1-dmin),Q1);
Theta22_max=blkdiag(O2,O2,O2,Q0*(1-dmax),Q1);
Theta_min=[Theta11,Theta12;Theta21,Theta22_min];
Theta_max=[Theta11,Theta12;Theta21,Theta22_max];


% stability condition
options=sdpsettings('solver',mysolver);

Contrainte_max=Theta_max+X*S_max+(X*S_max)';
Contrainte_min=Theta_min+X*S_min+(X*S_min)';
probleme=set(P>0)+set(Q0>0)+set(Q1>0)+set(Contrainte_min>0)+set(Contrainte_max>0);

solution=solvesdp(probleme,[],options);

primal=checkset(probleme);

if min(primal)<=0
    
    disp('---------');
    disp(' ');
    disp(' ');
    disp('infeasible problem, stability cannot be stated');
    disp(' ');
else
    disp('---------');
    disp(' ');
    disp(' ');
    disp('feasible problem, the time-delay system is stable');
    disp(' ');
end

end