[System Identification] Why NARX neural network and Hammerstein-Wiener model perform badly than simple sigmoid network nonlinearity estimator on any predictions?

I am currently working on dynamic modeling and exploring different techniques and algorithm to deploy a dynamic time-series black-box model.

My data looks like following: I have 7 inputs and 1 outputs (MISO system), the measurments contains delay and white noises already. I also separate the training and test/validation data into 80/20 from 1800 data.

I have also noticed that the order estimation for NARX modeling is very important. Thus I used arxstruc and selstruc functions to select model order, which are not very accurate. However, it is computational heavy to find any nonlinear regressors and narxnet cannot specify different model order in regard to multiple inputs, which makes the model prediction to suffer, even worst than normal neural network without timeseries.

Do you have any suggestions, how I can find accurate estimated model order, espically for nonlinear regressors? Are there any existing algorithm other than AIC and IV method, which I tried both of them?

I have tried using these kernels and methodologies to estimate the NARX model parameter. For validations I used compare function and NRMSE to display the goodness of the fit:

• Treepartition: around 53% NRMSE

Treepartition= nlarx(train_data,model_order, treepartition);

• Sigmoidnet: around 65% NRMSE

NL = sigmoidnet('NumberOfUnits',30);
opt = nlarxOptions('SearchMethod','gna');
opt.Regularization.Lambda = 5e-8;
Sigmoid = nlarx(train_data,model_order,NL,opt);

• NARX Feedforward Neural network, using nlarx function: around 45% NRMSE

setdemorandstream(491218381);
ff = feedforwardnet([10 20]);
ff.layers{2}.transferFcn = 'radbas';
ff.trainParam.epochs = 50;
ff.divideFcn = 'dividerand'  % Divide data randomly
ff.divideMode = 'sample';  % Divide up every sample
ff.divideParam.trainRatio = 0.7;
ff.divideParam.valRatio = 0.15;
ff.divideParam.testRatio = 0.15;
FFNN = nlarx(train_data,model_order, neuralnet(ff));

• Cascade Neural network: around 57% NRMSE

net = cascadeforwardnet([20]);
net_estimator = neuralnet(net);
Cascade_50 = nlarx(train_80, nn_50, net_estimator);

• Hammerstein-Wiener model: around 49% NRMSE

opt = nlhwOptions();
opt.SearchMethod = 'gna'; %Gauss-Newton search
opt.SearchOptions.MaxIterations = 20;
opt.Regularization.Lambda = 1e-8;
HWOptSigW = nlhw(train_data,model_order, sigmoidnet, wavenet, opt); % 51.7% NRMSE
HWOptSig = nlhw(train_data,model_order, sigmoidnet, sigmoidnet, opt); % 47% NRMSE

• Normal Neural network without timeseries: around 39% NRMSE

setdemorandstream(491218381);
trainFcn = 'trainlm';  % Levenberg-Marquardt backpropagation.
hiddenLayerSize = [10 20];
net_ANN = fitnet(hiddenLayerSize,trainFcn);
net_ANN.divideFcn = 'dividerand' %'divideind';  % Divide data randomly
net_ANN.divideMode = 'sample';  % Divide up every sample
net_ANN.divideParam.trainRatio = 70/100;
net_ANN.divideParam.valRatio = 15/100;
net_ANN.divideParam.testRatio = 15/100;
net_ANN.layers{2}.transferFcn = 'radbas';
net_ANN.performFcn = 'mse';
% Train the Network
[net_ANN,tr] = train(net_ANN,x,t);

• Timeseries NARX Neural network, using narxnet function: around 12% NRMSE

setdemorandstream(491218381);
trainFcn = 'trainlm';  % Levenberg-Marquardt backpropagation.
inputDelays = 1:2;
feedbackDelays = 1;
hiddenLayerSize = [10 20];
net_NARX = narxnet(inputDelays,feedbackDelays,hiddenLayerSize,'open',trainFcn);
net_NARX.divideFcn = 'dividerand';
net_NARX.divideMode = 'time'
net_NARX.divideParam.trainRatio = 0.7
net_NARX.divideParam.valRatio = 0.15
net_NARX.divideParam.testRatio = 0.15
net_NARX.layers{2}.transferFcn = 'radbas';
net_NARX.performFcn = 'mse'; 
% Train the Network
[net_NARX,tr] = train(net_NARX,x,t,xi,ai);

So in conclusion, the best NRMSE result is 65% from Sigmoidnet, it's R^2 is 0.88 and RMSE is 72.87. The timeseries NARX Neural network, using narxnet function performs miserably due to the fixed inputs orders. But other neural network models aren't perform better either. As well as Hammerstein-Wiener model although it contains sigmoid functions.

So here are my questions:

1. Do you have any suggestions on how I can improve the overall prediction accuracy of the dynamic models?
2. Why NARX neural network and Hammerstein-Wiener model perform so badly?
3. Do I need more data to train?
4. How I can find accurate estimated model order, espically for nonlinear regressors?
5. And do I have an option to deploy the matlab script into Simulink? As I realized that there were not much degree of freedom to configure the different inputs order using the system identification tool box and other neural network tool boxes.

Thank you for reading so patiently and I am looking forward to your kindly suggestions and helps :)