clear all
close all
clc

load heater

y1=ye-mean(ye);
u1=ue-mean(ue);
y2=yv-mean(yv);
u2=uv-mean(uv);

NetDef = ['HHHHHH';'L-----'];
NN=[2 2 4];
trparms=settrain;
trparms=settrain(trparms,'maxiter',500);

%[W1,W2]=nnarx(NetDef,NN,[],[],trparms,y1,u1);
[W1,W2]=nnoe(NetDef,NN,[],[],trparms,y1,u1);

a1=W1(:,1:end-1); 
b1=W1(:,end); 
a2=W2(1:end-1); 
b2=W2(end);

ys=zeros(size(y2));
for t=5:length(y2)
    % w=[y2(t);y2(t-1);u2(t-2);u2(t-3)]; % 1-step prediction
    w=[ys(t);ys(t-1);u2(t-3);u2(t-4)]; % simulation
    ys(t+1)=a2*tanh(a1*w+b1)+b2;
end

RMSE=norm(y2(101:1000)-ys(101:1000))/sqrt(900)

figure
subplot(2,1,1),hold on,box on
plot(y2)
plot(ys,'r')
xlim([0 500])
subplot(2,1,2)
plot(u2)
xlim([0 500])

