Problem with system identification toolbox example

Hello, i have a idnlarx model and i was trying to find a code that replicate the same computations done by the comand sim(model,inputdata,x0).

I found this example in the system identification toolbox user guide but i'm having some issues understanding it

%openExample('ident/SimulationAndPredictionOfSigmoidNetworkExample')

Simulation and Prediction of Sigmoid Network

This example shows how the software computes the simulated and predicted output of a nonlinear ARX model as a result of evaluating the output of its nonlinearity estimator for given regressor values.

Estimate nonlinear ARX model with sigmoid network nonlinearity.

load twotankdata
estData = iddata(y,u,0.2,'Tstart',0);
M = nlarx(estData,[1 1 0],'idSigmoidNetwork');

Explore Nonlinear ARX Model

Inspect the model properties and estimation result. This command provides information about input and output variables, regressors, and nonlinearity estimator.

present(M)

M =

Nonlinear ARX model with 1 output and 1 input
  Inputs: u1
  Outputs: y1

Regressors:
  Linear regressors in variables y1, u1

Output function: Sigmoid network with 10 units
Sample time: 0.2 seconds

Status:                                                      
Termination condition: Maximum number of iterations reached..
Number of iterations: 20, Number of function evaluations: 243
                                                             
Estimated using NLARX on time domain data "estData".         
Fit to estimation data: 96.31% (prediction focus)            
FPE: 4.804e-05, MSE: 4.666e-05                               
More information in model's "Report" property.

NL = M.OutputFcn;
class(NL) 

ans = 'idSigmoidNetwork'

NL

NL = 
Sigmoid Network
Inputs: y1(t-1), u1(t)
Output: y1(t)

 Nonlinear Function: Sigmoid network with 10 units
 Linear Function: initialized to [-0.161 -0.105]
 Output Offset: initialized to 0.00119

    NonlinearFcn: '<Sigmoid units and their parameters>'
       LinearFcn: '<Linear function parameters>'
          Offset: '<Offset parameters>'

Simulate Output

The model output is:

y1(t)=f(y1(t-1),u1(t))

where f is the sigmoid network function. The model regressors y1(t-1) and u1(t) are inputs to the nonlinearity estimator. Time t is a discrete variable representing kT , where k= 0, 1,.., and T is the sampling interval. In this example, T=0.2 second.

The simulated output is:

ys(t) = f(ys(t-1),u1_meas(t))

where ys(t) is the simulated value of the response at time t. The simulation equation is the same as the prediction equation, except that the past output value ys(t-1) results from the simulation at the previous time step, rather than the measured output value.

Computing the simulated response includes:

• Computing regressor values from input-output data using simulated output values.
• Evaluating the nonlinearity for given regressor values.

Specify zero initial states.The model has one state because there is only one delayed term y1(t-1). The number of states is equal to sum(getDelayInfo(M)).

x0 = 0;

Compute the simulated output at time t =0, ys(t=0)

RegValue = [0,estData.u(1)];
ys_0 = evaluate(NL,RegValue);

RegValue is the vector of regressors at t=0. ys(t=0)=f(y1(t=-1),u1_meas(t=0)). In terms of MATLAB variables, this output is f(0,estData.u(1)), where

• y1(t=-1) is the initial state x0 (=0).
• u1_meas(t=0) is the value of the input at t =0, which is the first input data sample estData.u

Then to do it for all the inputs i create a loop:

all_ys = zeros(length(estData.y), 1);
all_ys(1) = ys_0;

for t = 1:length(estData.y)-1

    RegValue_t = [all_ys(t),estData.u(t+1)];

    ys_t = evaluate(NL,RegValue_t);

    all_ys(t+1) = ys_t;
end

Now in the example it says that: "Unlike for output prediction, you cannot use getreg to compute regressor values for all time values. You must compute regressors values at each time sample separately because the output samples required for forming the regressor vector are available iteratively, one sample at a time.

These steps are equivalent to the simulated response computed in a single step using sim(idnlarx)."

ys = sim(M,estData,x0); %passing the whole iddata object or only the inputs doesn't change the results

My issue is that the two outputs are not the same.

plot(ys.InputData)
hold on
plot(all_ys)

I tried putting the optimization focus on simulation but the results are a similar, stil not comparable with the sim outputs. Am i interpreting wrong what is written? Thank you for your attention, hope the question is clear.

Answer 1

[Tianyu]

Hi Rocco,

The workflow you provided is correct, however, you missed the normalization step. The following script provides the desired result. Hope it helps.

ys = sim(M,estData,x0);

all_ys = zeros(length(estData.y), 1);
all_ys(1) = ys_0;
for t = 1:length(estData.y)-1
    RegValue_t = [all_ys(t),estData.u(t+1)];
    RegValue_t = (RegValue_t-M.Normalization.RegressorCenter)./M.Normalization.RegressorScale;
    ys_t = evaluate(NL,RegValue_t);
    ys_t = ys_t.*M.Normalization.OutputScale.'+M.Normalization.OutputCenter.';
    all_ys(t+1) = ys_t;
end
plot(ys.OutputData)
hold on
plot(all_ys)

[Rocco]

Thank you very much, it works!