% Example of identification with ANN - 1st order system

% process model: TANH process
% mathematical model of the process with all parameters
% y(k)=y(k-1)+Ts*(-tanh(y(k-1)+u(k-1)^3));

% simulation scheme, that enables rerun of data acquisition  (no masked parts);
% y(k)= dtanh([y(k-1) u(k-1)]);

clear all
close all
% demonstration of process nonlinearity with nonlinearity plane
Ts=0.5;
u=(-3:0.5:3);
y=(-3:0.5:3);
yp=[];
for i=1:length(y);
   for j=1:length(u);
  yp(j,i) = dtanh([y(i) u(j)]); 
end;
end;
surf(u,y,yp)
axis([-3 3 -3 3 -4 4])
xlabel('u(k)'),ylabel('y(k)'),zlabel('y(k+1)')

% identification input signal and corresponding response
load input
yi=0;
for k=1:length(ui)
yi(k+1)=dtanh([yi(k) ui(k)]);
end
t=(0:0.5:0.5*length(ui)-0.5);
yi1=yi(2:length(yi))';
yi=yi(1:length(ui))';
figure
subplot(211)
plot(t,ui)
title('Input signal for identification')
xlabel('Time'),grid
subplot(212)
plot(t,yi)
title('System response for identification')
xlabel('Time'),grid

% signal for validation and corresponding response
yv=0;
for k=1:length(uv)
yv(k+1)=dtanh([yv(k) uv(k)]);
end
t=(0:0.5:0.5*length(uv)-0.5)';
yv=yv(1:length(uv))';
figure
subplot(211)
plot(t,uv)
title('Input signal for validation')
xlabel('Time'),grid
subplot(212)
plot(t,yv)
title('System response for validation')
xlabel('Time'),grid

% histogram of magnitudes

figure
hist(ui,20)
title('Magnitude distribution of identification input signal')
xlabel('Signal magnitude')

figure
hist(uv,20)
title('Magnitude distribution of validation input signal')
xlabel('Signal magnitude')

% data distribution
figure
mesh(u,y,yp)
axis([-3 3 -3 3 -4 4])
xlabel('u(k)'),ylabel('y(k)'),zlabel('y(k+1)')

hold on
plot3(yi,ui,yi1,'*k')
hold off

% identification

% first step: scaling of identification signal to mean value  0 and variance 1, 
% afterwards scaling of validation signal with the same constants! 
[ui_s,usk]=dscale(ui');
[yi_s,ysk]=dscale(yi');

uv_s=dscale(uv',usk);
yv_s=dscale(yv',ysk);


NetDef=['HHHHH';'L----'];
NN=[1 1 1];
trparms = settrain;
trparms = settrain(trparms,'maxiter',500,'D',1e-2,'skip',10);

[w1,w2,critvec,iteration,lambda]=nnarx(NetDef,NN,ones(5,3),ones(1,6),trparms,yi_s,ui_s);


%rescaling of weights back on original signals
[W1,W2]=wrescale('nnarx',w1,w2,usk,ysk,NN);

% figure of the final network structure
figure
drawnet(W1,W2);

% matrices of weights
W1
W2

% validation of response residuals obtained with simulation on validation signal 
% comparison of planes
figure
ym=[];
[u,y] = meshgrid(-3:0.5:3, -3:0.5:3);
siz=size(y);
input=[];
for j=1:siz(2)
    for i=1:siz(1)
     input=[input;u(i,j),y(i,j)]; 
    end
end
[yhat,res,MSE] = nneval(NetDef,W1,W2,input',zeros(1,length(input)),1);
for k=1:length(yhat)
    [i,j]=ind2sub(siz,k);
    ym(i,j)=yhat(k);
end
err=[];
err=yp-ym;
surf(u,y,ym)
axis([-3 3 -3 3 -4 4])
xlabel('u(k)'),ylabel('y(k)'),zlabel('y(k+1)')
figure
surf(u,y,err)
axis([-3 3 -3 3 -1 1])
xlabel('u(k)'),ylabel('y(k)'),zlabel('Residuals')

figure
contour(u,y,err,16);
axis([-3 3 -3 3])
xlabel('u(k)'),ylabel('y(k)'),title('Residuals')

% process response obtained with simulation on validation signal
warning off MATLAB:break_outside_of_loop;
yodz=nnsimul('nnarx',NetDef,NN,W1,W2,yv',uv');
figure
plot(t,yv,'r',t(2:length(t)),yodz','b');
title('System response on validation signal - simulation')
xlabel('Time'),grid

% histogram of residuals
e=yv(1:length(yv)-1)'-yodz;
figure
hist(e,20)
title('Histogram of residuals')

% autocorrelation of residuals

coefEE = crossco(e,e,180);
Tau=(-90:0.5:90);
figure;
plot(Tau,coefEE);
grid on;
xlabel('Tau [s]');
title('Autocorrelation of residuals');

% cross-correlation between residuals and input signal
coefUE = crossco(uv_s(1:length(e)),e,180);
figure;
plot(Tau,coefUE);
grid on;
axis([-100 100 -1 1]);
xlabel('Tau [s]');
title('Cross-correlation between residuals and input signal');
