%   Script PHSENSOR.M
%
%  Überwachung eines pH-Sensors
%
% J. Lunze
% 24.5.2003
% für 2. Auflage: 27.12.2006
% für 3. Auflage: 1.12.2011
% für 4. Auflage: 5.1.2016
% für 5. Auflage: 27.10.2019
%
%
clear
close all
%                   Parameter
%
Ks1=-0.21;
T1=219;
Ks2=3.95;
T0=434;
T2=165;
%
n=2;
A=[-1/T1  0;
   1      -1/T2];
B=[1/T1; 0];
C=[Ks1        0;                        % y1=pH
   Ks2*T0/T2  Ks2*(1-T0/T2)/T2];    % y2=CcO2
D=[0; 0];
CpH=[Ks1  0];                           % y1=pH
CCO2=[Ks2*T0/T2  Ks2*(1-T0/T2)/T2]; % y2=CcO2
DCO2=0;
Bioreaktor=ss(A, B, C, D);

%%  Verhalten des Bioreaktors
%
figure(1)
[Y, T]=step(Bioreaktor);
Tend=1000;                 % bisher Tend=max(T);
%
subplot(4,1,[1,2]);
YFpH(:,1)=Y(:,1)+0.2*ones(size(Y(:,1)));    % fehlerhafte pH-Messung
YFCO2=Y(:,2);
UF=ones(size(YFpH(:,1)));                       % Zulauf
plot(T, Y(:,1), 'b', [-200 0],[0 0], 'b', T,YFpH(:,1),'b--',[-200 0],[0.2 0.2],'b--', [-200 Tend],[0 0],'k:');
latextitle('Bioreaktor bei Zulauferh\"ohung um 1 l/h');
latexylabel('$$pH$$');
rotateY;
xleer;
axis([-100 Tend -0.4 0.4]);
yticks([-0.4 -0.2 0 0.2 0.4]);
%
subplot(4,1,[3,4]);
plot(T, Y(:,2), 'b', [-200 0], [0 0], 'b', [-200 Tend], [0 0], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$c_{\rm CO2}$$ in \%');
rotateY;
axis([-100 Tend -1 5.5]);
yticks([0 2 4 6]);
xticks([0 200 400 600 800 1000]);


%%   Beobachterentwurf
%
figure(2)
%   Beobachter, der Messung CCO2 verwendet
LambdaB=2*[-1/200, -1/200];     % 200 ist Zeitkonst. des Reaktors
LTCO2=acker(A', CCO2', LambdaB);
LCO2=LTCO2';
AbCO2=[A zeros(n,n); LCO2*CCO2 A-LCO2*CCO2];
BbCO2=[B;B];
CbCO2=[CCO2 -CCO2];       % Differenz CCO2 - CCO2_B
DbCO2=[0];
SystemBeobachterCO2=ss(AbCO2, BbCO2, CbCO2, DbCO2);
x0b=[0.3; 0; 0; 0];   %  bei Beobachtungsbeginn pH-Wertdifferenz 0.3
[YbCO2,T,XbCO2]=initial(SystemBeobachterCO2, x0b, T);
%
subplot(6,1,[1,2]);
plot(T, YbCO2, 'b', [0 Tend], [0 0], 'k:');
latexylabel('$$c_{\rm CO2}-\hat{c}_{\rm CO2}$$ in \%');
rotateY;
xleer
latextitle('Beobachtungsfehler bei fehlerfreiem Sensor');
axis([0 Tend -1 3]);
yticks([0 3]);
%   rekonstruierter pH-Wert
subplot(6,1,[3,4]);
plot(T, Ks1*XbCO2(:,3), 'b', [0 Tend], [0 0], 'k:');
latexylabel('$$pH$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([0, Tend, -0.01, 0.03]);
yticks([0 0.03]);
xticks([0 200 400 600 800 1000]);

%% Fehlerdiagnose mit Beobachter
%
figure(3)
%
%    Beobachter für CCO2
ABCO2=A-LCO2*CCO2;
BU=B;           %  Eingang: Zulauf
BYCO2=LCO2;     %  Eingang: gemessener CCO2-Wert
CBCO2=-CCO2;    %  Residuum für pH-Sensor
DBCO2=[0 1];
BeobachterCO2=ss(ABCO2, [BU, BYCO2], CBCO2, DBCO2);
x0CO2=[0.5 0];    % Anfangszustand im pH-Wert um 0.5/Ks1 falsch
%  Beobachter erhält YFCO2 und UF als Eingänge
[YBCO2, T, XBCO2]=lsim(BeobachterCO2, [UF, YFCO2], T, x0CO2);
%
subplot(4,1,[1,2]);
%  Residuum wird aus Zustand x1 und Vergleich mit gemessenem pH-Wert gebildet
plot(T,YFpH(:,1)-Ks1*XBCO2(:,1), 'b', [-200, 0, 0], [0, 0, YFpH(1,1)-Ks1*XBCO2(1,1)], 'b', [0 Tend],[0, 0],'k:');
latexylabel('$$r$$');
rotateY;
axis([-100 Tend -0.05 0.4]);
yticks([0 0.2 0.4]);
xticks([0 200 400 600 800 1000]);
latexxlabel('$$t$$ in min');
latextitle('pH-Sensordiagnose am Biogasreaktor');



%%   Experimentdaten über Zeitraum -200 min... 1000 min
%
figure(4)
dt=10;
Tminus=200;
Tex=[-Tminus:dt:Tend];
feed=[zeros(Tminus/dt, 1); ones(Tend/dt+1, 1)];
feed=feed+0.05*randn(size(feed));
% Rekonstruktion des Zulaufexperiments (vgl. Riedel: Abschlussbericht)
feedbar=1.5;                % Arbeitspunkt Zulauf
feedApkt=feed+feedbar*ones(size(feed));
%
subplot(6,1,[1,2]);
plot(Tex, feedApkt, 'b', [-100, Tend], [0 0], 'k:');
xleer;
axis([-100 Tend -0.5+feedbar 1.5+feedbar]);
yticks([1 2 3]);
latexylabel('$$u$$ in $$\frac{\rm l}{\rm h}$$');
rotateY;
latextitle('Bioreaktor mit Experimentdaten');
%
Yex=lsim(Bioreaktor, feed, Tex);
Yex(:,1)=Yex(:,1)+0.01*randn(size(Yex(:,1)));
Yex(:,2)=Yex(:,2)+0.08*randn(size(Yex(:,2)));
phbar=7.0;              % Arbeitspunkt pH-Wert
YexApkt(:,1)=Yex(:,1)+phbar*ones(size(Yex(:,1)));
CO2bar=34.5;            % Arbeitspunkt CO2-Konzentration
YexApkt(:,2)=Yex(:,2)+CO2bar*ones(size(Yex(:,2)));
%
subplot(6,1,[3,4]);
YFpHex=Yex(:,1)+0.2*ones(size(Yex(:,1)));   % fehlerhafte pH-Messung
YFpHexApkt=YFpHex(:,1)+phbar*ones(size(YFpHex));
YFCO2ex=Yex(:,2);
UFex=feed;                      % Zulauf
plot(Tex, YexApkt(:,1), 'b', Tex,YFpHexApkt(:,1),'b--', [-200 Tend],[0 0],'k:');
latexylabel('$$pH$$');
rotateY;
xleer;
axis([-100 Tend 6.7 7.3]);
yticks([6.6 6.8 7.0 7.2 7.4]);
%
subplot(6,1,[5,6]);
plot(Tex, YexApkt(:,2), 'b', [-200 0], [0 0], [-200 Tend], [0 0], 'k:');
latexxlabel('$$t$$ in min');
latexylabel('$$c_{\rm CO2}$$ in \%');
rotateY;
axis([-100 Tend -1+CO2bar 5.5+CO2bar]);
yticks([34 36 38 40]);
xticks([0 200 400 600 800 1000]);

%% Fehlerdiagnose mit Beobachter mit Experimentdaten
%
figure(5)
Tbeob=[0:dt:Tend];
% Ausblenden der Daten für t>0
k1=Tminus/dt+1;
kend=length(Tex);
YFpHexred=YFpHex(k1:kend, 1);
Yexred=Yex(k1:kend, 1);
YFCO2exred=YFCO2ex(k1:kend, 1);
UFexred=feed(k1:kend, 1);                       % Zulauf

%
%    Beobachter für CCO2
x0CO2=[0.5 0];
%  Beobachter erhält YFCO2 und UF als Eingänge
[YBCO2ex, T, XBCO2ex]=lsim(BeobachterCO2, [UFexred, YFCO2exred], Tbeob, x0CO2);

%
subplot(6,1,[1,2]);
%  Residuum wird aus Zustand x1 und Vergleich mit gemessenem pH-Wert gebildet
plot(T,YFpHexred(:,1)-Ks1*XBCO2ex(:,1), 'b',...
          [-200, 0, 0], [0, 0, YFpHexred(1,1)-Ks1*XBCO2ex(1,1)], 'b',...
           [0 Tend],[0, 0],'k:');
latexylabel('$$r$$');
rotateY;
latextitle('pH-Sensordiagnose am Biogasreaktor mit Experimentdaten');
axis([-100 Tend -0.1 0.4]);
yticks([-0.2 0 0.2 0.4 0.6 0.8]);
xleer
%
%   Erklärung des Diagnoseergebnisses
%
%
subplot(6,1,[3,4]);
%   Wiederholung des mittleren Bildes aus der vorherigen Abbildung
plot(Tex, YexApkt(:,1),  'b',...
            Tex,YFpHexApkt(:,1),'b--', [-200 Tend],[0 0],'k:');
hold on
%   Darüberzeichnen des Beobachtungsergebnisses
plot(T,7+Ks1*XBCO2ex(:,1),'b--');
latexylabel('$$pH$$');
rotateY;
latexxlabel('$$t$$ in min');
axis([-100 Tend 6.7 7.3]);
yticks([6.6 6.8 7.0 7.2 7.4]);
xticks([0 200 400 600 800 1000]);



%%     Bilder in .EPS-Dateien
%
Epsbildc15('phsens');
