%   Script VERADemo.M
%   
%   Demonstrates the event-triggered control for VERA
%   Data from automatica and Control Engineering Practice
%
%   Figs. 2, 3, 4 are not tranferrable into eps
%
% J. Lunze  
% 29.10.2018
% version of 14.1.2019
% for 2nd edition: 1.5.2021
% 
echo off
clear
close all

%%    Simplified VERA model
%    x_1 - level            operation point: 40 cm
%    x_2 - temperature      operation point: 313 K
%    disturbance = cold water inflow [0, 1] (valve angle)
%    u_1 - inflow into TB  [0, 1] (valve angle)
%    u_2 - heating power  [0, 6]*1000 W
A=[-0.001 0; -0.00002 -0.002];
B=[0.02 0; -1 0.15];
C=eye(2,2);
E=[0.04; -0.16];
sigma=[-0.04, -0.03];
K=place(A, B, sigma);
%K=[5 -0.1; 10 2];
VERA=ss(A, E, C, 0);
VERAu=ss(A, B, C, 0);
controlledVERA=ss(A-B*K, E, C, 0);
x1bar=40;
x2bar=313;
u1bar=0.36;
u2bar=1.5;
%   discrete-time VERA
Tend=800;
Ts=0.5;
Time=0:Ts:Tend;
kend=length(Time);
VERAd=c2d(VERA, Ts);   % disturbance input
VERAud=c2d(VERAu, Ts); % control input
controlledVERAd=c2d(controlledVERA, Ts);
Ad=VERAd.A;
Ed=VERAd.B;
Bd=VERAud.B;
Cd=C;
Abar=controlledVERAd.A;
Ebar=controlledVERAd.B;
x0=zeros(2,1);
%  disturbance
t1=100;
kt1=ceil(t1/Ts);
t2=300;
kt2=ceil(t2/Ts);
t3=450;
kt3=ceil(t3/Ts);
kd=0.25/300;
Dist=zeros(kend,1);
for k1=kt1:kt2
    Dist(k1,1)=0.3;
end
for k1=kt3:kend
    Dist(k1,1)=kd*(k1-kt3)*Ts;
end

%%   Disturbance behaviour of the uncontrolled and controlled process
%
figure(1)
Y=lsim(controlledVERA,Dist,Time,zeros(2,1));
Yopen=lsim(VERA,Dist,Time,zeros(2,1));
%
subplot(6,1,[1 2])
plot(Time, Dist, 'b');
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([0 0.3]);
axis([0 Tend -0.02 0.4]);
latextitle('$${\tt VERADemo:}$$ Comparison of uncontrolled and controlled process');
%
subplot(6,1,[3 4])
plot(Time, Y(:,1), 'b', Time, Yopen(:,1), 'b--');
latexylabel('$$x_1$$ in cm');
rotateY;
xleer;
yticks([0 3]);
axis([0 Tend -0.1 3.5]);
%
subplot(6,1,[5 6])
plot(Time, Y(:,2), 'b', Time, Yopen(:,2), 'b--');
latexylabel('$$x_2$$ in K');
rotateY;
latexxlabel('$$t$$ in s');
yticks([-10 -5 0]);
xticks([0 200 400 600 800 1000]);
axis([0 Tend -12 2]);


%%  Event-triggered control
%   event threshold: ebar
figure(2);
ebar=0.4;
%   Determination of emax
%     vector norm = max |x_i|
%     induced matrix norm = max_i  sum_j |g_ij|
Errorsys=ss(A-B*K, B*K, eye(2,2), 0);
[G,TimeG]=impulse(Errorsys);
for k1=1:length(TimeG)
    Gnorm(k1,1)=max(sum(abs([G(k1,1,1) G(k1,1,2); G(k1,2,1) G(k1,2,2)])'));
end
TsG=TimeG(2)-TimeG(1);
emax=ebar*trapz(Gnorm)*TsG;
%   continuous feedback loop with tolerance band
subplot(9,1,[3,4,5]);
fill([Time flip(Time) 0], [Y(:,1)'+emax flip(Y(:,1)'-emax) Y(1,1)+emax], [0.9,0.9,0.9]);
hold on
subplot(9,1,[6,7,8]);
fill([Time flip(Time) 0], [Y(:,2)'+emax flip(Y(:,2)'-emax) Y(1,2)+emax],  [0.9,0.9,0.9]);
hold on
%
subplot(9,1,[1,2]);
plot(Time, Dist, 'b');
hold on
%
X=[];
Xs=[];
Dapprox=[];
X(1,:)=x0';
Xs(1,:)=x0';
Dapprox(1,1)=0;
told=0;
Eventtime=zeros(kend,1);
Eventtime(1,1)=1;
for k1=2:kend
    Dapprox(k1,1)=Dapprox(k1-1,1);
    Test(k1)=max([abs(X(k1-1,1)-Xs(k1-1,1)), abs(X(k1-1,2)-Xs(k1-1,2))]);
   if Test(k1)>ebar
       tnew=(k1-1)*Ts;   % new event time
       H=inv(A)*(expm(A*(tnew-told))-eye(2,2))*E;
       Hinv=inv(H'*H)*H';
       Dapprox(k1,1)=Dapprox(k1-1,1)+Hinv*(X(k1-1,:)-Xs(k1-1,:))';
       %Dapprox(k1,1)=0;
       told=tnew;
         % reset the state of the input generator
       Xs(k1,:)=Abar*X(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*X(k1-1,:)';
       if abs(abs(X(k1-1,1)-Xs(k1-1,1)))>ebar
          Eventtime(k1-1,1)=1;   % event generated by x_1
       else
          Eventtime(k1-1,1)=2;   % event generated by x_2
       end
   else
       Xs(k1,:)=Abar*Xs(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*Xs(k1-1,:)';
   end
   X(k1,:)=Ad*X(k1-1,:)' +Bd*U(k1-1,:)' +Ed*Dist(k1-1,1);
end
subplot(9,1,[1,2]);
plot(Time, Dapprox, 'b--');
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([0 0.3]);
axis([0 Tend -0.02 0.4]);
latextitle('$${\tt VERADemo:}$$ Event-triggered control');
hold off
subplot(9,1,[3,4,5])
plot(Time, X(:,1),'b', Time, Xs(:,1), 'b--');
latexylabel('$$x_1$$ in cm');
rotateY;
xleer;
yticks([-1 0 1]);
axis([0 Tend -1 1.2]);
hold off
subplot(9,1,[6,7,8])
plot(Time, X(:,2), 'b', Time, Xs(:,2), 'b--');
latexylabel('$$x_2$$ in K$$\;$$');
rotateY;
xleer;
yticks([-2 -1 0]);
axis([0 Tend -2.5 1]);
hold off
subplot(9,1,9)
for k1=1:kend
    if Eventtime(k1)==2
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
    if Eventtime(k1)==1
        plot(k1*Ts, 1, 'ro',...
            'MarkerEdgeColor','r',...
            'MarkerFaceColor','r');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'r');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off

%%   Periodic sampling with the same tolerance bands
%
figure(3)
Tsperiod=15;     % sampling time
%   discrete-time VERA
Timeperiod=0:Tsperiod:Tend;
kendperiod=length(Timeperiod);
VERAperiod=c2d(VERA, Tsperiod);   % disturbance input
VERAuperiod=c2d(VERAu, Tsperiod); % control input
Aperiod=VERAperiod.A;
Eperiod=VERAperiod.B;
Bperiod=VERAuperiod.B;
Cd=C;
x0=zeros(2,1);
%  disturbance
kt1=ceil(t1/Tsperiod);
kt2=ceil(t2/Tsperiod);
kt3=ceil(t3/Tsperiod);
Distperiod=zeros(kendperiod,1);
for k1=kt1:kt2
    Distperiod(k1,1)=0.3;
end
for k1=kt3:kendperiod
    Distperiod(k1,1)=kd*(k1-kt3)*Tsperiod;
end
subplot(9,1,[1,2]);
plot(Timeperiod, Distperiod, 'b');
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([0 0.3]);
axis([0 Tend -0.02 0.4]);
latextitle('$${\tt VERADemo:}$$ Sampled-data control');
hold off
%
X=[];
X(1,:)=x0';
Eventtime=zeros(kend,1);
Eventtime(1,1)=1;
for k1=2:kendperiod
    U(k1-1,:)=-K*X(k1-1,:)';
    X(k1,:)=Aperiod*X(k1-1,:)' +Bperiod*U(k1-1,:)' +Eperiod*Distperiod(k1-1,1);
    Eventtime(k1,1)=1;
end
subplot(9,1,[3,4,5]);
fill([Time flip(Time) 0], [Y(:,1)'+emax flip(Y(:,1)'-emax) Y(1,1)+emax],  [0.9,0.9,0.9]);
hold on
plot(Timeperiod, X(:,1), 'b');
latexylabel('$$x_1$$ in cm');
rotateY;
xleer;
yticks([-1 0 1]);
axis([0 Tend -1 1.2]);
hold off
subplot(9,1,[6,7,8]);
fill([Time flip(Time) 0], [Y(:,2)'+emax flip(Y(:,2)'-emax) Y(1,2)+emax],  [0.9,0.9,0.9]);
hold on
plot(Timeperiod, X(:,2),'b');
latexylabel('$$x_2$$ in K$$\;$$');
rotateY;
xleer;
yticks([-2 -1 0]);
axis([0 Tend -2.5 1]);
hold off
subplot(9,1,9)
for k1=1:kendperiod
    if Eventtime(k1)==1
        plot(k1*Tsperiod, 1, 'bo');
        hold on
        plot([k1*Tsperiod k1*Tsperiod],[1 0], 'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off


%%  Event-triggered control with scaled vector norm
%    w_1*x_1 or x_2 should satisfy the event condition
%    weighting matrices W_y = diag(wy1, 1)
%                       W_x = diag(wx1, 1)
figure(4);
wy1=4;
wx1=3.9;
%   Determination of emax
%     vector norm = max |Wx*x|
%     induced matrix norm = max_i  sum_j |Wy*G*Wx^-1|
Wy=diag([wy1 1]);
Wx=diag([wx1 1]);
Errorsys=ss(A-B*K, B*K, eye(2,2), 0);
[G,TimeG]=impulse(Errorsys);
for k1=1:length(TimeG)
    Gtilde=Wy*[G(k1,1,1) G(k1,1,2); G(k1,2,1) G(k1,2,2)]*inv(Wx);
    Gnorm(k1,1)=max(sum(abs(Gtilde)'));
end
TsG=TimeG(2)-TimeG(1);
emaxW=ebar*trapz(Gnorm)*TsG;
%   continuous feedback loop with tolerance band
subplot(9,1,[3,4,5]);
fill([Time flip(Time) 0], [Y(:,1)'+emaxW/wy1 flip(Y(:,1)'-emaxW/wy1) Y(1,1)+emaxW/wy1],  [0.9,0.9,0.9]);
hold on
subplot(9,1,[6,7,8]);
fill([Time flip(Time) 0], [Y(:,2)'+emaxW flip(Y(:,2)'-emaxW) Y(1,2)+emaxW],  [0.9,0.9,0.9]);
hold on
%
subplot(9,1,[1,2]);
plot(Time, Dist, 'b');
hold on
%
X=[];
Xs=[];
Dapprox=[];
X(1,:)=x0';
Xs(1,:)=x0';
Dapprox(1,1)=0;
told=0;
Eventtime=zeros(kend,1);
Eventtime(1,1)=1;
for k1=2:kend
    Dapprox(k1,1)=Dapprox(k1-1,1);
    Test(k1)=max([wx1*abs(X(k1-1,1)-Xs(k1-1,1)), abs(X(k1-1,2)-Xs(k1-1,2))]);
   if Test(k1)>ebar
       tnew=(k1-1)*Ts;   % new event time
       H=inv(A)*(expm(A*(tnew-told))-eye(2,2))*E;
       Hinv=inv(H'*H)*H';
       Dapprox(k1,1)=Dapprox(k1-1,1)+Hinv*(X(k1-1,:)-Xs(k1-1,:))';
       %Dapprox(k1,1)=Dapprox(k1-1,1);
       told=tnew;
         % reset the state of the input generator
       Xs(k1,:)=Abar*X(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*X(k1-1,:)';
       if abs(wx1*abs(X(k1-1,1)-Xs(k1-1,1)))>ebar
          Eventtime(k1-1,1)=1;   % event generated by x_1
       else
          Eventtime(k1-1,1)=2;   % event generated by x_2
       end
   else
       Xs(k1,:)=Abar*Xs(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*Xs(k1-1,:)';
   end
   X(k1,:)=Ad*X(k1-1,:)' +Bd*U(k1-1,:)' +Ed*Dist(k1-1,1);
end
subplot(9,1,[1,2]);
plot(Time, Dapprox, 'b--');
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([0 0.3]);
axis([0 Tend -0.02 0.4]);
latextitle('$${\tt VERADemo:}$$ Event-triggered control');
hold off
subplot(9,1,[3,4,5])
plot(Time, X(:,1),'b', Time, Xs(:,1), 'b--');
latexylabel('$$x_1$$ in cm');
rotateY;
xleer;
yticks([-1 0 1]);
axis([0 Tend -1 1.2]);
hold off
subplot(9,1,[6,7,8])
plot(Time, X(:,2), 'b', Time, Xs(:,2), 'b--');
latexylabel('$$x_2$$ in K$$\;$$');
rotateY;
xleer;
yticks([-2 -1 0]);
axis([0 Tend -2.5 1]);
hold off
subplot(9,1,9)
for k1=1:kend
    if Eventtime(k1)==2
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
    if Eventtime(k1)==1
        plot(k1*Ts, 1, 'ro',...
            'MarkerEdgeColor','r',...
            'MarkerFaceColor','r');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'r');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off


%%  Practical stability 
%   Controller as in figure(4)
%                       
figure(5);    
Tendlong=1200;
Timelong=0:Ts:Tendlong;
kendlong=length(Timelong);
%  determination of the ultimate bound
emaxI=max(sum(abs((Wy*inv(A-B*K)*B*K*inv(Wx))')))*ebar;  % ultimate bound
%  determination of the disturbance bound dbar
Plantdist=ss(A, E, eye(2,2), 0);
[Gd,TimeGd]=impulse(Plantdist);
for k1=1:length(TimeGd)
    Gtilde=Wy*[Gd(k1,1); Gd(k1,2)];
    Gnorm(k1,1)=max(abs(Gtilde));
end
TsG=TimeGd(2)-TimeGd(1);
dbar=ebar/trapz(Gnorm)/TsG;
%  determination of several trajectories ending in the set A
%  new disturbance: Dist2
t1=30;
kt1=ceil(t1/Ts);
t2=230;
kt2=ceil(t2/Ts);
Dist2=zeros(kendlong,1);
for k1=kt1:kt2
    Dist2(k1,1)=0.02;  % 0.01 -> 0.006
end
%
x0=[0.4; 1];
X=[];
Xs=[];
Dapprox=[];
X(1,:)=x0';
Xs(1,:)=x0';
Dapprox(1,1)=0;
told=0;
Eventtime=zeros(kendlong,1);
Eventtime(1,1)=1;
for k1=2:kendlong
    Dapprox(k1,1)=Dapprox(k1-1,1);
    Test(k1)=max([wx1*abs(X(k1-1,1)-Xs(k1-1,1)), abs(X(k1-1,2)-Xs(k1-1,2))]);
   if Test(k1)>ebar
       tnew=(k1-1)*Ts;   % new event time
       H=inv(A)*(expm(A*(tnew-told))-eye(2,2))*E;
       Hinv=inv(H'*H)*H';
       Dapprox(k1,1)=Dapprox(k1-1,1)+Hinv*(X(k1-1,:)-Xs(k1-1,:))';
       %Dapprox(k1,1)=Dapprox(k1-1,1);
       told=tnew;
         % reset the state of the input generator
       Xs(k1,:)=Abar*X(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*X(k1-1,:)';
       if abs(wx1*abs(X(k1-1,1)-Xs(k1-1,1)))>ebar
          Eventtime(k1-1,1)=1;   % event generated by x_1
       else
          Eventtime(k1-1,1)=2;   % event generated by x_2
       end
   else
       Xs(k1,:)=Abar*Xs(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*Xs(k1-1,:)';
   end
   X(k1,:)=Ad*X(k1-1,:)' +Bd*U(k1-1,:)' +Ed*Dist2(k1-1,1);
end
figure(5)
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(end,1), X(end,2), 'b*');
plot([-emaxI/wy1 emaxI/wy1 emaxI/wy1 -emaxI/wy1 -emaxI/wy1],...
     [-emaxI -emaxI emaxI emaxI -emaxI], 'k--');
plot(0, 0, 'b+');
%plot(Xs(:,1), Xs(:,2), 'b--');
axis([-0.2 0.2 -0.5 0.5]);
latexylabel('$$x_2$$ in K');
rotateY;
latexxlabel('$$x_1$$ in cm');
latextitle('$${\tt VERADemo:}$$ Practical stability');
yticks([-0.5 -0.25 0 0.25 0.5]);
xticks([-0.2 -0.1 0 0.1 0.2]);

Dapprox(end)

figure(6)
subplot(9,1,[1,2]);
plot(Timelong, Dist2, 'b', Timelong, Dapprox, 'b--');
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([-0.01 0 0.02]);
axis([0 Tendlong -0.003 0.025]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(9,1,[3,4,5])
plot(Timelong, X(:,1),'b', Timelong, Xs(:,1), 'b--');
latexylabel('$$x_1$$');  %  in cm
rotateY;
xleer;
yticks([ 0 0.3]);
axis([0 Tendlong -0.1 0.4]);
hold off
%
subplot(9,1,[6,7,8])
plot(Timelong, X(:,2), 'b', Timelong, Xs(:,2), 'b--');
latexylabel('$$x_2\;$$');  % in K
rotateY;
xleer;
yticks([ 0 1]);
axis([0 Tendlong -1 1]);
hold off
%
subplot(9,1,9)
for k1=1:kendlong
    if Eventtime(k1)==1
        plot(k1*Ts, 1, 'ro',...
            'MarkerEdgeColor','r',...
            'MarkerFaceColor','r');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'r');
    end
    if Eventtime(k1)==2
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tendlong 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000 1200]);
hold off


%%  Practical stability 
%   Draw several trajectories ending in different equilibria
%                       
figure(7);    
%
for kit=1:4
   if kit==1
      x0=[-1;-0.6];
      dk=-0.003;
    end
    if kit==2
        x0=[-1;-0.6];
        dk=0.002;
    end
    if kit==3
        x0=[-0.2;0.4];
        dk=-0.001;
    end
    if kit==4
        x0=[-1; -1];
        dk=0;
    end
X=[];
Xs=[];
X(1,:)=x0';
Xs(1,:)=x0';
Dapprox(1,1)=dk;
told=0;
Eventtime=zeros(kendlong,1);
Eventtime(1,1)=1;
for k1=2:kendlong
    Dapprox(k1,1)=Dapprox(k1-1,1);
    Test(k1)=max([wx1*abs(X(k1-1,1)-Xs(k1-1,1)), abs(X(k1-1,2)-Xs(k1-1,2))]);
   if Test(k1)>ebar
       tnew=(k1-1)*Ts;   % new event time
       H=inv(A)*(expm(A*(tnew-told))-eye(2,2))*E;
       Hinv=inv(H'*H)*H';
       Dapprox(k1,1)=Dapprox(k1-1,1)+Hinv*(X(k1-1,:)-Xs(k1-1,:))';
       %Dapprox(k1,1)=Dapprox(k1-1,1);
       told=tnew;
         % reset the state of the input generator
       Xs(k1,:)=Abar*X(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*X(k1-1,:)';
       if abs(wx1*abs(X(k1-1,1)-Xs(k1-1,1)))>ebar
          Eventtime(k1-1,1)=1;   % event generated by x_1
       else
          Eventtime(k1-1,1)=2;   % event generated by x_2
       end
   else
       Xs(k1,:)=Abar*Xs(k1-1,:)'+Ebar*Dapprox(k1-1,1);
       U(k1-1,:)=-K*Xs(k1-1,:)';
   end
   X(k1,:)=Ad*X(k1-1,:)' +Bd*U(k1-1,:)' +Ed*Dist2(k1-1,1);
end
plot(X(:,1), X(:,2), 'b');
hold on
plot(X(end,1), X(end,2), 'b*');
plot([-emaxI/wy1 emaxI/wy1 emaxI/wy1 -emaxI/wy1 -emaxI/wy1],...
     [-emaxI -emaxI emaxI emaxI -emaxI], 'k--');
plot(0, 0, 'b+');
end
axis([-0.2 0.2 -0.5 0.5]);
latexylabel('$$x_2$$ in K');
rotateY;
latexxlabel('$$x_1$$ in cm');
latextitle('$${\tt VERADemo:}$$ Practical stability');
yticks([-0.5 -0.25 0 0.25 0.5]);
xticks([-0.2 -0.1 0 0.1 0.2]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Disturbance behaviour of the decentralised control system
%
figure(8)
K1=2;
K2=6;
Kdec=diag([K1, K2]);
deccontrolledVERA=ss(A-B*Kdec, E, C, 0);
Ydec=lsim(deccontrolledVERA,Dist,Time,zeros(2,1));
Yopen=lsim(VERA,Dist,Time,zeros(2,1));
%
subplot(6,1,[1 2])
plot(Time, Dist, 'b');
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([0 0.5]);
axis([0 Tend -0.02 0.6]);
latextitle('$${\tt VERADemo:}$$ Uncontrolled and decentralised controlled process');
%
subplot(6,1,[3 4])
plot(Time, Ydec(:,1), 'b', Time, Yopen(:,1), 'b--');
latexylabel('$$x_1$$ in cm');
rotateY;
xleer;
yticks([0 3]);
axis([0 Tend -0.1 3.5]);
%
subplot(6,1,[5 6])
plot(Time, Ydec(:,2), 'b', Time, Yopen(:,2), 'b--');
latexylabel('$$x_2$$ in K');
rotateY;
latexxlabel('$$t$$ in s');
yticks([-10 -5 0]);
xticks([0 200 400 600 800 1000]);
axis([0 Tend -12 2]);


%%  Event-triggered decentralised control 
%    use the same weighting as in figure(4)
%    no disturbance estimation
%                      
figure(9);
%  decomposed discrete-time plant model
Controlloop1=ss(A(1,1)-B(1,1)*K1, E(1,1), 1, 0);
Controlloop1d=c2d(Controlloop1, Ts);
A1=Controlloop1d.A;
E1=Controlloop1d.B;
C1=Cd(1,1);
Controlloop2=ss(A(2,2)-B(2,2)*K2, E(2,1), 1, 0);
Controlloop2d=c2d(Controlloop2, Ts);
A2=Controlloop2d.A;
E2=Controlloop2d.B;
C2=Cd(1,1);
wx1=1.2;
wx2=0.1;
ebar1=ebar/wx1;
ebar2=ebar/wx2;
%  new disturbance
t1=100;
kt1=ceil(t1/Ts);
t2=300;
kt2=ceil(t2/Ts);
t3=450;
kt3=ceil(t3/Ts);
kd=0.2/300;
Distdec=zeros(kend,1);
for k1=kt1:kt2
    Distdec(k1,1)=0.1;
end
for k1=kt3:kend
    Distdec(k1,1)=kd*(k1-kt3)*Ts;
end
%
subplot(10,1,[1,2]);
plot(Time, Distdec, 'b');
hold on
%
X=[];
Xs1=[];
Xs2=[];
X(1,:)=x0';
Xs1(1)=x0(1,1);
Xs2(1)=x0(2,1);
Eventtime1=zeros(kend,1);
Eventtime2=zeros(kend,1);
Eventtime1(1,1)=1;
Eventtime2(1,1)=1;
for k1=2:kend
   %  feedback loop 1
   if abs(X(k1-1,1)-Xs1(k1-1))>ebar1
         % reset the state of the input generator
       Xs1(k1)=A1*X(k1-1,1);
       U1(k1-1)=-K1*X(k1-1,1);
       Eventtime1(k1-1,1)=1;   % event generated by x_1
    else
       Xs1(k1)=A1*Xs1(k1-1);
       U1(k1-1)=-K1*Xs1(k1-1);
   end
   %  feedback loop 2
   if abs(X(k1-1,2)-Xs2(k1-1))>ebar2
         % reset the state of the input generator
       Xs2(k1)=A2*X(k1-1,2);
       U2(k1-1)=-K2*X(k1-1,2);
       Eventtime2(k1-1,1)=1;   % event generated by x_1
    else
       Xs2(k1)=A2*Xs2(k1-1);
       U2(k1-1)=-K2*Xs2(k1-1);
   end
   X(k1,:)=Ad*X(k1-1,:)' +Bd*[U1(k1-1);U2(k1-1)] +Ed*Distdec(k1-1,1);
end
%
subplot(10,1,[1,2]);
latexylabel('$$d\;$$');
rotateY;
xleer;
yticks([0 0.3]);
axis([0 Tend -0.04 0.3]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(10,1,[3,4,5])
plot(Time, X(:,1),'b', Time, Xs1, 'b--');
latexylabel('$$x_1$$ in cm');
rotateY;
xleer;
yticks([-0.5 0 0.5]);
axis([0 Tend -0.2 0.7]);
hold off
%
subplot(10,1,6)
for k1=1:kend
    if Eventtime1(k1)==1
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
xleer;
hold off
%
subplot(10,1,[7,8,9])
plot(Time, X(:,2), 'b', Time, Xs2, 'b--');
latexylabel('$$x_2$$ in K$$\;$$');
rotateY;
xleer;
yticks([0 2 4]);
axis([0 Tend -2 5]);
hold off
%
subplot(10,1,10)
for k1=1:kend
    if Eventtime2(k1)==1
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Linearised model of the thermo-fluid process
%  from Stöcker 2014, pp. 188-189
%
%  Two reactors with two state variables each: [level; temperature]
%  two inputs at each reactor: u_11=u_T1, u_12=U_CU
%                              u_21=u_T3, u_22=u_H
%  linearised around the equilibrium l_1=0.33 m, theta_1=294.7 K
%                                    l_2=0.34 m, theta_2=300.2K
%                                    u_11=0.5, u_12=0.5
%                                    u_21=0.5, u_22=0.5
%  valves in the coupling pipes: u_BS=0.19, u_SB=0.22, u_BW=0.21, u_SW=0.29
%  
%  coupling strength alpha1, alpha2 simulates stronger/weaker interactions
alpha1=1.5;
alpha2=1.5;
A1=0.001*[-5.74 0; -34.5 -8.58];
B1=0.001*[2.3 0; 0 -38.9];
E1=0.001*[0; 169];
Es1=alpha1*0.001*[2.42 0; 43.9 5.44];
Cz1=eye(2,2);
C1=eye(2,2);    % outputs: y_11=x_11, y_12=x_12
SS1=ss(A1, B1, C1, 0);
SS1d=c2d(SS1, Ts);
A1d=SS1d.A;
B1d=SS1d.B;
SS1dist=ss(A1, Es1, C1, 0);
SS1distd=c2d(SS1dist, Ts);
Es1d=SS1distd.B;
%
A2=0.001*[-5.00 0; -39.2 -5.58];
B2=0.001*[2.59 0; 0 35.0];
E2=0.001*[1.16; -20.7];
Es2=alpha2*0.001*[2.85 0; -46.5 5.58];
Cz2=eye(2,2);
C2=eye(2,2);    % outputs: y_21=x_21, y_22=x_22
SS2=ss(A2, B2, C2, 0);
SS2d=c2d(SS2, Ts);
A2d=SS2d.A;
B2d=SS2d.B;
SS2dist=ss(A2, Es2, C2, 0);
SS2distd=c2d(SS2dist, Ts);
Es2d=SS2distd.B;
%  initial state
x0=[2; 2; 3; -2];
%
%  Overall system without controller
Null=zeros(2,2);
AGS=[A1 Es1*Cz2; Es2*Cz1 A2];
BGS=[B1 Null; Null B2];
EGS=[E1 zeros(2,1); zeros(2,1) E2];
CGS=[C1 Null; Null C2];
GS=ss(AGS, BGS, CGS, 0);
GSd=c2d(GS, Ts);
AGSd=GSd.A;
BGSd=GSd.B;
GSdist=ss(AGS, EGS, CGS, 0);
GSdistd=c2d(GSdist, Ts);
EGSd=GSdistd.B;
%  Decentralised controllers u_i= -K_i x_i 
K1=[10.5 0; 0.90 -0.05];     % cf. Stöcker p. 118
K2=[11.5 0; 1.1 0.40];
A1bar=A1-B1*K1;
A2bar=A2-B2*K2;
%Ks1=inv(B1'*B1)*B1'*Es1;
%Ks2=inv(B2'*B2)*B2'*Es2;
KGS=[K1 Null; Null K2];
AGSbar=AGS-BGS*KGS;
ContrOverallSyst=ss(AGSbar, EGS, CGS, 0);
%  control input generators
I1=ss(A1bar, Es1, C1, 0);
I1d=c2d(I1, Ts);
A1bard=I1d.A;
Es1d=I1d.B;
I2=ss(A2bar, Es2, C2, 0);
I2d=c2d(I2, Ts);
A2bard=I2d.A;
Es2d=I2d.B;



%%  Disturbance behaviour of the continuous control loop
%
figure(10)
%  disturbance
t1=250;
kt1=ceil(t1/Ts);
t2=500;
kt2=ceil(t2/Ts);
Dist1=zeros(kend,1);
for k1=kt1:kt2
    Dist1(k1,1)=0.2;
end
t3=450;
kt3=ceil(t3/Ts);
t4=700;
kt4=ceil(t4/Ts);
Dist2=zeros(kend,1);
for k1=kt3:kt4
    Dist2(k1,1)=-0.4;
end
DistGS=[Dist1 Dist2];
%   behaviour of the overall system
X=lsim(ContrOverallSyst, DistGS, Time, x0);
%
subplot(9,1,[1,2]);
plot(Time, DistGS(:,1), 'b', Time, DistGS(:,2), 'b--');
hold on
latexylabel('$$d_1, d_2\;$$');
rotateY;
xleer;
yticks([-0.2 0 0.2]);
axis([0 Tend -0.45 0.25]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(9,1,[3,4,5])
plot(Time, X(:,1),'b', Time, X(:,3), 'b--');
latexylabel('$$x_{11}, x_{21}$$ in cm');
rotateY;
xleer;
yticks([-1 0 3]);
axis([0 Tend -0.5 3.2]);
hold off
%
subplot(9,1,[6,7,8])
plot(Time, X(:,2), 'b', Time, X(:,4), 'b--');
latexylabel('$$x_{12}, x_{22}$$ in K$$\;$$');
rotateY;
xleer;
yticks([-5 0 5]);
axis([0 Tend -7 6]);
hold off





%%  Event-triggered decentralised controller
%
figure(11)
ebar=0.8;
wx1=25;
wx2=20;
ebar1=ebar;
ebar2=ebar;
%
X=[];
Xs1=[];
Xs2=[];
U1=[];
U2=[];
Test1=[];
Test2=[];
s1hat=[];
s2hat=[];
X(1,:)=x0';
Xs1(1,:)=x0(1:2,1);
Xs2(1,:)=x0(3:4,1);
Eventtime1=zeros(kend,1);
Eventtime2=zeros(kend,1);
Eventtime1(1,1)=1;
Eventtime2(1,1)=1;
s1hat(1,:)=x0(3:4,1)';
s2hat(1,:)=x0(1:2,1)';
for k1=2:kend
   %  feedback loop 1
    Test1(k1)=max([wx1*abs(X(k1-1,1)-Xs1(k1-1,1)), abs(X(k1-1,2)-Xs1(k1-1,2))]);
   if Test1(k1)>ebar1
         % reset the coupling input estimate s1hat
       s1hat(k1-1,:)=X(k1-1,3:4);
         % reset the state of the input generator  
       Xs1(k1,:)=A1bard*X(k1-1,1:2)'+Es1d*s1hat(k1-1,:)';
       U1(k1-1,:)=-K1*X(k1-1,1:2)';
       Eventtime1(k1-1,1)=1;   % event generated by x_1
       s1hat(k1,:)=s1hat(k1-1,:);
   else
       Xs1(k1,:)=A1bard*Xs1(k1-1,:)'+Es1d*s1hat(k1-1,:)';
       U1(k1-1,:)=-K1*Xs1(k1-1,:)';
       s1hat(k1,:)=s1hat(k1-1,:);
   end
   %  feedback loop 2
    Test2(k1)=max([wx2*abs(X(k1-1,3)-Xs2(k1-1,1)), abs(X(k1-1,4)-Xs2(k1-1,2))]);
   if Test2(k1)>ebar2
         % reset the coupling input estimate s2hat
       s2hat(k1-1,:)=X(k1-1,1:2);
         % reset the state of the input generator
       Xs2(k1,:)=A2bard*X(k1-1,3:4)'+Es2d*s2hat(k1-1,:)';
       U2(k1-1,:)=-K2*X(k1-1,3:4)';
       Eventtime2(k1-1,1)=1;   % event generated by x_1
       s2hat(k1,:)=s2hat(k1-1,:);
    else
       Xs2(k1,:)=A2bard*Xs2(k1-1,:)'+Es2d*s2hat(k1-1,:)';
       U2(k1-1,:)=-K2*Xs2(k1-1,:)';
       s2hat(k1,:)=s2hat(k1-1,:);
   end
   X(k1,:)=AGSd*X(k1-1,:)' +BGSd*[U1(k1-1,:) U2(k1-1,:)]' +EGSd*DistGS(k1-1,:)';
end
%
subplot(9,1,[1,2]);
fill([t1 t1 t2 t2 t1],[-0.35 0.25 0.25 -0.35 -0.354],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
plot(Time, DistGS(:,1), 'b');
latexylabel('$$d_1\;$$');
rotateY;
xleer;
yticks([-0.2 0 0.2]);
axis([0 Tend -0.35 0.25]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(9,1,[3,4,5])
fill([t1 t1 t2 t2 t1],[-0.5 2.2 2.2 -0.5 -0.5],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
plot(Time, X(:,1),'b', Time, Xs1(:,1), 'b--');
latexylabel('$$x_{11}$$ in cm');
rotateY;
xleer;
yticks([-1 0 2]);
axis([0 Tend -0.5 2.2]);
hold off
%
subplot(9,1,[6,7,8])
fill([t1 t1 t2 t2 t1],[-6 6 6 -6 -6],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
plot(Time, X(:,2), 'b', Time, Xs1(:,2), 'b--');
latexylabel('$$x_{12}$$ in K$$\;$$');
rotateY;
xleer;
yticks([-5 0 5]);
axis([0 Tend -6 6]);
hold off
%
subplot(9,1,9)
fill([t1 t1 t2 t2 t1],[0 1.5 1.5 0 0],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
for k1=1:kend
    if Eventtime1(k1)==1
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off
%
%  Subsystem 2
%
figure(12)
subplot(9,1,[1,2]);
fill([t3 t3 t4 t4 t3],[-0.45 0.25 0.25 -0.45 -0.45],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
plot(Time, DistGS(:,2), 'b');
latexylabel('$$d_2\;$$');
rotateY;
xleer;
yticks([-0.2 0 0.2]);
axis([0 Tend -0.45 0.25]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(9,1,[3,4,5])
fill([t3 t3 t4 t4 t3],[-0.5 3.2 3.2 -0.5 -0.5],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
plot(Time, X(:,3),'b', Time, Xs2(:,1), 'b--');
latexylabel('$$x_{21}$$ in cm');
rotateY;
xleer;
yticks([-1 0 3]);
axis([0 Tend -0.5 3.2]);
hold off
%
subplot(9,1,[6,7,8])
fill([t3 t3 t4 t4 t3],[-7 6 6 -7 -7],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
plot(Time, X(:,4), 'b', Time, Xs2(:,2), 'b--');
latexylabel('$$x_{22}$$ in K$$\;$$');
rotateY;
xleer;
yticks([-5 0 5]);
axis([0 Tend -7 6]);
hold off
%
subplot(9,1,9)
fill([t3 t3 t4 t4 t3],[0 1.5 1.5 0 0],...
        [0.9 0.9 0.9], 'EdgeColor', 'none');
hold on
for k1=1:kend
    if Eventtime2(k1)==1
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off


%%  Draw the coupling signals and the estimates
%
figure(13)
subplot(6,1,[1,2])
plot(Time, X(:,2), 'b', Time, s2hat(:,2), 'b--');
latexylabel('$$x_{12}, \hat{s}_{21}$$ in K');
rotateY;
xleer;
yticks([ 0 5]);
axis([0 Tend -1 5]);
%
subplot(6,1,[3,4])
plot(Time, X(:,4), 'b', Time, s1hat(:,2), 'b--');
latexylabel('$$x_{22}, \hat{s}_{12}$$ in K');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
yticks([-10 -5 0 5]);
axis([0 Tend -6 6]);


%%  For comparison: Event-triggered decentralised controller without 
%   the communication of the coupling signals
%
figure(14)
%
X=[];
Xs1=[];
Xs2=[];
Test1=[];
Test2=[];
X(1,:)=x0';
Xs1(1,:)=x0(1:2,1);
Xs2(1,:)=x0(3:4,1);
Eventtime1=zeros(kend,1);
Eventtime2=zeros(kend,1);
Eventtime1(1,1)=1;
Eventtime2(1,1)=1;
s1hat=zeros(2,1);
s2hat=zeros(2,1);
for k1=2:kend
   %  feedback loop 1
    Test1(k1)=max([wx1*abs(X(k1-1,1)-Xs1(k1-1,1)), abs(X(k1-1,2)-Xs1(k1-1,2))]);
   if Test1(k1)>ebar1
%          % reset the coupling input estimate s1hat
%        s1hat=X(k1-1,3:4)';
         % reset the state of the input generator  
       Xs1(k1,:)=A1bard*X(k1-1,1:2)'+Es1d*s1hat;
       U1(k1-1,:)=-K1*X(k1-1,1:2)';
       Eventtime1(k1-1,1)=1;   % event generated by x_1
    else
       Xs1(k1,:)=A1bard*Xs1(k1-1,:)'+Es1d*s1hat;
       U1(k1-1,:)=-K1*Xs1(k1-1,:)';
   end
   %  feedback loop 2
    Test2(k1)=max([wx2*abs(X(k1-1,3)-Xs2(k1-1,1)), abs(X(k1-1,4)-Xs2(k1-1,2))]);
   if Test2(k1)>ebar2
%          % reset the coupling input estimate s2hat
%        s2hat=X(k1-1,1:2)';
         % reset the state of the input generator
       Xs2(k1,:)=A2bard*X(k1-1,3:4)'+Es2d*s2hat;
       U2(k1-1,:)=-K2*X(k1-1,3:4)';
       Eventtime2(k1-1,1)=1;   % event generated by x_1
    else
       Xs2(k1,:)=A2bard*Xs2(k1-1,:)'+Es2d*s2hat;
       U2(k1-1,:)=-K2*Xs2(k1-1,:)';
   end
   X(k1,:)=AGSd*X(k1-1,:)' +BGSd*[U1(k1-1,:) U2(k1-1,:)]' +EGSd*DistGS(k1-1,:)';
end
%
subplot(9,1,[1,2]);
plot(Time, DistGS(:,1), 'b');
hold on
latexylabel('$$d_1\;$$');
rotateY;
xleer;
yticks([-0.2 0 0.2]);
axis([0 Tend -0.35 0.25]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(9,1,[3,4,5])
plot(Time, X(:,1),'b', Time, Xs1(:,1), 'b--');
latexylabel('$$x_{11}$$ in cm');
rotateY;
xleer;
yticks([-1 0 2]);
axis([0 Tend -0.5 2.2]);
hold off
%
subplot(9,1,[6,7,8])
plot(Time, X(:,2), 'b', Time, Xs1(:,2), 'b--');
latexylabel('$$x_{12}$$ in K$$\;$$');
rotateY;
xleer;
yticks([-5 0 5]);
axis([0 Tend -6 6]);
hold off
%
subplot(9,1,9)
for k1=1:kend
    if Eventtime1(k1)==1
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off
%
%  Subsystem 2
%
figure(15)
subplot(9,1,[1,2]);
plot(Time, DistGS(:,2), 'b');
hold on
latexylabel('$$d_2\;$$');
rotateY;
xleer;
yticks([-0.2 0 0.2]);
axis([0 Tend -0.45 0.25]);
latextitle('$${\tt VERADemo:}$$ Event-triggered decentralised control');
hold off
%
subplot(9,1,[3,4,5])
plot(Time, X(:,3),'b', Time, Xs2(:,1), 'b--');
latexylabel('$$x_{21}$$ in cm');
rotateY;
xleer;
yticks([-1 0 3]);
axis([0 Tend -0.5 3.2]);
hold off
%
subplot(9,1,[6,7,8])
plot(Time, X(:,4), 'b', Time, Xs2(:,2), 'b--');
latexylabel('$$x_{22}$$ in K$$\;$$');
rotateY;
xleer;
yticks([-5 0 5]);
axis([0 Tend -7 6]);
hold off
%
subplot(9,1,9)
for k1=1:kend
    if Eventtime2(k1)==1
        plot(k1*Ts, 1, 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','w');
        hold on
        plot([k1*Ts k1*Ts],[1 0],'b');
    end
end
axis([0 Tend 0 1.5]);
yleer;
yticks([]);
latexylabel('events');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800 1000]);
hold off







%%   Figures
epsfigc15('VERADemo');
