function primal = drstab_tvds(A,Ad,hmin,hmax,d,mysolver)
% DRSTAB_TVDS   delay-range stability for time-varying delay systems
%
% primal=drstab_tvds(A,Ad,hmin,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 < hmin < 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. Note that the delay-free system is note required to be 
% stable. Indeed, this stability test is able to detect some pockets of 
% stability.
%
% This code corresponds to Theorem 2 from the paper
% Y. Ariba, F. Gouaisbaut and K.H. Johansson "Stability interval for time 
% varying delay systems", IEEE Conference on Decision and Control (CDC), 
% 2010.
%
%
% Exemple: let define
%
%    A=[-2 0;0 -0.9]   and   Ad = [-1 0 ; -1 -1]
%
% then
%     >> drstab_tvds(A,Ad,0.01,4.081,0.1,'sedumi');
%
%        feasible problem, the time-delay system is stable
%
% and
%    >> drstab_tvds(A,Ad,0.01,4.082,0.1,'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.01 < h(t) < 4.081    and     |dh/dt| < 0.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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


% augmented system
I=eye(n);
O=zeros(n,n);
I2=eye(2*n);
O2=zeros(2*n,2*n);
E=[I2;O I];
Abar=[A O;O A;I O];
Adbar_min=[Ad*(1-dmin) O;O Ad;O O];
Adbar_max=[Ad*(1-dmax) O;O Ad;O O];


% building the interconnection
Smin=[A -I Ad -I*hmin O O;
       I O  -I  O    -I*hmin O;
       A  O Ad  O  -I  -I];
Smax=[A -I Ad -I*hmax O O;
       I O  -I  O    -I*hmax O;
       A  O Ad  O  -I  -I];
Nmin=[A*A  Ad*(1-dmin)  A*Ad  O  O  O;
       A        O        Ad   O  O  O;
       A        O        Ad   O  O  O;
       I        O        O    O  O  O;
       A*A  Ad*(1-dmin)  A*Ad  O  O  O;
       A        O        Ad   O  O  O;
       A*A  Ad*(1-dmin)  A*Ad  O  O  O;
       A        O        Ad   O  O  O;
       I        O        O    O  O  O;
       O        I        O    O  O  O;
       O        O        I    O  O  O;
       O        O        O    I  O  O;
       O        O        O    O  I  O;
       O        O        O    O  O  I];
Nmax=[A*A  Ad*(1-dmax)  A*Ad  O  O  O;
       A        O        Ad   O  O  O;
       A        O        Ad   O  O  O;
       I        O        O    O  O  O;
      A*A  Ad*(1-dmax)  A*Ad  O  O  O;
       A        O        Ad   O  O  O;
      A*A  Ad*(1-dmax)  A*Ad  O  O  O;
       A        O        Ad   O  O  O;
       I        O        O    O  O  O;
       O        I        O    O  O  O;
       O        O        I    O  O  O;
       O        O        O    I  O  O;
       O        O        O    O  I  O;
       O        O        O    O  O  I];

   
% separator
P=sdpvar(2*n,2*n);
Q0=sdpvar(2*n,2*n);
Q1=sdpvar(2*n,2*n);
Q2=sdpvar(n,n);
X=sdpvar(6*n,3*n,'full');
Theta11=blkdiag(O2,-Q0,-Q1*hm,-Q2*hm^2/2);
Theta12=blkdiag(-P,O2,O2,O);
Theta21=Theta12';
Theta22maxmin=blkdiag(O2,Q0*(1-dmax),Q1*hmin,2*Q2);
Theta22minmax=blkdiag(O2,Q0*(1-dmin),Q1*hmax,2*Q2);
Theta22maxmax=blkdiag(O2,Q0*(1-dmax),Q1*hmax,2*Q2);
Theta22minmin=blkdiag(O2,Q0*(1-dmin),Q1*hmin,2*Q2);

Thetaminmin=[Theta11,Theta12;Theta21,Theta22minmin];
Thetamaxmin=[Theta11,Theta12;Theta21,Theta22maxmin];
Thetaminmax=[Theta11,Theta12;Theta21,Theta22minmax];
Thetamaxmax=[Theta11,Theta12;Theta21,Theta22maxmax];


% stability condition
options=sdpsettings('solver',mysolver);
   
Contrainteminmin=Nmin'*Thetaminmin*Nmin+X*Smin+(X*Smin)';
Contraintemaxmin=Nmax'*Thetamaxmin*Nmax+X*Smin+(X*Smin)';
Contrainteminmax=Nmin'*Thetaminmax*Nmin+X*Smax+(X*Smax)';
Contraintemaxmax=Nmax'*Thetamaxmax*Nmax+X*Smax+(X*Smax)';

probleme=set(P>0)+set(Q0>0)+set(Q1>0)+set(Q2>0)+set(Contraintemaxmin>0)+set(Contrainteminmax>0)+set(Contraintemaxmax>0)+set(Contrainteminmin>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

