%   Script BatchTanks.M 
%   
%    Batchprozess mit 2 gekoppelten Tanks
%    mit ereignisdiskretem Modell
%
% J. Lunze  
% für 3. Auflage: 17.1.2012
% für 5. Auflage: 25.10.2019
%
%
clear 
close all
BspTanks;
%
Tmax=850;   % s
dt=0.1;
T=0:dt:Tmax;
kmax=length(T);  % Anzahl der Zeitpunkte, beginnend bei k=1: t=0;
%
%   Anfangszustand
h10=0;
h20=0;
q0=1;    %  Zustand der Steuerung
k2=0;    %  Ereigniszähler
%
%   Modell des geregelten Systems
h1(1)=h10;
h2(1)=h20;
q(1)=q0;
u1(1)=1;
u2(1)=0;
u3(1)=0;
for k1=1:kmax-1
    q1(k1)=q1max*u1(k1);
    q2(k1)=s2*sqrt(h1(k1))*u2(k1);
    q3(k1)=s3*sqrt(h2(k1))*u3(k1);
    if h1(k1)<hv
        if h2(k1)<hv
            q4(k1)=0;
        else
            q4(k1)=-s4*sqrt(h2(k1)-hv);
        end
    else if h2(k1)<hv
            q4(k1)=s4*sqrt(h1(k1)-hv);
        else
            q4(k1)=s4*sign(h1(k1)-h2(k1))*sqrt(abs(h1(k1)-h2(k1)));
        end
    end
    h1(k1+1)=h1(k1)+dt/A1*(q1(k1)-q2(k1)-q4(k1));
    h2(k1+1)=h2(k1)+dt/A2*(q4(k1)-q3(k1));
    if h1(k1+1)<0
        h1(k1+1)=0;
    end
    if h2(k1+1)<0
        h2(k1+1)=0;
    end
    %      Ablaufsteuerung
    if q(k1)==1
        % Test, ob beide Behälter vollständig gefüllt sind
        if h1(k1)>=L1high && h2(k1)>=L2high
            k2=k2+1;
            e(k2)=k1;              % Ereigniszeitpunkte speichern
            q(k1+1)=2;             % neuer Reglerzustand
            u1(k1+1)=0;
            u2(k1+1)=0;
            u3(k1+1)=1;
        else
            q(k1+1)=1;
            u1(k1+1)=1;
            u2(k1+1)=0;
            u3(k1+1)=0;
        end
    end
    if q(k1)==2
        % Test, ob rechter Behälter leer ist
        if h2(k1)<=L2low
            k2=k2+1;
            e(k2)=k1;              % Ereigniszeitpunkte speichern
            q(k1+1)=3;             % neuer Reglerzustand
            u1(k1+1)=0;
            u2(k1+1)=1;
            u3(k1+1)=0;
        else
            q(k1+1)=2;
            u1(k1+1)=0;
            u2(k1+1)=0;
            u3(k1+1)=1;
        end
    end
    if q(k1)==3
        % Test, ob linker Behälter leer ist
        if h1(k1)<=L1low
            k2=k2+1;
            e(k2)=k1;              % Ereigniszeitpunkte speichern
            q(k1+1)=1;             % neuer Reglerzustand
            u1(k1+1)=1;
            u2(k1+1)=0;
            u3(k1+1)=0;
        else
            q(k1+1)=3;
            u1(k1+1)=0;
            u2(k1+1)=1;
            u3(k1+1)=0;
        end
    end
end
%   Angleichen der Länge der Signalvektoren
k1=kmax;
q1(k1)=q1max*u1(k1);
q2(k1)=s2*sqrt(h1(k1))*u2(k1);
q3(k1)=s3*sqrt(h2(k1))*u3(k1);
if h1(k1)<hv
    if h2(k1)<hv
        q4(k1)=0;
    else
        q4(k1)=-s4*sqrt(h2(k1)-hv);
    end
else if h2(k1)<hv
        q4(k1)=s4*sqrt(h1(k1)-hv);
    else
        q4(k1)=s4*sign(h1(k1)-h2(k1))*sqrt(abs(h1(k1)-h2(k1)));
    end
end
figure(1)
subplot(4,1,[1, 2])
latextitle('Gesteuerter Batchprozess');
plot(T, h1, 'b');
hold on
plot([0 Tmax],[L1low, L1low],'k:');
%plot([0 Tmax],[L1med, L1med],'k:');
plot([0 Tmax],[L1high, L1high],'k:');
%     Markieren der Ereigniszeitpunkte
for k3=1:k2
    te=(e(k3)-1)*dt;
    plot([te, te],[0 0.7],'k:');
end
xleer;
latexylabel('$$h_1$$ in m');
rotateY;
%
plot(T, h2,'r');
hold on
plot([0 Tmax],[L2low, L2low],'r--');
%plot([0 Tmax],[L2med, L2med],'r--');
plot([0 Tmax],[L2high, L2high],'r--');
latextext(Tmax+10, L1low, '$$l_1, l_3$$');
latextext(Tmax+10, L1high, '$$l_2$$');
latextext(Tmax+10, L2high, '$$l_4$$');
xleer;
yleer;
latextext(220, 0.36, '$$h_1$$');
latextext(150, 0.2, '$$h_2$$');
%latexylabel('$$h1, h_2$$'); 
rotateY;
axis([0 Tmax 0 0.68]);
hold off
%
%
%    Bestimmung des diskreten Zustands bei zwei Sensoren
for k1=1:length(h1)
    if h1(k1)<L1low
        if h2(k1)<L2low
            hqual(k1)=1;
        else if h2(k1)<L2high
                hqual(k1)=4;
            else 
                hqual(k1)=7;
            end
        end
    end
    if h1(k1)>= L1low && h1(k1)<L1high
        if h2(k1)<L2low
            hqual(k1)=2;
        else if h2(k1)<L2high
                hqual(k1)=5;
            else 
                hqual(k1)=8;
            end
        end
    end
    if h1(k1)>= L1high
        if h2(k1)<L2low
            hqual(k1)=3;
        else if h2(k1)<L2high
                hqual(k1)=6;
            else 
                hqual(k1)=9;
            end
        end
    end
end
%
subplot(4,1,3)
plot(T, hqual, 'b')
hold on
for k3=1:k2
    te=(e(k3)-1)*dt;
    plot([te, te],[0 10],'k:');
end
xleer
latexylabel('$$z$$');
rotateY;
yticks([1, 9]);
axis([0 Tmax 0 10]);
hold off;
%
subplot(4,1,4)
plot(T, q, 'b')
hold on
for k3=1:k2
    te=(e(k3)-1)*dt;
    plot([te, te],[0 4],'k:');
end
latexxlabel('$$t$$ in s');
xticks([0 200 400 600 800]);
latexylabel('$$z_{\rm s}$$');
rotateY;
yticks([1,  3]);
axis([0 Tmax 0.5 3.5]);
xticks([0 200 400 600 800]);
hold off;

%%
%
Epsbildc15('Batchtanks');

% % %    Bestimmung des diskreten Zustands bei drei Sensoren
% % for k1=1:length(h1)
% %     if h1(k1)<L1low
% %         if h2(k1)<L2low
% %             hqual(k1)=1;
% %         else if h2(k1)<L2med
% %                 hqual(k1)=5;
% %             else if h2(k1)<L2high
% %                 hqual(k1)=9;
% %                 else
% %                     hqual(k1)=13;
% %                 end
% %             end
% %         end
% %     end
% %     if h1(k1)>= L1low && h1(k1)<L1med
% %         if h2(k1)<L2low
% %             hqual(k1)=2;
% %         else if h2(k1)<L2med
% %                 hqual(k1)=6;
% %             else if h2(k1)<L2high
% %                 hqual(k1)=10;
% %                 else
% %                     hqual(k1)=14;
% %                 end
% %             end
% %         end
% %     end
% %     if h1(k1)>= L1med && h1(k1)<L1high
% %         if h2(k1)<L2low
% %             hqual(k1)=3;
% %         else if h2(k1)<L2med
% %                 hqual(k1)=7;
% %             else if h2(k1)<L2high
% %                 hqual(k1)=11;
% %                 else
% %                     hqual(k1)=15;
% %                 end
% %             end
% %         end
% %      end
% %     if h1(k1)>= L1high
% %         if h2(k1)<L2low
% %             hqual(k1)=4;
% %         else if h2(k1)<L2med
% %                 hqual(k1)=8;
% %             else if h2(k1)<L2high
% %                 hqual(k1)=12;
% %                 else
% %                     hqual(k1)=16;
% %                 end
% %             end
% %         end
% %     end
% % end

            
