% DEMO_TMP
% Demonstrate the D-SR Toolbox functions, DSR.M, DSRSIM.M and DSROPT.M.
% The data is from a thermo mechanical pulping refiner at Union Co., Skien.

load ytmp.dat;                                % Get data from ascii file.
load utmp.dat;
Y=ytmp;
U=utmp;
[N,m]=size(Y); [N,r]=size(U);

L=3;
L=dread('Choose state prediction horizon L',L);

[a,b,d,e,c,f,x0]=dsr(Y,U,L);                  % Identify state space model.

yp_d=dsrsim(a,b,d,e,U,x0);                    % Compute the simulated output, yp_d.

figure(1), clf                                % Plot simulated outputs vs. actual.
subplot(211), plot([Y(:,1) yp_d(:,1)])
grid, title('y1: actual and simulated output')
subplot(212), plot([Y(:,2) yp_d(:,2)])
grid, title('y2: actual and simulated output')

ed=Y-yp_d;                                    % Simulated error 
Vd = ed'*ed/(N-1)                             % Simulated error criterion.
trVd=trace(Vd)

yp_s=dsropt(a,b,d,e,c,f,Y,U,x0);              % Compute 1-step ahead predictions.

figure(2), clf                                % Plot predictions vs. actual outputs.
subplot(211), plot([Y(:,1) yp_s(:,1)])
grid, title('y1: actual and optimal 1-step ahead prediction')
subplot(212), plot([Y(:,2) yp_s(:,2)])
grid, title('y2: actual and optimal 1-step ahead prediction')

es=Y-yp_s;                                    % Prediction error 
Vs = es'*es/(N-1)                             % Prediction error criterion.
trVs=trace(Vs)

%
figure(3)
subplot(311), plot(U(:,1)),
grid, ylabel('u_1')
title('u_1: Actual dosage screw speed [rpm].')

subplot(312), plot(U(:,2))
grid, ylabel('u_2')
title('u_2: actual dilution water [kg/s].')

subplot(313), plot(U(:,3))
grid, xlabel('Time [samples]'), ylabel('u_3')
title('u_3: actual plate gap [mm].')