# Online State Estimation Using Identified Models - Nonlinear Models

This example describes how to generate the state-transition and measurement functions for online state and output estimation using nonlinear system identification techniques.

System Identification Toolbox™ provides several recursive algorithms for state estimation such as Kalman Filter, Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF), and Particle Filter (PF). Use of a state estimator requires specification of state and measurement dynamics:

• How the system advances the states, represented by the state-transition equation$\mathit{X}\left(\mathit{k}+1\right)=\mathit{f}\left(\mathit{X}\left(\mathit{k}\right),\dots \right)$, where $\mathit{X}\left(\mathit{k}\right)$ is the state value at the time step $\mathit{k}$; the function$\mathit{f}\left(\cdot \right)\text{\hspace{0.17em}}$can accept additional input arguments to accommodate any dependency on exogenous inputs and the time instant.

• How the system computes the output or measurement $\mathit{y}\left(\mathit{k}\right)$ as a function of $\mathit{X}\left(\mathit{k}\right)$ and possibly some additional inputs; $\mathit{y}\left(\mathit{k}\right)=\mathit{g}\left(\mathit{X}\left(\mathit{k}\right),\dots \right)$

In addition, you are required to provide the covariance values of the various disturbances affecting the state and measurement dynamics.

This example shows how you can generate the required state-transition and measurement equations by using nonlinear (black-box or grey-box) identification techniques. Four different modeling techniques are explored - Nonlinear ARX, Hammerstein-Wiener, Neural State Space and Nonlinear Grey Box.

This example describes:

• How to extract, and separate out, the state-transition and measurement functions from identified nonlinear models of various types.

• How to separate out noise dynamics in case the nonlinear model uses a nontrivial (that is, not an Output-Error structure) disturbance component; this is true only for the Nonlinear ARX models.

• How to extend the nonlinear models of Output-Error structure in order to prescribe the state/process noise covariance required by the state estimators such as EKF.

#### Definition of an Output-Error Model Structure

A generic discrete-time dynamic model in the System Identification Toolbox has the following representation:

`$\mathit{y}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{y}\left(\mathit{t}-1\right),\mathit{y}\left(\mathit{t}-2\right),...,\mathit{u}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}-1\right),....,\mathit{e}\left(\mathit{t}\right),\mathit{e}\left(\mathit{t}-1\right),...\right)$`

where the function $\mathit{f}\left(.\right)$ relates the past values of the output $\mathit{y}\left(\mathit{t}\right)$, the current/past values of the input $\mathit{u}\left(\mathit{t}\right)$, and the disturbances current/past values $\mathit{e}\left(\mathit{t}\right)$ to the current value of the output. In the linear case, this can be written as a stochastic relationship using additive noise:

`$\mathit{y}\left(\mathit{t}\right)=\mathit{G}\left(\mathit{q}\right)\mathit{u}\left(\mathit{t}\right)+\mathit{H}\left(\mathit{q}\right)\mathit{e}\left(\mathit{t}\right)$`

where $\mathit{G}\left(\mathit{q}\right)$, $\mathit{H}\left(\mathit{q}\right)$ are linear filters described by rational functions. When the noise filter $\mathit{H}\left(\mathit{q}\right)\equiv 1$, the model is called an Output-Error model. This is because when $\mathit{H}\left(\mathit{q}\right)$=1, the governing equation simplifies to:

`$\mathit{y}\left(\mathit{t}\right)=\mathit{G}\left(\mathit{q}\right)\mathit{u}\left(\mathit{t}\right)+\mathit{e}\left(\mathit{t}\right)$`

Here, the output at time $\mathit{t}$ is not affected by errors at previous time instants. In state-space form, this translates to a model of the form:

`$\begin{array}{l}\mathit{x}\left(\mathit{t}+1\right)=\mathrm{Ax}\left(\mathit{t}\right)+\mathit{B}\mathit{u}\left(\mathit{t}\right)\\ \mathit{y}\left(\mathit{t}\right)=\mathrm{Cx}\left(\mathit{t}\right)+\mathrm{Du}\left(\mathit{t}\right)+\mathit{e}\left(\mathit{t}\right)\end{array}$`

This is a system with no process noise. A state estimator based on this model cannot correct the state predictions using current output measurements. Hence an Output-Error model needs to extended in some way to extract the process noise information.

Most nonlinear models in the System Identification Toolbox are of Output-Error structure. The only exception is Nonlinear ARX models, where the model structure is:

`$\mathit{y}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{y}\left(\mathit{t}-1\right),\mathit{y}\left(\mathit{t}-2\right),...,\mathit{u}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}-1\right)\right)+\mathit{e}\left(\mathit{t}\right)$`

where $\mathit{y}\left(\mathit{t}\right)$ are measured (noise-corrupted) outputs values. The influence of past errors $\mathit{e}\left(\mathit{t}-1\right),\mathit{e}\left(\mathit{t}-2\right),...$can be seen more clearly by writing the equation in terms of the true (noise free) output ${\mathit{y}}_{0}\left(\mathit{t}\right)$, so that $\mathit{y}\left(\mathit{t}\right)={\mathit{y}}_{0}\left(\mathit{t}\right)+\mathit{e}\left(\mathit{t}\right)$:

`${\mathit{y}}_{0}\left(\mathit{t}\right)=\mathit{f}\left({\mathit{y}}_{0}\left(\mathit{t}-1\right)+\mathit{e}\left(\mathit{t}-1\right),{\mathit{y}}_{0}\left(\mathit{t}-2\right)+\mathit{e}\left(\mathit{t}-2\right),...,\mathit{u}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}-1\right)\right)$`

### Example System: An Industrial Robot Arm

Consider a robot arm is described by a nonlinear three-mass flexible model. The input to the robot is the applied torque $u\left(t\right)=\tau \left(t\right)$ generated by the electrical motor, and the resulting angular velocity of the motor $y\left(t\right)=\frac{d}{dt}{q}_{m}\left(t\right)$ is the measured output.

It is possible to derive the equations of motion for this system. For example a nonlinear state-space description of this system uses 5 states. For more information on this system, see the example titled "`Modeling an Industrial Robot Arm`".

This system is excited with inputs of different profiles and the resulting angular velocity of the motor recorded. A periodic excitation signal with approximately 10 seconds duration was employed to generate a reference speed for the controller. The chosen sampling frequency was `2 kHz` (sample time, `Ts `= 0.0005 seconds). For estimation (training) data, a multisine input signal with a flat amplitude spectrum in the frequency interval 1-40 Hz with a peak value of 16 rad/s was used. The multisine signal is superimposed on a filtered square wave with amplitude 20 rad/s and cut-off frequency 1 Hz. For validation data, a multisine input signal containing frequencies 0.1, 0.3, and 0.5 Hz was used, with peak value 40 rad/s.

Both datasets are downsampled 10 times for the purposes of the current example.

```load robotarmdata.mat % Prepare estimation data eData = iddata(ye,ue,5e-4,'InputName','Torque','OutputName','Angular Velocity','Tstart',0); eData = idresamp(eData,[1 10]); eData.Name = 'estimation data'```
```eData = Time domain data set with 1984 samples. Sample time: 0.005 seconds Name: estimation data Outputs Unit (if specified) Angular Velocity Inputs Unit (if specified) Torque Data Properties ```
```% Prepare validation data vData = iddata(yv3,uv3,5e-4,'InputName','Torque','OutputName','Angular Velocity','Tstart',0); vData = idresamp(vData,[1 10]); vData.Name = 'validation data'; idplot(eData, vData) legend show```

### Nonlinear Black Box Model

Suppose you have prior knowledge that the underlying process is nonlinear (or that the linear black box model does not provide a good result), but you do not have detailed knowledge of the underlying physics. In that case you can identify nonlinear black-box models of different structures. Then use the identified models to build up the state estimators.

#### Nonlinear ARX Model

A nonlinear ARX model is a discrete-time model that has the following structure:

`$\mathit{y}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{y}\left(\mathit{t}-1\right),\mathit{y}\left(\mathit{t}-2\right),...,\mathit{u}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}-1\right)\right)+\mathit{e}\left(\mathit{t}\right)$`

That is, the output at time $t$ is computed as a nonlinear function of past outputs and present and past inputs. At any time step the disturbance $e\left(t\right)$ gets added to the output. All future values $\mathit{y}\left(\mathit{t}+1\right),\mathit{y}\left(\mathit{t}+2\right),\dots$ would typically contain the effect of this disturbance since the past outputs are used to compute the future ones. In a sense, the effect of disturbances accumulate as time progresses. As a result, long-term predictions are usually harder with these models especially if they have been trained in open loop (for example, by using `Focus = 'prediction'` with the `nlarx` command).

#### IDNLARX Model States

The states of this model are defined as the lagged input/output variables. So a nonlinear ARX model using the equation:

`$\mathit{y}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{y}\left(\mathit{t}-1\right),\text{\hspace{0.17em}}\mathit{u}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}-2\right)\right)+\mathit{e}\left(\mathit{t}\right)$`

can be expressed in state-space form by using the state variables:

`$\begin{array}{l}{\mathit{x}}_{1}\left(\mathit{t}\right)=\mathit{y}\left(\mathit{t}-1\right)\\ {\mathit{x}}_{2}\left(\mathit{t}\right)=\mathit{u}\left(\mathit{t}-1\right)\\ {\mathit{x}}_{3}\left(\mathit{t}\right)=\mathit{u}\left(\mathit{t}-2\right)\end{array}$`

leading to the following state-space representation:

`${\mathit{x}}_{1}\left(\mathit{t}+1\right)=\mathit{f}\left({\mathit{x}}_{1}\left(\mathit{t}\right),{\mathit{x}}_{3}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}\right)\right)+\mathit{e}\left(\mathit{t}\right)$`

`${\mathit{x}}_{2}\left(\mathit{t}+1\right)=\mathit{u}\left(\mathit{t}\right)$`

`${\mathit{x}}_{3}\left(\mathit{t}+1\right)={\mathit{x}}_{2}\left(\mathit{t}\right)$`

`$\mathit{y}\left(\mathit{t}\right)={\mathit{x}}_{1}\left(\mathit{t}+1\right)=\mathit{f}\left({\mathit{x}}_{1}\left(\mathit{t}\right),{\mathit{x}}_{3}\left(\mathit{t}\right),\mathit{u}\left(\mathit{t}\right)\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\mathit{e}\left(\mathit{t}\right)$`

Define $X\left(t\right)=\left[{x}_{1}\left(t\right),{x}_{2}\left(t\right),{x}_{3}\left(t\right){\right]}^{T}$. Then in the matrix form you have the following representation:

`$\begin{array}{l}X\left(t+1\right)=F\left(X\left(t\right),u\left(t\right)\right)+Ke\left(t\right)\\ y\left(t\right)=f\left(X\left(t\right),u\left(t\right)\right)+e\left(t\right)\end{array}$`

where $K=\left[1,0,0{\right]}^{T}$ is a constant matrix, and the vector function $\mathit{F}\left(.\right)={\left[\mathit{f}\left(.\right),\text{\hspace{0.17em}}\mathit{u}\left(\mathit{t}\right),{\text{\hspace{0.17em}}\mathit{x}}_{2}\left(\mathit{t}\right)\right]}^{\mathit{T}}$. Note that in this form, the process and measurement noises are correlated.

#### Identification of Wavelet Network Based Nonlinear ARX Model

Identify a Nonlinear ARX model for the robot arm data.

```vars = {'Angular Velocity','Torque'}; L = linearRegressor(vars,1:2:8); P = polynomialRegressor(vars,1:2:8,3); opt = nlarxOptions(SearchMethod='gna',Focus='simulation'); nlarx1 = nlarx(eData, [L;P], idWaveletNetwork(3), opt)```
```nlarx1 = Nonlinear ARX model with 1 output and 1 input Inputs: Torque Outputs: Angular Velocity Regressors: 1. Linear regressors in variables Angular Velocity, Torque 2. Order 3 regressors in variables Angular Velocity, Torque List of all regressors Output function: Wavelet network with 3 units Sample time: 0.005 seconds Status: Termination condition: No improvement along the search direction with line search.. Number of iterations: 6, Number of function evaluations: 129 Estimated using NLARX on time domain data "estimation data". Fit to estimation data: 78.18% (simulation focus) FPE: 6.334, MSE: 20.65 More information in model's "Report" property. Model Properties ```
```% Perform 10 step ahead prediction for the output data in the validation dataset. compare(vData, nlarx1, 10);```

#### Using the Model in the Extended Kalman Filter (EKF) Block

In order to use this model in the EKF block, you need to specify the state-transition and measurement functions as Simulink functions. For the current example, these are implemented using MATLAB Function Blocks. The `nlarxPredictFcn` function performs the core computation; this function is called by the MATLAB Function blocks implementing the state transition and measurement updates.

`type nlarxPredictFcn`
```function [y,xnext] = nlarxPredictFcn(x,u) % Predict the states and output of a nonlinear ARX model. % The identified model ("nlarx1") is assumed to be available in the MATLAB base workspace. % Copyright 2022 The MathWorks, Inc. sys = evalin('base','nlarx1'); % Step 1: Compute model regressors as a function of states. Use GETERG to see the list of % regressors. Regressors are non-tunable, static functions of model states and inputs. I = [1 3 5 7 8 10 12 14]; R = [x(I,:); x(I,:).^3].'; % Step 2: Map the regressors to the model output. This is required since x1(t+1) = y(t) for this % model. Note that output y is a static function (represented by the wavelet network) of the model's % regressors (R). N = sys.Normalization; NL = sys.OutputFcn; R = (R-N.RegressorCenter)./N.RegressorScale; y = evaluate(NL, R); y = (y.*N.OutputScale + N.OutputCenter).'; if nargout>1 % Compose the output xnext := x(t+1) xnext = x; % initialization xnext(1,:) = y; xnext(2:7,:) = x(1:6,:); xnext(8,:) = u; xnext(9:14,:) = xnext(8:13,:); end ```

The state estimation scenario based on the identified model (`nlarx1`) is implemented in the Simulink model `StateEstimation_Using_NonlinearARX.slx`. The scenario considered is 10-step ahead prediction realized by enabling the multi-rate operation; the measurements are used for correction of the state estimates only once for every 10 successive predictions.

```% Prepare some data for the EKF Block NV = nlarx1.NoiseVariance; % measurement noise variance K = zeros(14,1); % (there are 14 states) K(1) = 1; StateNoiseCovariance = K*NV*K'; Ts = nlarx1.Ts; % prediction time interval Ts2 = 10*Ts; % correction time interval open_system('StateEstimation_Using_NonlinearARX')```

#### Using the Model in the Unscented Kalman Filter (UKF) Block

The UKF block is configured in the same manner as the EKF block. The Simulink functions `stateTransitionFcn` and `measurementFcn` in the model `StateEstimation_Using_NonlinearARX.slx` are used by both the EKF and UKF blocks.

#### Using the Model in the Particle Filter (PF) Block

For the Particle Filter block, you need to specify the transition function for the particles and the value of the measurement likelihood as a probability distribution.

`NumParticles = 1000; % number of particles used by the Particle filter`

#### 10 Step Ahead State Estimation Results

`sim('StateEstimation_Using_NonlinearARX');`

For this example, the three state estimators perform similarly. Note that the system output is same as the first state value. The filters underestimate the outputs at certain regions of rapid change. This observation is consistent with the direct 10-steap ahead prediction result shown in the `compare` plot.

#### Hammerstein Wiener Model

Another popular type of nonlinear dynamic model is a Hammerstein-Wiener model. It connects static nonlinear functions in series with a linear transfer function. The states of this model are the states of its linear block.

```% Estimate a model that uses piecewise linear functions for its input and output nonlinearities nlhw1 = nlhw(eData,[7 7 1])```
```nlhw1 = Hammerstein-Wiener model with 1 output and 1 input Linear transfer function corresponding to the orders nb = 7, nf = 7, nk = 1 Input nonlinearity: Piecewise linear with 10 break-points Output nonlinearity: Piecewise linear with 10 break-points Sample time: 0.005 seconds Status: Termination condition: Maximum number of iterations reached.. Number of iterations: 20, Number of function evaluations: 116 Estimated using NLHW on time domain data "estimation data". Fit to estimation data: 89.91% FPE: 4.674, MSE: 4.413 More information in model's "Report" property. Model Properties ```
`compare(vData, nlhw1) % validate the model quality on an independent dataset`

A Hammerstein Wiener model is an output-error model; there is no process noise. However, one can consider extending the identified model by prescribing process noise of non-zero covariance (covariance value determination would require a separate analysis). For this example, this extension is not explored.

```uNL = nlhw1.InputNonlinearity; yNL = nlhw1.OutputNonlinearity; linTF = nlhw1.LinearModel; [A,B,C,D] = ssdata(linTF); Ts = nlhw1.Ts;```

For using this model for state estimation, you can use the intermediate signal `w(t)` as the input signal since the I/O nonlinearities do not contain any states.

```t = vData.SamplingInstants; % the time vector for the validation data u = vData.InputData; % the input signal u(t) w = evaluate(uNL, u); % the output of the "Input Nonlinearity f" block w(t)```

The model states can be obtained by simulating the linear block using the signal w(t) and any prescribed initial conditions. For example:

```[~,~,x] = sim(nlhw1, w); plot(t,x) xlabel('Time (seconds)') ylabel('States (X)')```

shows the state trajectories generated for zero initial conditions.

Since there is no process noise, there is no correction step required in the state predictors. The entire state prediction exercise is same as a pure simulation. The state estimators such as EKF are required only when there is an opportunity to correct the estimates using the measurements. This is explored more in the context of state-space models.

### Nonlinear State Space Models

#### Neural State Space Model

A neural state space model is a continuous- or discrete-time nonlinear model where you specify the state transition and measurement functions as multi-layer feedforward neural networks. This modeling technique requires Deep Learning Toolbox™. Neural state space models are output-error models - they assume the states are not affected by process noise.

Training a neural state space model requires measurements of the states as well as the outputs. Assume that the state definitions are not available. Hence define states in terms of lagged inputs and outputs. You can generate data for the lagged variables by using the `getreg` command that is called on a template Nonlinear ARX model that uses the regressors defined as the lagged variables.

```eData2 = idresamp(eData,[1 2]); % reduce data for quicker results % normalize the data eM = mean([eData2.OutputData,eData2.InputData]); eSD = std([eData2.OutputData,eData2.InputData]); eData2.OutputData = (eData2.OutputData-eM(1))/eSD(1); eData2.InputData = (eData2.InputData-eM(2))/eSD(2); Ts = eData2.Ts; vars = {'Angular Velocity','Torque'}; nx = 7; L = linearRegressor(vars,1:nx); n0 = idnlarx(vars(1),vars(2),L,'Ts',Ts); % View the state variables States = getreg(n0); % there are 8 states disp(strjoin(States,'\n'))```
```Angular Velocity(t-1) Angular Velocity(t-2) Angular Velocity(t-3) Angular Velocity(t-4) Angular Velocity(t-5) Angular Velocity(t-6) Angular Velocity(t-7) Torque(t-1) Torque(t-2) Torque(t-3) Torque(t-4) Torque(t-5) Torque(t-6) Torque(t-7) ```
```% get training data corresponding to these states XData = getreg(n0, eData2); % convert data into an iddata object; add the actual output (Angular Velocity(t)) as the last % output signal ymeas = eData2.OutputData; xmeas = XData{:,:}; XYData = iddata([xmeas,ymeas],eData2.InputData,... 'OutputName',[States; vars(1)],'InputName',vars(2),'Ts',Ts,'Tstart',0); % Discard first 4 samples that contain effects of unknown initial conditions XYData = XYData(nx+1:end);```

Create a Neural State Space model that uses 14 states and one input. Since the states are considered measured, the model has 15 outputs, one corresponding to each state, plus the actual output.

`nss = idNeuralStateSpace(numel(States),NumInputs=1,NumOutputs=numel(States)+1,Ts=Ts);`

Model the state-transition function as a 1 hidden layer network using `sigmoid` activation layer.

```rng(5) % for reproducibility nss.StateNetwork = createMLPNetwork(nss,'state', ... LayerSizes=numel(States), ... Activations="sigmoid", ... WeightsInitializer="glorot", ... BiasInitializer="zeros"); summary(nss.StateNetwork)```
``` Initialized: true Number of learnables: 434 Inputs: 1 'x[k]' 14 features 2 'u[k]' 1 features ```

The first 14 outputs are the states themselves Hence nss.OutputNetwork(1) is a trivial network representing the equation ${\mathit{y}}_{1}\left(\mathit{t}\right)=\mathit{x}\left(\mathit{t}\right)$. The 15th output of this model represents the true measured output ${\mathit{y}}_{2}\left(\mathit{t}\right)=\mathit{y}\left(\mathit{t}\right)$. For this output, create a separate, single hidden layer feedforward network using tanh activation layer.

```nss.OutputNetwork(2) = createMLPNetwork(nss,'output', ... LayerSizes=1, ... Activations="sigmoid", ... WeightsInitializer="ones", ... BiasInitializer="zeros"); summary(nss.OutputNetwork(2))```
``` Initialized: true Number of learnables: 17 Inputs: 1 'x[k]' 14 features 2 'u[k]' 1 features ```

Prepare training data by splitting it into 4 overlapping segments.

```NumBatches = 4; Ns = size(XYData,1); Nsb = 300; Data = cell(1,NumBatches); pt = 0; for i = 1:NumBatches if pt+Nsb<=Ns Data{i} = XYData(pt+(1:Nsb)); else Data{i} = XYData(pt+1:end); end Data{i}.Tstart = 0; pt = pt + Nsb-20; % 20 sample overlap end Data = merge(Data{:});```

Prepare training options. Choose ADAM Solver and increase the learning rate to 0.002. Also set other options such as maximum number of epochs to run and the frequency of validation plot update.

```opt = nssTrainingOptions('adam'); opt.MaxEpochs = 1500; opt.MiniBatchSize = Nsb; opt.LearnRate = 0.02;```

Train the model, where you use the last experiment in Data for in-estimation cross validation.

```tic nss = nlssest(Data,nss,opt,UseLastExperimentForValidation=true)```
```Training in progress (completed epoch/max epochs): 0/ 1500 1/ 1500 2/ 1500 3/ 1500 4/ 1500 5/ 1500 6/ 1500 7/ 1500 8/ 1500 9/ 1500 10/ 1500 11/ 1500 12/ 1500 13/ 1500 14/ 1500 15/ 1500 16/ 1500 17/ 1500 18/ 1500 19/ 1500 20/ 1500 21/ 1500 22/ 1500 23/ 1500 24/ 1500 25/ 1500 26/ 1500 27/ 1500 28/ 1500 29/ 1500 30/ 1500 31/ 1500 32/ 1500 33/ 1500 34/ 1500 35/ 1500 36/ 1500 37/ 1500 38/ 1500 39/ 1500 40/ 1500 41/ 1500 42/ 1500 43/ 1500 44/ 1500 45/ 1500 46/ 1500 47/ 1500 48/ 1500 49/ 1500 50/ 1500 51/ 1500 52/ 1500 53/ 1500 54/ 1500 55/ 1500 56/ 1500 57/ 1500 58/ 1500 59/ 1500 60/ 1500 61/ 1500 62/ 1500 63/ 1500 64/ 1500 65/ 1500 66/ 1500 67/ 1500 68/ 1500 69/ 1500 70/ 1500 71/ 1500 72/ 1500 73/ 1500 74/ 1500 75/ 1500 76/ 1500 77/ 1500 78/ 1500 79/ 1500 80/ 1500 81/ 1500 82/ 1500 83/ 1500 84/ 1500 85/ 1500 86/ 1500 87/ 1500 88/ 1500 89/ 1500 90/ 1500 91/ 1500 92/ 1500 93/ 1500 94/ 1500 95/ 1500 96/ 1500 97/ 1500 98/ 1500 99/ 1500 100/ 1500 101/ 1500 102/ 1500 103/ 1500 104/ 1500 105/ 1500 106/ 1500 107/ 1500 108/ 1500 109/ 1500 110/ 1500 111/ 1500 112/ 1500 113/ 1500 114/ 1500 115/ 1500 116/ 1500 117/ 1500 118/ 1500 119/ 1500 120/ 1500 121/ 1500 122/ 1500 123/ 1500 124/ 1500 125/ 1500 126/ 1500 127/ 1500 128/ 1500 129/ 1500 130/ 1500 131/ 1500 132/ 1500 133/ 1500 134/ 1500 135/ 1500 136/ 1500 137/ 1500 138/ 1500 139/ 1500 140/ 1500 141/ 1500 142/ 1500 143/ 1500 144/ 1500 145/ 1500 146/ 1500 147/ 1500 148/ 1500 149/ 1500 150/ 1500 151/ 1500 152/ 1500 153/ 1500 154/ 1500 155/ 1500 156/ 1500 157/ 1500 158/ 1500 159/ 1500 160/ 1500 161/ 1500 162/ 1500 163/ 1500 164/ 1500 165/ 1500 166/ 1500 167/ 1500 168/ 1500 169/ 1500 170/ 1500 171/ 1500 172/ 1500 173/ 1500 174/ 1500 175/ 1500 176/ 1500 177/ 1500 178/ 1500 179/ 1500 180/ 1500 181/ 1500 182/ 1500 183/ 1500 184/ 1500 185/ 1500 186/ 1500 187/ 1500 188/ 1500 189/ 1500 190/ 1500 191/ 1500 192/ 1500 193/ 1500 194/ 1500 195/ 1500 196/ 1500 197/ 1500 198/ 1500 199/ 1500 200/ 1500 201/ 1500 202/ 1500 203/ 1500 204/ 1500 205/ 1500 206/ 1500 207/ 1500 208/ 1500 209/ 1500 210/ 1500 211/ 1500 212/ 1500 213/ 1500 214/ 1500 215/ 1500 216/ 1500 217/ 1500 218/ 1500 219/ 1500 220/ 1500 221/ 1500 222/ 1500 223/ 1500 224/ 1500 225/ 1500 226/ 1500 227/ 1500 228/ 1500 229/ 1500 230/ 1500 231/ 1500 232/ 1500 233/ 1500 234/ 1500 235/ 1500 236/ 1500 237/ 1500 238/ 1500 239/ 1500 240/ 1500 241/ 1500 242/ 1500 243/ 1500 244/ 1500 245/ 1500 246/ 1500 247/ 1500 248/ 1500 249/ 1500 250/ 1500 251/ 1500 252/ 1500 253/ 1500 254/ 1500 255/ 1500 256/ 1500 257/ 1500 258/ 1500 259/ 1500 260/ 1500 261/ 1500 262/ 1500 263/ 1500 264/ 1500 265/ 1500 266/ 1500 267/ 1500 268/ 1500 269/ 1500 270/ 1500 271/ 1500 272/ 1500 273/ 1500 274/ 1500 275/ 1500 276/ 1500 277/ 1500 278/ 1500 279/ 1500 280/ 1500 281/ 1500 282/ 1500 283/ 1500 284/ 1500 285/ 1500 286/ 1500 287/ 1500 288/ 1500 289/ 1500 290/ 1500 291/ 1500 292/ 1500 293/ 1500 294/ 1500 295/ 1500 296/ 1500 297/ 1500 298/ 1500 299/ 1500 300/ 1500 301/ 1500 302/ 1500 303/ 1500 304/ 1500 305/ 1500 306/ 1500 307/ 1500 308/ 1500 309/ 1500 310/ 1500 311/ 1500 312/ 1500 313/ 1500 314/ 1500 315/ 1500 316/ 1500 317/ 1500 318/ 1500 319/ 1500 320/ 1500 321/ 1500 322/ 1500 323/ 1500 324/ 1500 325/ 1500 326/ 1500 327/ 1500 328/ 1500 329/ 1500 330/ 1500 331/ 1500 332/ 1500 333/ 1500 334/ 1500 335/ 1500 336/ 1500 337/ 1500 338/ 1500 339/ 1500 340/ 1500 341/ 1500 342/ 1500 343/ 1500 344/ 1500 345/ 1500 346/ 1500 347/ 1500 348/ 1500 349/ 1500 350/ 1500 351/ 1500 352/ 1500 353/ 1500 354/ 1500 355/ 1500 356/ 1500 357/ 1500 358/ 1500 359/ 1500 360/ 1500 361/ 1500 362/ 1500 363/ 1500 364/ 1500 365/ 1500 366/ 1500 367/ 1500 368/ 1500 369/ 1500 370/ 1500 371/ 1500 372/ 1500 373/ 1500 374/ 1500 375/ 1500 376/ 1500 377/ 1500 378/ 1500 379/ 1500 380/ 1500 381/ 1500 382/ 1500 383/ 1500 384/ 1500 385/ 1500 386/ 1500 387/ 1500 388/ 1500 389/ 1500 390/ 1500 391/ 1500 392/ 1500 393/ 1500 394/ 1500 395/ 1500 396/ 1500 397/ 1500 398/ 1500 399/ 1500 400/ 1500 401/ 1500 402/ 1500 403/ 1500 404/ 1500 405/ 1500 406/ 1500 407/ 1500 408/ 1500 409/ 1500 410/ 1500 411/ 1500 412/ 1500 413/ 1500 414/ 1500 415/ 1500 416/ 1500 417/ 1500 418/ 1500 419/ 1500 420/ 1500 421/ 1500 422/ 1500 423/ 1500 424/ 1500 425/ 1500 426/ 1500 427/ 1500 428/ 1500 429/ 1500 430/ 1500 431/ 1500 432/ 1500 433/ 1500 434/ 1500 435/ 1500 436/ 1500 437/ 1500 438/ 1500 439/ 1500 440/ 1500 441/ 1500 442/ 1500 443/ 1500 444/ 1500 445/ 1500 446/ 1500 447/ 1500 448/ 1500 449/ 1500 450/ 1500 451/ 1500 452/ 1500 453/ 1500 454/ 1500 455/ 1500 456/ 1500 457/ 1500 458/ 1500 459/ 1500 460/ 1500 461/ 1500 462/ 1500 463/ 1500 464/ 1500 465/ 1500 466/ 1500 467/ 1500 468/ 1500 469/ 1500 470/ 1500 471/ 1500 472/ 1500 473/ 1500 474/ 1500 475/ 1500 476/ 1500 477/ 1500 478/ 1500 479/ 1500 480/ 1500 481/ 1500 482/ 1500 483/ 1500 484/ 1500 485/ 1500 486/ 1500 487/ 1500 488/ 1500 489/ 1500 490/ 1500 491/ 1500 492/ 1500 493/ 1500 494/ 1500 495/ 1500 496/ 1500 497/ 1500 498/ 1500 499/ 1500 500/ 1500 501/ 1500 502/ 1500 503/ 1500 504/ 1500 505/ 1500 506/ 1500 507/ 1500 508/ 1500 509/ 1500 510/ 1500 511/ 1500 512/ 1500 513/ 1500 514/ 1500 515/ 1500 516/ 1500 517/ 1500 518/ 1500 519/ 1500 520/ 1500 521/ 1500 522/ 1500 523/ 1500 524/ 1500 525/ 1500 526/ 1500 527/ 1500 528/ 1500 529/ 1500 530/ 1500 531/ 1500 532/ 1500 533/ 1500 534/ 1500 535/ 1500 536/ 1500 537/ 1500 538/ 1500 539/ 1500 540/ 1500 541/ 1500 542/ 1500 543/ 1500 544/ 1500 545/ 1500 546/ 1500 547/ 1500 548/ 1500 549/ 1500 550/ 1500 551/ 1500 552/ 1500 553/ 1500 554/ 1500 555/ 1500 556/ 1500 557/ 1500 558/ 1500 559/ 1500 560/ 1500 561/ 1500 562/ 1500 563/ 1500 564/ 1500 565/ 1500 566/ 1500 567/ 1500 568/ 1500 569/ 1500 570/ 1500 571/ 1500 572/ 1500 573/ 1500 574/ 1500 575/ 1500 576/ 1500 577/ 1500 578/ 1500 579/ 1500 580/ 1500 581/ 1500 582/ 1500 583/ 1500 584/ 1500 585/ 1500 586/ 1500 587/ 1500 588/ 1500 589/ 1500 590/ 1500 591/ 1500 592/ 1500 593/ 1500 594/ 1500 595/ 1500 596/ 1500 597/ 1500 598/ 1500 599/ 1500 600/ 1500 601/ 1500 602/ 1500 603/ 1500 604/ 1500 605/ 1500 606/ 1500 607/ 1500 608/ 1500 609/ 1500 610/ 1500 611/ 1500 612/ 1500 613/ 1500 614/ 1500 615/ 1500 616/ 1500 617/ 1500 618/ 1500 619/ 1500 620/ 1500 621/ 1500 622/ 1500 623/ 1500 624/ 1500 625/ 1500 626/ 1500 627/ 1500 628/ 1500 629/ 1500 630/ 1500 631/ 1500 632/ 1500 633/ 1500 634/ 1500 635/ 1500 636/ 1500 637/ 1500 638/ 1500 639/ 1500 640/ 1500 641/ 1500 642/ 1500 643/ 1500 644/ 1500 645/ 1500 646/ 1500 647/ 1500 648/ 1500 649/ 1500 650/ 1500 651/ 1500 652/ 1500 653/ 1500 654/ 1500 655/ 1500 656/ 1500 657/ 1500 658/ 1500 659/ 1500 660/ 1500 661/ 1500 662/ 1500 663/ 1500 664/ 1500 665/ 1500 666/ 1500 667/ 1500 668/ 1500 669/ 1500 670/ 1500 671/ 1500 672/ 1500 673/ 1500 674/ 1500 675/ 1500 676/ 1500 677/ 1500 678/ 1500 679/ 1500 680/ 1500 681/ 1500 682/ 1500 683/ 1500 684/ 1500 685/ 1500 686/ 1500 687/ 1500 688/ 1500 689/ 1500 690/ 1500 691/ 1500 692/ 1500 693/ 1500 694/ 1500 695/ 1500 696/ 1500 697/ 1500 698/ 1500 699/ 1500 700/ 1500 701/ 1500 702/ 1500 703/ 1500 704/ 1500 705/ 1500 706/ 1500 707/ 1500 708/ 1500 709/ 1500 710/ 1500 711/ 1500 712/ 1500 713/ 1500 714/ 1500 715/ 1500 716/ 1500 717/ 1500 718/ 1500 719/ 1500 720/ 1500 721/ 1500 722/ 1500 723/ 1500 724/ 1500 725/ 1500 726/ 1500 727/ 1500 728/ 1500 729/ 1500 730/ 1500 731/ 1500 732/ 1500 733/ 1500 734/ 1500 735/ 1500 736/ 1500 737/ 1500 738/ 1500 739/ 1500 740/ 1500 741/ 1500 742/ 1500 743/ 1500 744/ 1500 745/ 1500 746/ 1500 747/ 1500 748/ 1500 749/ 1500 750/ 1500 751/ 1500 752/ 1500 753/ 1500 754/ 1500 755/ 1500 756/ 1500 757/ 1500 758/ 1500 759/ 1500 760/ 1500 761/ 1500 762/ 1500 763/ 1500 764/ 1500 765/ 1500 766/ 1500 767/ 1500 768/ 1500 769/ 1500 770/ 1500 771/ 1500 772/ 1500 773/ 1500 774/ 1500 775/ 1500 776/ 1500 777/ 1500 778/ 1500 779/ 1500 780/ 1500 781/ 1500 782/ 1500 783/ 1500 784/ 1500 785/ 1500 786/ 1500 787/ 1500 788/ 1500 789/ 1500 790/ 1500 791/ 1500 792/ 1500 793/ 1500 794/ 1500 795/ 1500 796/ 1500 797/ 1500 798/ 1500 799/ 1500 800/ 1500 801/ 1500 802/ 1500 803/ 1500 804/ 1500 805/ 1500 806/ 1500 807/ 1500 808/ 1500 809/ 1500 810/ 1500 811/ 1500 812/ 1500 813/ 1500 814/ 1500 815/ 1500 816/ 1500 817/ 1500 818/ 1500 819/ 1500 820/ 1500 821/ 1500 822/ 1500 823/ 1500 824/ 1500 825/ 1500 826/ 1500 827/ 1500 828/ 1500 829/ 1500 830/ 1500 831/ 1500 832/ 1500 833/ 1500 834/ 1500 835/ 1500 836/ 1500 837/ 1500 838/ 1500 839/ 1500 840/ 1500 841/ 1500 842/ 1500 843/ 1500 844/ 1500 845/ 1500 846/ 1500 847/ 1500 848/ 1500 849/ 1500 850/ 1500 851/ 1500 852/ 1500 853/ 1500 854/ 1500 855/ 1500 856/ 1500 857/ 1500 858/ 1500 859/ 1500 860/ 1500 861/ 1500 862/ 1500 863/ 1500 864/ 1500 865/ 1500 866/ 1500 867/ 1500 868/ 1500 869/ 1500 870/ 1500 871/ 1500 872/ 1500 873/ 1500 874/ 1500 875/ 1500 876/ 1500 877/ 1500 878/ 1500 879/ 1500 880/ 1500 881/ 1500 882/ 1500 883/ 1500 884/ 1500 885/ 1500 886/ 1500 887/ 1500 888/ 1500 889/ 1500 890/ 1500 891/ 1500 892/ 1500 893/ 1500 894/ 1500 895/ 1500 896/ 1500 897/ 1500 898/ 1500 899/ 1500 900/ 1500 901/ 1500 902/ 1500 903/ 1500 904/ 1500 905/ 1500 906/ 1500 907/ 1500 908/ 1500 909/ 1500 910/ 1500 911/ 1500 912/ 1500 913/ 1500 914/ 1500 915/ 1500 916/ 1500 917/ 1500 918/ 1500 919/ 1500 920/ 1500 921/ 1500 922/ 1500 923/ 1500 924/ 1500 925/ 1500 926/ 1500 927/ 1500 928/ 1500 929/ 1500 930/ 1500 931/ 1500 932/ 1500 933/ 1500 934/ 1500 935/ 1500 936/ 1500 937/ 1500 938/ 1500 939/ 1500 940/ 1500 941/ 1500 942/ 1500 943/ 1500 944/ 1500 945/ 1500 946/ 1500 947/ 1500 948/ 1500 949/ 1500 950/ 1500 951/ 1500 952/ 1500 953/ 1500 954/ 1500 955/ 1500 956/ 1500 957/ 1500 958/ 1500 959/ 1500 960/ 1500 961/ 1500 962/ 1500 963/ 1500 964/ 1500 965/ 1500 966/ 1500 967/ 1500 968/ 1500 969/ 1500 970/ 1500 971/ 1500 972/ 1500 973/ 1500 974/ 1500 975/ 1500 976/ 1500 977/ 1500 978/ 1500 979/ 1500 980/ 1500 981/ 1500 982/ 1500 983/ 1500 984/ 1500 985/ 1500 986/ 1500 987/ 1500 988/ 1500 989/ 1500 990/ 1500 991/ 1500 992/ 1500 993/ 1500 994/ 1500 995/ 1500 996/ 1500 997/ 1500 998/ 1500 999/ 1500 1000/ 1500 1001/ 1500 1002/ 1500 1003/ 1500 1004/ 1500 1005/ 1500 1006/ 1500 1007/ 1500 1008/ 1500 1009/ 1500 1010/ 1500 1011/ 1500 1012/ 1500 1013/ 1500 1014/ 1500 1015/ 1500 1016/ 1500 1017/ 1500 1018/ 1500 1019/ 1500 1020/ 1500 1021/ 1500 1022/ 1500 1023/ 1500 1024/ 1500 1025/ 1500 1026/ 1500 1027/ 1500 1028/ 1500 1029/ 1500 1030/ 1500 1031/ 1500 1032/ 1500 1033/ 1500 1034/ 1500 1035/ 1500 1036/ 1500 1037/ 1500 1038/ 1500 1039/ 1500 1040/ 1500 1041/ 1500 1042/ 1500 1043/ 1500 1044/ 1500 1045/ 1500 1046/ 1500 1047/ 1500 1048/ 1500 1049/ 1500 1050/ 1500 1051/ 1500 1052/ 1500 1053/ 1500 1054/ 1500 1055/ 1500 1056/ 1500 1057/ 1500 1058/ 1500 1059/ 1500 1060/ 1500 1061/ 1500 1062/ 1500 1063/ 1500 1064/ 1500 1065/ 1500 1066/ 1500 1067/ 1500 1068/ 1500 1069/ 1500 1070/ 1500 1071/ 1500 1072/ 1500 1073/ 1500 1074/ 1500 1075/ 1500 1076/ 1500 1077/ 1500 1078/ 1500 1079/ 1500 1080/ 1500 1081/ 1500 1082/ 1500 1083/ 1500 1084/ 1500 1085/ 1500 1086/ 1500 1087/ 1500 1088/ 1500 1089/ 1500 1090/ 1500 1091/ 1500 1092/ 1500 1093/ 1500 1094/ 1500 1095/ 1500 1096/ 1500 1097/ 1500 1098/ 1500 1099/ 1500 1100/ 1500 1101/ 1500 1102/ 1500 1103/ 1500 1104/ 1500 1105/ 1500 1106/ 1500 1107/ 1500 1108/ 1500 1109/ 1500 1110/ 1500 1111/ 1500 1112/ 1500 1113/ 1500 1114/ 1500 1115/ 1500 1116/ 1500 1117/ 1500 1118/ 1500 1119/ 1500 1120/ 1500 1121/ 1500 1122/ 1500 1123/ 1500 1124/ 1500 1125/ 1500 1126/ 1500 1127/ 1500 1128/ 1500 1129/ 1500 1130/ 1500 1131/ 1500 1132/ 1500 1133/ 1500 1134/ 1500 1135/ 1500 1136/ 1500 1137/ 1500 1138/ 1500 1139/ 1500 1140/ 1500 1141/ 1500 1142/ 1500 1143/ 1500 1144/ 1500 1145/ 1500 1146/ 1500 1147/ 1500 1148/ 1500 1149/ 1500 1150/ 1500 1151/ 1500 1152/ 1500 1153/ 1500 1154/ 1500 1155/ 1500 1156/ 1500 1157/ 1500 1158/ 1500 1159/ 1500 1160/ 1500 1161/ 1500 1162/ 1500 1163/ 1500 1164/ 1500 1165/ 1500 1166/ 1500 1167/ 1500 1168/ 1500 1169/ 1500 1170/ 1500 1171/ 1500 1172/ 1500 1173/ 1500 1174/ 1500 1175/ 1500 1176/ 1500 1177/ 1500 1178/ 1500 1179/ 1500 1180/ 1500 1181/ 1500 1182/ 1500 1183/ 1500 1184/ 1500 1185/ 1500 1186/ 1500 1187/ 1500 1188/ 1500 1189/ 1500 1190/ 1500 1191/ 1500 1192/ 1500 1193/ 1500 1194/ 1500 1195/ 1500 1196/ 1500 1197/ 1500 1198/ 1500 1199/ 1500 1200/ 1500 1201/ 1500 1202/ 1500 1203/ 1500 1204/ 1500 1205/ 1500 1206/ 1500 1207/ 1500 1208/ 1500 1209/ 1500 1210/ 1500 1211/ 1500 1212/ 1500 1213/ 1500 1214/ 1500 1215/ 1500 1216/ 1500 1217/ 1500 1218/ 1500 1219/ 1500 1220/ 1500 1221/ 1500 1222/ 1500 1223/ 1500 1224/ 1500 1225/ 1500 1226/ 1500 1227/ 1500 1228/ 1500 1229/ 1500 1230/ 1500 1231/ 1500 1232/ 1500 1233/ 1500 1234/ 1500 1235/ 1500 1236/ 1500 1237/ 1500 1238/ 1500 1239/ 1500 1240/ 1500 1241/ 1500 1242/ 1500 1243/ 1500 1244/ 1500 1245/ 1500 1246/ 1500 1247/ 1500 1248/ 1500 1249/ 1500 1250/ 1500 1251/ 1500 1252/ 1500 1253/ 1500 1254/ 1500 1255/ 1500 1256/ 1500 1257/ 1500 1258/ 1500 1259/ 1500 1260/ 1500 1261/ 1500 1262/ 1500 1263/ 1500 1264/ 1500 1265/ 1500 1266/ 1500 1267/ 1500 1268/ 1500 1269/ 1500 1270/ 1500 1271/ 1500 1272/ 1500 1273/ 1500 1274/ 1500 1275/ 1500 1276/ 1500 1277/ 1500 1278/ 1500 1279/ 1500 1280/ 1500 1281/ 1500 1282/ 1500 1283/ 1500 1284/ 1500 1285/ 1500 1286/ 1500 1287/ 1500 1288/ 1500 1289/ 1500 1290/ 1500 1291/ 1500 1292/ 1500 1293/ 1500 1294/ 1500 1295/ 1500 1296/ 1500 1297/ 1500 1298/ 1500 1299/ 1500 1300/ 1500 1301/ 1500 1302/ 1500 1303/ 1500 1304/ 1500 1305/ 1500 1306/ 1500 1307/ 1500 1308/ 1500 1309/ 1500 1310/ 1500 1311/ 1500 1312/ 1500 1313/ 1500 1314/ 1500 1315/ 1500 1316/ 1500 1317/ 1500 1318/ 1500 1319/ 1500 1320/ 1500 1321/ 1500 1322/ 1500 1323/ 1500 1324/ 1500 1325/ 1500 1326/ 1500 1327/ 1500 1328/ 1500 1329/ 1500 1330/ 1500 1331/ 1500 1332/ 1500 1333/ 1500 1334/ 1500 1335/ 1500 1336/ 1500 1337/ 1500 1338/ 1500 1339/ 1500 1340/ 1500 1341/ 1500 1342/ 1500 1343/ 1500 1344/ 1500 1345/ 1500 1346/ 1500 1347/ 1500 1348/ 1500 1349/ 1500 1350/ 1500 1351/ 1500 1352/ 1500 1353/ 1500 1354/ 1500 1355/ 1500 1356/ 1500 1357/ 1500 1358/ 1500 1359/ 1500 1360/ 1500 1361/ 1500 1362/ 1500 1363/ 1500 1364/ 1500 1365/ 1500 1366/ 1500 1367/ 1500 1368/ 1500 1369/ 1500 1370/ 1500 1371/ 1500 1372/ 1500 1373/ 1500 1374/ 1500 1375/ 1500 1376/ 1500 1377/ 1500 1378/ 1500 1379/ 1500 1380/ 1500 1381/ 1500 1382/ 1500 1383/ 1500 1384/ 1500 1385/ 1500 1386/ 1500 1387/ 1500 1388/ 1500 1389/ 1500 1390/ 1500 1391/ 1500 1392/ 1500 1393/ 1500 1394/ 1500 1395/ 1500 1396/ 1500 1397/ 1500 1398/ 1500 1399/ 1500 1400/ 1500 1401/ 1500 1402/ 1500 1403/ 1500 1404/ 1500 1405/ 1500 1406/ 1500 1407/ 1500 1408/ 1500 1409/ 1500 1410/ 1500 1411/ 1500 1412/ 1500 1413/ 1500 1414/ 1500 1415/ 1500 1416/ 1500 1417/ 1500 1418/ 1500 1419/ 1500 1420/ 1500 1421/ 1500 1422/ 1500 1423/ 1500 1424/ 1500 1425/ 1500 1426/ 1500 1427/ 1500 1428/ 1500 1429/ 1500 1430/ 1500 1431/ 1500 1432/ 1500 1433/ 1500 1434/ 1500 1435/ 1500 1436/ 1500 1437/ 1500 1438/ 1500 1439/ 1500 1440/ 1500 1441/ 1500 1442/ 1500 1443/ 1500 1444/ 1500 1445/ 1500 1446/ 1500 1447/ 1500 1448/ 1500 1449/ 1500 1450/ 1500 1451/ 1500 1452/ 1500 1453/ 1500 1454/ 1500 1455/ 1500 1456/ 1500 1457/ 1500 1458/ 1500 1459/ 1500 1460/ 1500 1461/ 1500 1462/ 1500 1463/ 1500 1464/ 1500 1465/ 1500 1466/ 1500 1467/ 1500 1468/ 1500 1469/ 1500 1470/ 1500 1471/ 1500 1472/ 1500 1473/ 1500 1474/ 1500 1475/ 1500 1476/ 1500 1477/ 1500 1478/ 1500 1479/ 1500 1480/ 1500 1481/ 1500 1482/ 1500 1483/ 1500 1484/ 1500 1485/ 1500 1486/ 1500 1487/ 1500 1488/ 1500 1489/ 1500 1490/ 1500 1491/ 1500 1492/ 1500 1493/ 1500 1494/ 1500 1495/ 1500 1496/ 1500 1497/ 1500 1498/ 1500 1499/ 1500 ```

``` 1500/ 1500 ```

```Generating estimation report...done. Discrete-time Neural State-Space Model with 15 outputs, 14 states, and 1 inputs x(t+1) = f(x(t),u(t)) y_1(t) = x(t) + e_1(t) y_2(t) = g(x(t),u(t)) + e_2(t) y(t) = [y_1(t); y_2(t)] f(.) network: Deep network with 1 fully connected, hidden layers Activation function: Sigmoid g(.) network: Deep network with 1 fully connected, hidden layers Activation function: Sigmoid Inputs: Torque Outputs: Angular Velocity(t-1), Angular Velocity(t-2), Angular Velocity(t-3), Angular Velocity(t-4), Angular Velocity(t-5), Angular Velocity(t-6), Angular Velocity(t-7), Torque(t-1), Torque(t-2), Torque(t-3), Torque(t-4), Torque(t-5), Torque(t-6), Torque(t-7), Angular Velocity States: x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14 Sample time: 0.01 seconds Status: Estimated using NLSSEST on time domain data "Data". Fit to estimation data: [65.19 58.06 63.49 63.41 62.98 65.14 64.52 31.35 12.49 8.173....% FPE: 2.292e-15, MSE: [6.141 6.324 6.006] More information in model's "Report" property. Model Properties ```
`toc`
```Elapsed time is 1516.352434 seconds. ```
```[yp,fit] = compare(XYData, nss); plot(XYData(:,end,[]),yp(:,end,[])) legend('measured',['model (',num2str(fit(end)),'%)'])```

This estimation takes a long time. For reference, the results are available in the data file `IdentifiedNeuralSSModel.mat`.

Just like the Hammerstein-Wiener model, a Neural State Space model is an Output Error model. This means that there is no process noise considered as part of the state transition. However, the identified model can be extended post-estimation with a process noise component.

```% Simulate the state trajectory [~,~,Xp] = sim(nss,XYData); % Fetch measured values of the states X_measured = XYData.OutputData(:,1:end-1); % Compute covariace of dX := Xp-X_measured dX = Xp-X_measured; XCov = cov(dX); NV = nss.NoiseVariance; Xsd = chol(XCov); NVc = chol(NV); NumParticles = 1000; x0 = X_measured(1,:)'```
```x0 = 14×1 -1.6494 -1.0918 -1.3285 -0.5497 -0.5271 -1.0293 -0.8808 -1.1322 -0.3761 -0.5067 ⋮ ```

Use the `generateMATLABFunction` command to generate MATLAB functions that implement the state-transition and the measurement functions. The functions are generated as M files. An accompanying MAT file is also generated and it must be present together with the M files for the interface to work properly in downstream applications.

Run the following command in a folder with write privileges to generate the Jacobian files `nssStateFcn.m, nssMeasurementFcn.m`, and the accompanying files.

`generateMATLABFunction(nss, "nssStateFcn", "nssMeasurementFcnOutput15");`

The generated functions - `nssStateFcn` and `nssMeasurementFcnOutput15 use persistent variables which are not allowed to be used inside MATLAB System blocks (the MATLAB System blocks are used to implement the predictor and corrector functions). Hence we modify these functions, post generation, to remove the use of the persistent variables. The modified functions are called nssStateFcn_2` and `nssMeasurementFcnOutput15_2 respectively.`

`type nssStateFcn_2`
```function x1 = nssStateFcn(x,u) %% auto-generated state function of neural state space system %# codegen persistent StateNetwork MATname = 'nssStateFcnData'; if isempty(StateNetwork) StateNetwork = coder.load(MATname); end out = [x;u]; % hidden layer #1 out = StateNetwork.fc1.Weights*out + StateNetwork.fc1.Bias; out = deep.internal.coder.sigmoid(out); % output layer dx = StateNetwork.output.Weights*out + StateNetwork.output.Bias; x1 = x + dx; ```
`type nssMeasurementFcnOutput15_2`
```function y = nssMeasurementFcnOutput15(x,u) %% auto-generated output function of neural state space system %# codegen persistent OutputNetwork MATname = 'nssMeasurementFcnOutput15Data'; if isempty(OutputNetwork) OutputNetwork = coder.load(MATname); end out = x; % hidden layer #1 out = OutputNetwork.fc1.Weights*out + OutputNetwork.fc1.Bias; out = deep.internal.coder.sigmoid(out); % output layer y = OutputNetwork.output.Weights*out + OutputNetwork.output.Bias; ```

Note that the measurement function file represents only the extra outputs (output numbers `nx+1:ny`). This is output signal number 15 in `XYData`. For use in EKF/UKF filters, we need a measurement function that describes all the outputs. Hence we wrap the `nssMeasurementFcnOutput15` function in another function called `nssMeasurementFcnAllOutputs` which appends the model states as its first 14 outputs.

`type nssMeasurementFcnAllOutputs`
```function y = nssMeasurementFcnAllOutputs(x,u) %% output function of neural state space system % all outputs %# codegen y15 = nssMeasurementFcnOutput15(x,u); y = [x; y15]; ```

The state-transition function (`nssStateFcn.m`) and the measurement function (`nssMeasurementFcnAllOutputs.m`) along with the process noise covariance (`XCov`) and measurement noise covariance (`NV`) can be used to set up online state estimators. This is implemented in `StateEstimation_Using_NeuralSS.slx`.

```Ts = nss.Ts; Ts2 = Ts*10; % correction time interval open_system('StateEstimation_Using_NeuralSS')```

The 10-step ahead prediction results are shown in the 3 scope plots (one for each type of filter), obtained for `XYData`.

`sim('StateEstimation_Using_NeuralSS')`
```ans = Simulink.SimulationOutput: tout: [99101x1 double] SimulationMetadata: [1x1 Simulink.SimulationMetadata] ErrorMessage: [0x0 char] ```

#### Nonlinear Grey Box Model

Finally, we explore the modeling approach that requires the maximum amount of physical knowledge regarding the dynamics. This results in a 5th order nonlinear differential equation. In state-space form these equations are:

`$\begin{array}{l}\mathbf{State}\text{\hspace{0.17em}}\mathbf{defintion}:\\ {\mathit{x}}_{1}\left(\mathit{t}\right)={\mathit{q}}_{\mathit{m}}\left(\mathit{t}\right)-{\mathit{q}}_{\mathit{g}}\left(\mathit{t}\right)\\ {\mathit{x}}_{2}\left(\mathit{t}\right)={\mathit{q}}_{\mathit{g}}\left(\mathit{t}\right)-{\mathit{q}}_{\mathit{a}}\left(\mathit{t}\right)\\ {\mathit{x}}_{3}\left(\mathit{t}\right)={\stackrel{˙}{\mathit{q}}}_{\mathit{m}}\left(\mathit{t}\right)\\ {\mathit{x}}_{4}\left(\mathit{t}\right)={\stackrel{˙}{\mathit{q}}}_{\mathit{g}}\left(\mathit{t}\right)\\ {\mathit{x}}_{5}\left(\mathit{t}\right)={\stackrel{˙}{\mathit{q}}}_{\mathit{a}}\left(\mathit{t}\right)\text{\hspace{0.17em}}\end{array}$`

`$\begin{array}{l}\mathbf{Equations}\text{\hspace{0.17em}}\mathbf{of}\text{\hspace{0.17em}}\mathbf{motion}:\\ {\stackrel{˙}{\mathit{x}}}_{1}\left(\mathit{t}\right)={\mathit{x}}_{3}\left(\mathit{t}\right)-{\mathit{x}}_{4}\left(\mathit{t}\right)\\ {\stackrel{˙}{\mathit{x}}}_{2}\left(\mathit{t}\right)={\mathit{x}}_{4}\left(\mathit{t}\right)-{\mathit{x}}_{5}\left(\mathit{t}\right)\\ {\stackrel{˙}{\mathit{x}}}_{3}\left(\mathit{t}\right)=\frac{1}{{\mathit{J}}_{\mathit{m}}}\left(-{\tau }_{\mathit{s}}\left(\mathit{t}\right)-{\mathit{d}}_{\mathit{g}}\left({\mathit{x}}_{3}\left(\mathit{t}\right)-{\mathit{x}}_{4}\left(\mathit{t}\right)\right)-{\tau }_{\mathit{f}}\left(\mathit{t}\right)+\mathit{u}\left(\mathit{t}\right)\right)\\ {\stackrel{˙}{\mathit{x}}}_{4}\left(\mathit{t}\right)=\frac{1}{{\mathit{J}}_{\mathit{g}}}\left({\tau }_{\mathit{s}}\left(\mathit{t}\right)+{\mathit{d}}_{\mathit{g}}\left({\mathit{x}}_{3}\left(\mathit{t}\right)-{\mathit{x}}_{4}\left(\mathit{t}\right)\right)-{\mathit{k}}_{\mathit{a}}{\mathit{x}}_{2}\left(\mathit{t}\right)-{\mathit{d}}_{\mathit{a}}\left({\mathit{x}}_{4}\left(\mathit{t}\right)-{\mathit{x}}_{5}\left(\mathit{t}\right)\right)\right)\\ {\stackrel{˙}{\mathit{x}}}_{5}\left(\mathit{t}\right)=\frac{1}{{\mathit{J}}_{\mathit{a}}}\left({\mathit{k}}_{\mathit{a}}{\mathit{x}}_{2}\left(\mathit{t}\right)+{\mathit{d}}_{\mathit{a}}\left({\mathit{x}}_{4}\left(\mathit{t}\right)-{\mathit{x}}_{5}\left(\mathit{t}\right)\right)\right)\\ \mathit{y}\left(\mathit{t}\right)={\mathit{x}}_{3}\left(\mathit{t}\right)\end{array}$`

where:

• ${\mathit{J}}_{\mathit{m}}.{\text{\hspace{0.17em}}\mathit{J}}_{\mathit{g}},{\text{\hspace{0.17em}}\mathit{J}}_{\mathit{a}}$ are the moments of inertia of the motor, the gear-box and the arm structure, respectively

• ${\mathit{d}}_{\mathit{g}},{\mathit{d}}_{\mathit{a}}$ are damping parameters

• ${\mathit{k}}_{\mathit{a}}$ is the stiffness of the arm structure.

• The gear-box friction torque ${\tau }_{\mathit{f}}\left(\mathit{t}\right)={\mathit{F}}_{\mathit{v}}{\mathit{x}}_{3}\left(\mathit{t}\right)+\left({\mathit{F}}_{\mathit{c}}+{\mathit{F}}_{\mathrm{cs}}\mathrm{sech}\left(\alpha {\mathit{x}}_{3}\left(\mathit{t}\right)\right)\right)\text{\hspace{0.17em}}\mathrm{tanh}\left(\beta {\mathit{x}}_{3}\left(\mathit{t}\right)\right)$, where ${\mathit{F}}_{\mathit{v}}$ and ${\mathit{F}}_{\mathit{c}}$ are the viscous and the Coulomb friction coefficients, ${\mathit{F}}_{\mathrm{cs}}$ and $\alpha$ are coefficients for reflecting the Stribeck effect, and $\beta$ is a parameter used to obtain a smooth transition from negative to positive velocities of ${\mathit{x}}_{3}\left(\mathit{t}\right)$.

• The torque of the spring ${\tau }_{\mathit{s}}\left(\mathit{t}\right)={\mathit{k}}_{\mathrm{g1}}{\mathit{x}}_{1}\left(\mathit{t}\right)+{\mathit{k}}_{\mathrm{g3}}{\mathit{x}}_{1}^{3}\left(\mathit{t}\right)$, where ${\mathit{k}}_{\mathrm{g1}},{\text{\hspace{0.17em}}\mathit{k}}_{\mathrm{g3}}$ are stiffness parameters of the gear-box spring.

The equations can be described by a Nonlinear Grey Box model (the `idnlgrey` object). The model structure for a nonlinear grey-box model is:

`$\begin{array}{l}\stackrel{˙}{\mathit{x}}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{t},\mathit{x},\mathit{u},\theta \right)\\ \mathit{y}\left(\mathit{t}\right)=\mathit{g}\left(\mathit{t},\mathit{x},\mathit{u},\theta \right)+\mathit{e}\left(\mathit{t}\right)\end{array}$`

where $\theta$ denotes the set of parameters used in the functions $\mathit{f}\left(.\right),\mathit{g}\left(.\right)$. The identification task is to estimate $\theta$ in order to minimize the difference between the simulated response of the model and its measured value. The procedure is similar to the one described above for the linear grey box case - the main task is to write an ODE file that returns the state derivatives $\stackrel{˙}{\mathit{x}}\left(\mathit{t}\right)$, and output $\mathit{y}\left(\mathit{t}\right)$ as a function of $\mathit{t},\mathit{u}\left(\mathit{t}\right),\mathit{x}\left(\mathit{t}\right)$ and the parameters. This exercise is the focus of the example "`Modeling an Industrial Robot Arm`". For current purposes, load the estimation results from this example.

```load RobotArmNonlinearGreyBoxModel nlgr nlgr.Name = 'Original Model (no process noise)'; nlgr```
```nlgr = Continuous-time nonlinear grey-box model defined by 'robotarm_c' (MEX-file): dx/dt = F(t, u(t), x(t), p1, ..., p13) y(t) = H(t, u(t), x(t), p1, ..., p13) + e(t) with 1 input(s), 5 state(s), 1 output(s), and 7 free parameter(s) (out of 13). Name: Original Model (no process noise) Status: Estimated using Solver: ode45; Search: lsqnonlin on time domain data "estimation data". Fit to estimation data: 85.52% FPE: 9.184, MSE: 9.092 Model Properties ```
`compare(vData, nlgr)`

Similar to the linear case, you can embed the identified nonlinear grey box model into a more general stochastic framework by adding process noise to the equations for $\stackrel{˙}{\mathit{x}}\left(\mathit{t}\right)$, that is, pose $\stackrel{˙}{\mathit{x}}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{t},\mathit{x},\mathit{u},\theta \right)+\mathrm{Gw}\left(\mathit{t}\right)$, where $\mathrm{G}$ is the process noise gain matrix. However, there is no direct way to compute the value of $\mathrm{G}$ under the nonlinear grey-box estimation framework. This is because the Nonlinear Grey Box methodology only handles Output-Error models, that is, it does not consider the influence of process noise.

An option is to pose a simple extension $\stackrel{˙}{\mathit{x}}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{t},\mathit{x},\mathit{u},\theta \right)+\mathrm{Ke}\left(\mathit{t}\right)$, where $\mathit{e}\left(\mathit{t}\right)$ is the output prediction error. This is similar to the innovations form of the linear state-space model, where $\stackrel{˙}{\mathit{x}}\left(\mathit{t}\right)=\mathrm{Ax}\left(\mathit{t}\right)+\mathrm{Bu}\left(\mathit{t}\right)+\mathrm{Ke}\left(\mathit{t}\right)$

For prediction of model output, the nonlinear predictor therefore takes the form:

`$\begin{array}{l}{\stackrel{˙}{\mathit{x}}}_{\mathit{p}}\left(\mathit{t}\right)=\mathit{f}\left(\mathit{t},\mathit{x},\mathit{u},\theta \right)+\mathit{K}\left({\mathit{y}}_{\mathrm{measured}}-{\mathit{y}}_{\mathit{p}}\left(\mathit{t}\right)\right)\\ {\mathit{y}}_{\mathit{p}}\left(\mathit{t}\right)=\mathit{g}\left(\mathit{t},\mathit{x},\mathit{u},\theta \right)\end{array}$`

In order to make this extension, add $\mathit{K}$as a free parameter (`5x1` vector) to the model `nlgr`. Also, the ODE function needs to compute ${\mathit{y}}_{\mathrm{measured}}-{\mathit{y}}_{\mathit{p}}\left(\mathit{t}\right)$. Replace the original ODE function employed by `nlgr` (`robot_c`) with a new one called `robotARMNonlinearODE`. For this, pass a gridded interpolator as a function argument.

`type robotARMNonlinearODE`
```function [dx, y] = robotARMNonlinearODE(t, x, u, Fv, Fc, Fcs, alpha, beta, J, am, ... ag, kg1, kg3, dg, ka, da, K, Interpolator) %robotARMNonlinearODE A physically parameterized robot arm. % t: time instant (seconds) % x: state value at time t (5 x 1 vector) % u: input at time t (scalar) % Fv, Fc, Fcs, alpha, beta, J, am, ag, kg1, kg3, dg, ka, da: parameters related to the model % dynamics % K: process noise gain matrix (free parameter) % Interpolator: griddedInterpolant object to fetch measured output at a time t % Copyright 2005-2022 The MathWorks, Inc. % Output equation. y = x(3); % Rotational velocity of the motor. % Intermediate quantities (gear box). tauf = Fv*x(3)+(Fc+Fcs/cosh(alpha*x(3)))*tanh(beta*x(3)); % Gear friction torque. taus = kg1*x(1)+kg3*x(1)^3; % Spring torque. % State equations. dx = [x(3)-x(4); ... % Rotational velocity difference between the motor and the gear-box. x(4)-x(5); ... % Rotational velocity difference between the gear-box and the arm. 1/(J*am)*(-taus-dg*(x(3)-x(4))-tauf+u); ... % Rotational velocity of the motor. 1/(J*ag)*(taus+dg*(x(3)-x(4))-ka*x(2)-da*(x(4)-x(5))); ... % Rotational velocity after the gear-box. 1/(J*(1-am-ag))*(ka*x(2)+da*(x(4)-x(5))) ... % Rotational velocity of the robot arm. ]; ym = Interpolator{1}(t); dx = dx + K*(ym-y); ```

The last line of code adds a correction term to the state derivatives based on the output prediction error. In order to compute the measured output value at a given time instant, we can simply interpolate the values in `eData.OutputData`. To enable this, we pass a gridded interpolant as an extra input argument to the ODE function (see the last input argument "`Interpolator`"). When setting up the nonlinear grey-box model (`idnlgrey` object), you can specify any extra arguments to the ODE function using the `FileArgument` property.

```% Change the ODE function of nlgr from "motor_c" to "robotARMNonlinearODE". robotARMNonlinearODE is % mostly identical to motor_c except for the addition of the correction to dx induced by the process % noise. PredictorModel = nlgr; PredictorModel.Name = 'Predictor form'; PredictorModel.FileName = 'robotARMNonlinearODE'; % Create a new parameter to represent the process noise gain K % Set 0 as initial guess for K entries NoiseGain = struct('Name','K','Unit','','Value',zeros(5,1),... 'Minimum',-Inf(5,1),'Maximum',Inf(5,1),'Fixed',false(5,1)); % Add this parameter to the model PredictorModel.Parameters(end+1) = NoiseGain; % Create a linear interpolator t = eData.SamplingInstants; y = eData.OutputData; % measured output Interpolator = griddedInterpolant(t,y); Interpolator.ExtrapolationMethod = 'nearest'; PredictorModel.FileArgument = {Interpolator};```

Perform parameter estimation, Because of the process noise correction, the estimation is implicitly minimizing the 1-step ahead prediction error (rather than the default behavior which is to minimize the simulation error)

```opt = nlgreyestOptions('Display', 'on','SearchMethod','lm'); opt.SearchOptions.MaxIterations = 20; opt.EstimateCovariance = false; % to speed up estimation % Keep the parameter estimates close to their initial guesses. This is because we already have a % good guess of their values from previous estimation as carried out in the example "Modeling an % Industrial Robot Arm" opt.Regularization.Lambda = .1; opt.Regularization.Nominal = 'model'; % Estimate the model parameters PredictorModel = nlgreyest(eData, PredictorModel, opt)```
```PredictorModel = Continuous-time nonlinear grey-box model defined by 'robotARMNonlinearODE' (MATLAB file): dx/dt = F(t, u(t), x(t), p1, ..., p14, FileArgument) y(t) = H(t, u(t), x(t), p1, ..., p14, FileArgument) + e(t) with 1 input(s), 5 state(s), 1 output(s), and 12 free parameter(s) (out of 18). Name: Predictor form Status: Estimated using Solver: ode45; Search: lm on time domain data "estimation data". Fit to estimation data: 90.87% FPE: 3.674, MSE: 3.618 Model Properties ```
`compare(eData, nlgr, PredictorModel)`

```% Fetch the computed process noise gain matrix K = PredictorModel.Parameters(end).Value```
```K = 5×1 0.3661 -0.4356 44.7611 -7.1629 -0.3057 ```
```% Fetch the output noise variance NV = PredictorModel.NoiseVariance```
```NV = 0.0182 ```

Validate the model.

```% Check 1-step ahead prediction quality on validation dataset tv = vData.SamplingInstants; yv = vData.OutputData; % measured output Interpolator2 = griddedInterpolant(tv,yv); Interpolator2.ExtrapolationMethod = 'nearest'; PredictorModel.FileArgument = {Interpolator2}; compare(vData, nlgr, PredictorModel)```

As seen from both the comparisons to the training data (`eData`) as well as the validation data (`vData`), the correction applied to the state-transition using the noise gain K has a beneficial effect on 1-step ahead predictions. Note that the original model (`nlgr`) has no process noise component. Hence for this model, 1-step ahead prediction is same as pure simulation.

You now have all the information available to implement nonlinear state estimators. Use the functions $\mathit{f}\left(.\right),\text{\hspace{0.17em}}\mathit{g}\left(.\right)$, as implemented in `robotARMNonlinearODE` to specify the state transition and measurement functions for the filters. The process noise is zero-mean Gaussian with a variance equal to $\mathit{E}\left\{{\mathrm{ww}}^{\prime }\right\}=\mathit{K}\Sigma {\mathit{K}}^{\mathit{T}}$, where $\Sigma$`=NV` is the measurement covariance.

Compute 10-step ahead prediction of the `vData` output using EKF, UKF, and PF state estimators based on `PredictorModel`.
```% Gather block data Ts = eData.Ts; % prediction interval Ts2 = 10*Ts; % correction interval % Parameter variables for use in the filter blocks [Fv, Fc, Fcs, alpha, beta, J, am, ag, kg1, kg3, dg, ka, da, K] = deal(PredictorModel.Parameters.Value); NumParticles = 1000; % number of particles used by the Particle filter open_system('StateEstimation_Using_NonlinearGreyBox')```
The Simulink model is set up to use a fixed step discrete solver with a fixed-step size of `Ts/100 = 50 ms`. The state transition equation used by the filters is implemented in the `nlgreyPredictFcn` function. This model is setup in a manner very similar to the one for Nonlinear ARX model based state prediction (see `StateEstimation_Using_NonlinearARX.slx`)
`type nlgreyPredictFcn`
```function xnext = nlgreyPredictFcn(x,u,Fv, Fc, Fcs, alpha, beta, J, am, ... ag, kg1, kg3, dg, ka, da, Ts) % Predict the state update for the nonlinear grey box model of robot arm. % State transition function for the UKF/EKF block. % Ts: data sample time. The model has been configured to use a fixed step size of Ts/100. % Copyright 2022 The MathWorks, Inc. % Euler integration of continuous-time dynamics x'=f(x) with sample time dt xnext = x; for i = 1:size(x,2) xnext(:,i) = x(:,i) + continuousTimeDerivative(x(:,i),u,... Fv,Fc,Fcs,alpha,beta,J,am,ag,kg1,kg3,dg,ka,da)*Ts/100; end end %--------------------------------------------------------------------------------------------------- function dx = continuousTimeDerivative(x,u,Fv, Fc, Fcs, alpha, beta, J, am, ... ag, kg1, kg3, dg, ka, da) % Robot arm state derivatives. tauf = Fv*x(3)+(Fc+Fcs/cosh(alpha*x(3)))*tanh(beta*x(3)); % Gear friction torque. taus = kg1*x(1)+kg3*x(1)^3; % Spring torque. % State equations. dx = [x(3)-x(4); ... % Rotational velocity difference between the motor and the gear-box. x(4)-x(5); ... % Rotational velocity difference between the gear-box and the arm. 1/(J*am)*(-taus-dg*(x(3)-x(4))-tauf+u); ... % Rotational velocity of the motor. 1/(J*ag)*(taus+dg*(x(3)-x(4))-ka*x(2)-da*(x(4)-x(5))); ... % Rotational velocity after the gear-box. 1/(J*(1-am-ag))*(ka*x(2)+da*(x(4)-x(5))) ... % Rotational velocity of the robot arm. ]; end ```
`sim('StateEstimation_Using_NonlinearGreyBox');`