% System Identification example using data from a quadruple tank process.
% Data in files Y.dat and U.dat. Sampling time Dt=0.1 Sec
% Valve parameters: gamma_1=0.8 and gamma_2=0.7

clear all

load Y.dat
load U.dat

N=length(U)

% Plot the data
figure(1)
subplot(211), plot(U(:,1)), grid, ylabel('u_1')
title('u_1: Volt signal to pump 1')

subplot(212), plot(U(:,2)), grid
title('u_2: Volt signal to pump 2')
xlabel('Samples'), ylabel('u_2')

figure(2)
subplot(211), plot(Y(:,1)), grid, ylabel('y_1')
title('y_1: Level in lower left tank 1')

subplot(212), plot(Y(:,2)), grid, ylabel('y_2')
title('y_2: Level in lower right tank 2')
xlabel('Samples'), ylabel('y_2')

% Calculate an estimate of the nominal outputs corresponding to
% nominal inputs u1=63 and u2=63
u01=mean(U(2000:11000,1))
u02=mean(U(2000:11000,2))
y01=mean(Y(2000:11000,1))   % Nominal level 1
y02=mean(Y(2000:11000,2))   % Nominal level 2

% Create identification data using data centered around the nominal values
Yid=[Y(:,1)-y01 Y(:,2)-y02];
Uid=[U(:,1)-u01 U(:,2)-u02];
Nend=N;
Yid0=Yid(1000:Nend,:); Uid0=Uid(1000:Nend,:); % Cut off first 999 samples


% Using DSR method
L=2;
figure(3)
[A,B,D,E,C,F,x0]=dsr(Yid0,Uid0,L,0);

% Simulate model from input u to output y
Ym=dsrsim(A,B,D,E,Uid,x0);
figure(4)
subplot(211), plot([Yid(:,1) Ym(:,1)]), title('DSR model outputs against real outputs')
subplot(212), plot([Yid(:,2) Ym(:,2)])

er_dsr=Yid-Ym;
L_dsr=er_dsr'*er_dsr/N

% Calculate the optimal predictions from the Kalman filter
Yp=dsropt(A,B,D,E,C,F,Yid,Uid,x0);
figure(5)
subplot(211), plot([Yid(:,1) Yp(:,1)])
subplot(212), plot([Yid(:,2) Yp(:,2)])

% Using PEM function in the IDENT Toolbox
disp('Using PEM')
dat=iddata(Yid0,Uid0); % Pack data into an object
th=pem(dat,4);

%Simulate PEM model
Ym_pem=dsrsim(th.a,th.b,th.c,th.d,Uid,th.x0);
figure(6)
subplot(211), plot([Yid(:,1) Ym_pem(:,1)]), title('PEM model outputs against real outputs')
subplot(212), plot([Yid(:,2) Ym_pem(:,2)])

er_pem=Yid-Ym_pem;
L_pem=er_pem'*er_pem/N




