% tmp_demo2
% Demo som viser systemidentifikasjon av effekt og konsistens 
% på en TMP raffinør ved Norske Skog Union.
% Vi studerer her sammenhengen fra settpunkt til de lokale
% regulatorene til utgangsvariablene (raffinøreffekt og
% konsistens).

clear all
load r1sd1710.txt

%% Inngangsdata
% Aktuelle paadragsvariable (process values)
u1=r1sd1710(:,2);   % rpm
u2=r1sd1710(:,6);   % dilution 
u3=r1sd1710(:,11);

% Settpunkt paadragsvariable
u1s=r1sd1710(:,23)*4.75;
u2s=r1sd1710(:,22);
u3s=r1sd1710(:,24);

%% Utgangsdata
y1=r1sd1710(:,7);    % power [MW]
y2=r1sd1710(:,10);   % consistency [1]

Lf=0.0529*3.6*u1;    % Production [tonn/h]
y3=y1./Lf;           % Specific energy [MWh/tonn]

%%%%%%%%%%%%%%
idata=1;
if idata==1
   U=[u1s u2s u3s];
   Y=[y1 y2];
elseif idata==0
   Utmp=[u1s u2s u3s];Ytmp=[y1 y2];  % data for øving4-> utmp2.dat og ytmp2.dat
   u1s0=mean(u1s(1:25)); u2s0=mean(u2s(1:25)); u3s0=mean(u3s(1:25)); % Trends
   y10=mean(y1(1:25)); y20=mean(y2(1:25));
   U=[u1s-u1s0 u2s-u2s0 u3s-u3s0];
   Y=[y1-y10 y2-y20];
   Uid=U(1:1500,:);
   Yid=Y(1:1500,:);
end

[a,b,d,e,c,f,x0]=dsr(Yid,Uid,1);
Ys=dsrsim(a,b,d,e,U,x0);

%% Plot results and data.
figure(1), clf
subplot(211), plot([Y(:,1)+y10 Ys(:,1)+y10])
grid, xlabel('Time [samples]'), ylabel('y_1')
title('Motor power: Actual and simulated output #1: First 1500 samples used by DSR.')
subplot(212), plot([Y(:,2)+y20 Ys(:,2)+y20])
grid, xlabel('Time [samples]'), ylabel('y_2')
title('Cons.: Actual and simulated output #2: First 1500 samples used by DSR.')

figure(2), clf
subplot(321), plot(u1s)
grid, ylabel('u_1')
title('u_1: dosage screw speed setpoint [rpm].')
axis([0 2400 50 60]) 

subplot(322), plot(u1)
grid, ylabel('u_1^pv')
title('Actual dosage screw speed.')
axis([0 2400 50 60]) 

subplot(323), plot(u2s)
grid, ylabel('u_2')
title('u_2: dilution water setpoint [kg/s].')
axis([0 2400 6.5 8.5]) 

subplot(324), plot(u2)
grid, ylabel('u_2^pv')
title('Actual dilution water.')
axis([0 2400 6.5 8.5]) 

subplot(325), plot(u3s)
grid, xlabel('Time [samples]'), ylabel('u_3')
title('u_3: plate gap setpoint [mm].')
axis([0 2400 0.4 0.8]) 

subplot(326), plot(u3)
grid, xlabel('Time [samples]'), ylabel('u_3^pv')
title('Actual plate gap.')
axis([0 2400 0.4 0.8]) 
