Laboratory 6 - Estimation, filtering and system identification - Prof. M. Taragna

Exercise: model identification for a water heater, using nonlinear models (real data)

Introduction
Procedure
Setup
Neural Network model class
Choosing the best model

Introduction

The program code may be splitted in sections using the characters "%%". Each section can run separately with the command "Run Section" (in the Editor toolbar, just to the right of the "Run" button). You can do the same thing by highlighting the code you want to run and by using the button function 9 (F9). This way, you can run only the desired section of your code, saving your time. This script can be considered as a reference example.

clear all, close all

Procedure

Load the .mat file containing input-output data
Remove the mean value from the data
Plot input-output data
For each input/output delay nk from 1 to 5, and for each order na=nb (NNARX) from 1 to 5
Identify the NNARX model using the estimation dataset
Compute the predicted output on the validation dataset
Compute RMSE and MDL criteria in order to select the most suitable model

Setup

% Step 1: load of data

load heater

% Step 2: remove the mean value

ue_0mean = ue - mean(ue);
ye_0mean = ye - mean(ye);
uv_0mean = uv - mean(uv);
yv_0mean = yv - mean(yv);

% Step 3: plot of input-output data

Ne = length(ue); % number of estimation data
Nv = length(uv); % number of validation data

figure, plot(ue_0mean), ylabel('ue\_0mean'),
figure, plot(ye_0mean), ylabel('ye\_0mean'),
figure, plot(uv_0mean), ylabel('uv\_0mean'),
figure, plot(yv_0mean), ylabel('yv\_0mean')

Neural Network model class

Define the model structure of the neural network having:

r>=5 hyperbolic tangent neurons (or nodes or units) in the hidden layer
1 linear neuron (or node or unit) in the output layer

NetDef=['HHHHHH';'L-----'];

% For any Neural Network ARX model NNARX(na,nb,nk), define its orders:
% - na = # of past measurements used in the regressor phi
% - nb = # of past inputs used in the regressor phi
% - nk = minimum input-output delay

% Step 4: definition of NNARX parameters

na=1;
nb=na;
nk=1;

Set the maximum number of iterations to 500

trparms=settrain; trparms=settrain(trparms,'maxiter',500);

Perform the NNARX identification, using

as ouput signal Y, a row vector with dim(Y) = [1 * (# of data)]
as input signal U, a matrix with dim(U) = [(# of inputs) * (# of data)] In the SISO case, both Y and U are row vectors

% Step 5: estimation of NNARX parameters

[W1,W2]=nnarx(NetDef,[na,nb,nk],[],[],trparms,ye_0mean',ue_0mean');

Network training started at 14. 9.49
iteration # 1   W = 9.962e+00iteration # 2   W = 5.950e+00iteration # 3   W = 5.755e+00iteration # 4   W = 5.661e+00iteration # 5   W = 5.140e+00iteration # 6   W = 5.011e+00iteration # 7   W = 4.832e+00iteration # 8   W = 4.542e+00iteration # 9   W = 4.190e+00iteration # 10   W = 3.863e+00iteration # 11   W = 3.861e+00iteration # 12   W = 3.578e+00iteration # 13   W = 3.508e+00iteration # 14   W = 3.464e+00iteration # 15   W = 3.428e+00iteration # 16   W = 3.365e+00iteration # 17   W = 3.248e+00iteration # 18   W = 3.032e+00iteration # 19   W = 2.673e+00iteration # 20   W = 2.129e+00iteration # 21   W = 1.565e+00iteration # 22   W = 1.115e+00iteration # 23   W = 8.778e-01iteration # 24   W = 7.516e-01iteration # 25   W = 6.907e-01iteration # 26   W = 6.567e-01iteration # 27   W = 6.387e-01iteration # 28   W = 6.302e-01iteration # 29   W = 6.222e-01iteration # 30   W = 6.206e-01iteration # 31   W = 6.189e-01iteration # 32   W = 6.184e-01iteration # 33   W = 6.178e-01iteration # 34   W = 6.177e-01iteration # 35   W = 6.173e-01iteration # 36   W = 6.172e-01iteration # 37   W = 6.171e-01iteration # 38   W = 6.168e-01iteration # 39   W = 6.165e-01iteration # 40   W = 6.159e-01iteration # 41   W = 6.157e-01iteration # 42   W = 6.155e-01iteration # 43   W = 6.152e-01iteration # 44   W = 6.148e-01iteration # 45   W = 6.140e-01iteration # 46   W = 6.140e-01iteration # 47   W = 6.124e-01iteration # 48   W = 6.118e-01iteration # 49   W = 6.111e-01iteration # 50   W = 6.107e-01iteration # 51   W = 6.102e-01iteration # 52   W = 6.100e-01iteration # 53   W = 6.095e-01iteration # 54   W = 6.093e-01iteration # 55   W = 6.091e-01iteration # 56   W = 6.088e-01iteration # 57   W = 6.085e-01iteration # 58   W = 6.082e-01iteration # 59   W = 6.080e-01iteration # 60   W = 6.079e-01iteration # 61   W = 6.079e-01iteration # 62   W = 6.077e-01iteration # 63   W = 6.077e-01iteration # 64   W = 6.076e-01iteration # 65   W = 6.075e-01iteration # 66   W = 6.075e-01iteration # 67   W = 6.074e-01iteration # 68   W = 6.074e-01iteration # 69   W = 6.073e-01iteration # 70   W = 6.073e-01iteration # 71   W = 6.073e-01iteration # 72   W = 6.072e-01iteration # 73   W = 6.072e-01iteration # 74   W = 6.072e-01iteration # 75   W = 6.071e-01iteration # 76   W = 6.071e-01iteration # 77   W = 6.071e-01iteration # 78   W = 6.071e-01iteration # 79   W = 6.070e-01iteration # 80   W = 6.070e-01iteration # 81   W = 6.070e-01iteration # 82   W = 6.070e-01iteration # 83   W = 6.069e-01iteration # 84   W = 6.069e-01iteration # 85   W = 6.069e-01iteration # 86   W = 6.068e-01iteration # 87   W = 6.068e-01iteration # 88   W = 6.068e-01iteration # 89   W = 6.067e-01iteration # 90   W = 6.067e-01iteration # 91   W = 6.066e-01iteration # 92   W = 6.066e-01iteration # 93   W = 6.065e-01iteration # 94   W = 6.065e-01iteration # 95   W = 6.064e-01iteration # 96   W = 6.063e-01iteration # 97   W = 6.063e-01iteration # 98   W = 6.061e-01iteration # 99   W = 6.061e-01iteration # 100   W = 6.060e-01iteration # 101   W = 6.059e-01iteration # 102   W = 6.058e-01iteration # 103   W = 6.057e-01iteration # 104   W = 6.055e-01iteration # 105   W = 6.054e-01iteration # 106   W = 6.053e-01iteration # 107   W = 6.053e-01iteration # 108   W = 6.053e-01iteration # 109   W = 6.053e-01iteration # 110   W = 6.053e-01iteration # 111   W = 6.053e-01iteration # 112   W = 6.053e-01iteration # 113   W = 6.053e-01iteration # 114   W = 6.053e-01iteration # 115   W = 6.053e-01iteration # 116   W = 6.053e-01iteration # 117   W = 6.053e-01iteration # 118   W = 6.053e-01iteration # 119   W = 6.053e-01iteration # 120   W = 6.053e-01iteration # 121   W = 6.053e-01iteration # 122   W = 6.052e-01iteration # 123   W = 6.052e-01iteration # 124   W = 6.052e-01iteration # 125   W = 6.052e-01iteration # 126   W = 6.052e-01iteration # 127   W = 6.052e-01iteration # 128   W = 6.052e-01iteration # 129   W = 6.052e-01iteration # 130   W = 6.052e-01iteration # 131   W = 6.052e-01iteration # 132   W = 6.052e-01iteration # 133   W = 6.052e-01iteration # 134   W = 6.052e-01iteration # 135   W = 6.051e-01iteration # 136   W = 6.051e-01iteration # 137   W = 6.051e-01iteration # 138   W = 6.051e-01iteration # 139   W = 6.051e-01iteration # 140   W = 6.051e-01iteration # 141   W = 6.051e-01iteration # 142   W = 6.051e-01iteration # 143   W = 6.051e-01iteration # 144   W = 6.051e-01iteration # 145   W = 6.051e-01iteration # 146   W = 6.051e-01iteration # 147   W = 6.051e-01iteration # 148   W = 6.051e-01iteration # 149   W = 6.051e-01iteration # 150   W = 6.051e-01iteration # 151   W = 6.051e-01iteration # 152   W = 6.051e-01iteration # 153   W = 6.051e-01iteration # 154   W = 6.051e-01iteration # 155   W = 6.051e-01iteration # 156   W = 6.050e-01iteration # 157   W = 6.050e-01iteration # 158   W = 6.050e-01iteration # 159   W = 6.050e-01iteration # 160   W = 6.050e-01iteration # 161   W = 6.050e-01iteration # 162   W = 6.050e-01iteration # 163   W = 6.050e-01iteration # 164   W = 6.050e-01iteration # 165   W = 6.050e-01iteration # 166   W = 6.050e-01iteration # 167   W = 6.050e-01iteration # 168   W = 6.050e-01iteration # 169   W = 6.050e-01iteration # 170   W = 6.050e-01iteration # 171   W = 6.050e-01iteration # 172   W = 6.050e-01iteration # 173   W = 6.050e-01iteration # 174   W = 6.050e-01iteration # 175   W = 6.050e-01iteration # 176   W = 6.050e-01iteration # 177   W = 6.050e-01iteration # 178   W = 6.050e-01iteration # 179   W = 6.050e-01iteration # 180   W = 6.050e-01iteration # 181   W = 6.050e-01iteration # 182   W = 6.050e-01iteration # 183   W = 6.050e-01iteration # 184   W = 6.050e-01iteration # 185   W = 6.050e-01iteration # 186   W = 6.050e-01iteration # 187   W = 6.050e-01iteration # 188   W = 6.050e-01iteration # 189   W = 6.050e-01iteration # 190   W = 6.050e-01iteration # 191   W = 6.050e-01iteration # 192   W = 6.050e-01iteration # 193   W = 6.050e-01iteration # 194   W = 6.050e-01iteration # 195   W = 6.050e-01iteration # 196   W = 6.050e-01iteration # 197   W = 6.050e-01iteration # 198   W = 6.050e-01iteration # 199   W = 6.050e-01iteration # 200   W = 6.050e-01iteration # 201   W = 6.050e-01iteration # 202   W = 6.050e-01iteration # 203   W = 6.050e-01iteration # 204   W = 6.050e-01iteration # 205   W = 6.050e-01iteration # 206   W = 6.050e-01iteration # 207   W = 6.050e-01iteration # 208   W = 6.050e-01iteration # 209   W = 6.050e-01iteration # 210   W = 6.050e-01iteration # 211   W = 6.050e-01iteration # 212   W = 6.050e-01iteration # 213   W = 6.050e-01iteration # 214   W = 6.050e-01iteration # 215   W = 6.050e-01iteration # 216   W = 6.050e-01iteration # 217   W = 6.050e-01iteration # 218   W = 6.050e-01iteration # 219   W = 6.050e-01iteration # 220   W = 6.050e-01iteration # 221   W = 6.050e-01iteration # 222   W = 6.050e-01iteration # 223   W = 6.050e-01iteration # 224   W = 6.050e-01iteration # 225   W = 6.050e-01iteration # 226   W = 6.050e-01iteration # 227   W = 6.050e-01iteration # 228   W = 6.050e-01iteration # 229   W = 6.050e-01iteration # 230   W = 6.050e-01iteration # 231   W = 6.050e-01iteration # 232   W = 6.050e-01iteration # 233   W = 6.050e-01iteration # 234   W = 6.050e-01iteration # 235   W = 6.050e-01iteration # 236   W = 6.050e-01iteration # 237   W = 6.050e-01iteration # 238   W = 6.050e-01iteration # 239   W = 6.050e-01iteration # 240   W = 6.050e-01iteration # 241   W = 6.050e-01iteration # 242   W = 6.050e-01iteration # 243   W = 6.050e-01iteration # 244   W = 6.050e-01iteration # 245   W = 6.050e-01iteration # 246   W = 6.050e-01iteration # 247   W = 6.050e-01iteration # 248   W = 6.050e-01iteration # 249   W = 6.050e-01iteration # 250   W = 6.050e-01iteration # 251   W = 6.050e-01iteration # 252   W = 6.050e-01iteration # 253   W = 6.050e-01iteration # 254   W = 6.050e-01iteration # 255   W = 6.050e-01iteration # 256   W = 6.050e-01iteration # 257   W = 6.050e-01iteration # 258   W = 6.050e-01iteration # 259   W = 6.050e-01iteration # 260   W = 6.050e-01iteration # 261   W = 6.050e-01iteration # 262   W = 6.050e-01iteration # 263   W = 6.050e-01iteration # 264   W = 6.050e-01iteration # 265   W = 6.050e-01iteration # 266   W = 6.050e-01iteration # 267   W = 6.050e-01iteration # 268   W = 6.050e-01iteration # 269   W = 6.050e-01iteration # 270   W = 6.050e-01iteration # 271   W = 6.050e-01iteration # 272   W = 6.050e-01iteration # 273   W = 6.050e-01iteration # 274   W = 6.050e-01iteration # 275   W = 6.050e-01iteration # 276   W = 6.050e-01iteration # 277   W = 6.050e-01iteration # 278   W = 6.050e-01iteration # 279   W = 6.050e-01iteration # 280   W = 6.050e-01iteration # 281   W = 6.050e-01iteration # 282   W = 6.050e-01iteration # 283   W = 6.050e-01iteration # 284   W = 6.050e-01iteration # 285   W = 6.050e-01iteration # 286   W = 6.050e-01iteration # 287   W = 6.050e-01iteration # 288   W = 6.050e-01iteration # 289   W = 6.050e-01iteration # 290   W = 6.050e-01iteration # 291   W = 6.050e-01iteration # 292   W = 6.050e-01iteration # 293   W = 6.050e-01iteration # 294   W = 6.050e-01iteration # 295   W = 6.050e-01iteration # 296   W = 6.050e-01iteration # 297   W = 6.050e-01iteration # 298   W = 6.050e-01iteration # 299   W = 6.050e-01iteration # 300   W = 6.050e-01iteration # 301   W = 6.050e-01iteration # 302   W = 6.050e-01iteration # 303   W = 6.050e-01iteration # 304   W = 6.050e-01iteration # 305   W = 6.050e-01iteration # 306   W = 6.050e-01iteration # 307   W = 6.050e-01iteration # 308   W = 6.050e-01iteration # 309   W = 6.050e-01iteration # 310   W = 6.050e-01iteration # 311   W = 6.050e-01iteration # 312   W = 6.050e-01iteration # 313   W = 6.050e-01iteration # 314   W = 6.050e-01iteration # 315   W = 6.050e-01iteration # 316   W = 6.050e-01iteration # 317   W = 6.050e-01iteration # 318   W = 6.050e-01iteration # 319   W = 6.050e-01iteration # 320   W = 6.050e-01iteration # 321   W = 6.050e-01iteration # 322   W = 6.050e-01iteration # 323   W = 6.050e-01iteration # 324   W = 6.050e-01iteration # 325   W = 6.050e-01iteration # 326   W = 6.050e-01iteration # 327   W = 6.050e-01iteration # 328   W = 6.050e-01iteration # 329   W = 6.050e-01iteration # 330   W = 6.050e-01iteration # 331   W = 6.049e-01iteration # 332   W = 6.049e-01iteration # 333   W = 6.049e-01iteration # 334   W = 6.049e-01iteration # 335   W = 6.049e-01iteration # 336   W = 6.049e-01iteration # 337   W = 6.049e-01iteration # 338   W = 6.049e-01iteration # 339   W = 6.049e-01iteration # 340   W = 6.049e-01iteration # 341   W = 6.049e-01iteration # 342   W = 6.049e-01iteration # 343   W = 6.049e-01iteration # 344   W = 6.049e-01iteration # 345   W = 6.049e-01iteration # 346   W = 6.049e-01iteration # 347   W = 6.049e-01iteration # 348   W = 6.049e-01iteration # 349   W = 6.049e-01iteration # 350   W = 6.049e-01iteration # 351   W = 6.049e-01iteration # 352   W = 6.049e-01iteration # 353   W = 6.049e-01iteration # 354   W = 6.049e-01iteration # 355   W = 6.049e-01iteration # 356   W = 6.049e-01iteration # 357   W = 6.049e-01iteration # 358   W = 6.049e-01iteration # 359   W = 6.049e-01iteration # 360   W = 6.049e-01iteration # 361   W = 6.049e-01iteration # 362   W = 6.049e-01iteration # 363   W = 6.049e-01iteration # 364   W = 6.049e-01iteration # 365   W = 6.049e-01iteration # 366   W = 6.049e-01iteration # 367   W = 6.049e-01iteration # 368   W = 6.049e-01iteration # 369   W = 6.049e-01iteration # 370   W = 6.049e-01iteration # 371   W = 6.049e-01iteration # 372   W = 6.049e-01iteration # 373   W = 6.049e-01iteration # 374   W = 6.049e-01iteration # 375   W = 6.049e-01iteration # 376   W = 6.049e-01iteration # 377   W = 6.049e-01iteration # 378   W = 6.049e-01iteration # 379   W = 6.049e-01iteration # 380   W = 6.049e-01iteration # 381   W = 6.049e-01iteration # 382   W = 6.049e-01iteration # 383   W = 6.049e-01iteration # 384   W = 6.049e-01iteration # 385   W = 6.049e-01iteration # 386   W = 6.049e-01iteration # 387   W = 6.049e-01iteration # 388   W = 6.049e-01iteration # 389   W = 6.049e-01iteration # 390   W = 6.049e-01iteration # 391   W = 6.049e-01iteration # 392   W = 6.049e-01iteration # 393   W = 6.049e-01iteration # 394   W = 6.049e-01iteration # 395   W = 6.049e-01iteration # 396   W = 6.049e-01iteration # 397   W = 6.049e-01iteration # 398   W = 6.049e-01iteration # 399   W = 6.049e-01iteration # 400   W = 6.049e-01iteration # 401   W = 6.049e-01iteration # 402   W = 6.049e-01iteration # 403   W = 6.049e-01iteration # 404   W = 6.049e-01iteration # 405   W = 6.049e-01iteration # 406   W = 6.049e-01iteration # 407   W = 6.049e-01iteration # 408   W = 6.049e-01iteration # 409   W = 6.049e-01iteration # 410   W = 6.049e-01iteration # 411   W = 6.049e-01iteration # 412   W = 6.049e-01iteration # 413   W = 6.049e-01iteration # 414   W = 6.049e-01iteration # 415   W = 6.049e-01iteration # 416   W = 6.049e-01iteration # 417   W = 6.049e-01iteration # 418   W = 6.049e-01iteration # 419   W = 6.049e-01iteration # 420   W = 6.049e-01iteration # 421   W = 6.049e-01iteration # 422   W = 6.049e-01iteration # 423   W = 6.049e-01iteration # 424   W = 6.049e-01iteration # 425   W = 6.049e-01iteration # 426   W = 6.049e-01iteration # 427   W = 6.049e-01iteration # 428   W = 6.049e-01iteration # 429   W = 6.049e-01iteration # 430   W = 6.049e-01iteration # 431   W = 6.049e-01iteration # 432   W = 6.049e-01iteration # 433   W = 6.049e-01iteration # 434   W = 6.049e-01iteration # 435   W = 6.049e-01iteration # 436   W = 6.049e-01iteration # 437   W = 6.049e-01iteration # 438   W = 6.049e-01iteration # 439   W = 6.049e-01iteration # 440   W = 6.049e-01iteration # 441   W = 6.049e-01iteration # 442   W = 6.049e-01iteration # 443   W = 6.049e-01iteration # 444   W = 6.049e-01iteration # 445   W = 6.049e-01iteration # 446   W = 6.049e-01iteration # 447   W = 6.049e-01iteration # 448   W = 6.049e-01iteration # 449   W = 6.049e-01iteration # 450   W = 6.049e-01iteration # 451   W = 6.049e-01iteration # 452   W = 6.049e-01iteration # 453   W = 6.049e-01iteration # 454   W = 6.049e-01iteration # 455   W = 6.049e-01iteration # 456   W = 6.049e-01iteration # 457   W = 6.049e-01iteration # 458   W = 6.049e-01iteration # 459   W = 6.049e-01iteration # 460   W = 6.049e-01iteration # 461   W = 6.049e-01iteration # 462   W = 6.049e-01iteration # 463   W = 6.049e-01iteration # 464   W = 6.049e-01iteration # 465   W = 6.049e-01iteration # 466   W = 6.049e-01iteration # 467   W = 6.049e-01iteration # 468   W = 6.049e-01iteration # 469   W = 6.049e-01iteration # 470   W = 6.049e-01iteration # 471   W = 6.049e-01iteration # 472   W = 6.049e-01iteration # 473   W = 6.049e-01iteration # 474   W = 6.049e-01iteration # 475   W = 6.049e-01iteration # 476   W = 6.049e-01iteration # 477   W = 6.049e-01iteration # 478   W = 6.049e-01iteration # 479   W = 6.049e-01iteration # 480   W = 6.049e-01iteration # 481   W = 6.049e-01iteration # 482   W = 6.049e-01iteration # 483   W = 6.049e-01iteration # 484   W = 6.049e-01iteration # 485   W = 6.049e-01iteration # 486   W = 6.049e-01iteration # 487   W = 6.049e-01iteration # 488   W = 6.049e-01iteration # 489   W = 6.049e-01iteration # 490   W = 6.049e-01iteration # 491   W = 6.049e-01iteration # 492   W = 6.049e-01iteration # 493   W = 6.049e-01iteration # 494   W = 6.049e-01iteration # 495   W = 6.049e-01iteration # 496   W = 6.049e-01iteration # 497   W = 6.049e-01iteration # 498   W = 6.049e-01iteration # 499   W = 6.049e-01iteration # 500   W = 6.049e-01
Network training ended at 14. 9.50

Plot the structure of the NNARX model

figure, drawnet(W1,W2)

Extract the estimated parameters of the NNARX model

alfa=W2(1:end-1);
alfa_0=W2(end);
beta=W1(:,1:end-1);
beta_0=W1(:,end);

Compute the predicted output of the NNARX model using the validation dataset, defining the regressor vector phi as a column

% Step 6: computation of predicted output

t_min=max([na+1,nb+nk]);
for t=t_min:Nv,
    phi= [yv_0mean(t-1:-1:t-na); uv_0mean(t-nk:-1:t-nk-nb+1)];
    yp(t)=alfa*tanh(beta*phi+beta_0)+alfa_0;
end

figure, plot(1:Nv,yv_0mean,'g', 1:Nv,yp,'r-.')

Compute the NNARX performances in terms of RMSE and MDL

% Step 7: computation of RMSE and MDL

N0=10; na,nb,nk
RMSE=sqrt(1/(Nv-N0))*norm(yv_0mean(N0+1:end)-yp(N0+1:end)')
n=na+nb;
MDL=n*log(Nv-N0)/(Nv-N0)+log(RMSE^2)

na =
     1
nb =
     1
nk =
     1
RMSE =
    1.0166
MDL =
    0.0469

Choosing the best model

The RMSE is minimized over 10 trials for each model (with N0=10):