Nonlinear ARX Model of SI Engine Torque Dynamics

This example describes modeling the nonlinear torque dynamics of a spark-ignition (SI) engine as a nonlinear ARX model. The identified model can be used for hardware-in-the-loop (HIL) testing, powertrain control, diagnostic, and training algorithm design. For example, you can use the model for aftertreatment control and diagnostics algorithm development. For more information on nonlinear ARX models, see Nonlinear ARX Models.

You use measurements of the system inputs and outputs to identify the model. This data can come from measurements on a real engine or a high-fidelity model such as one created using the Powertrain Blockset™ SI Reference application.

Data Preparation

Load and display the engine data timetable.

```load SIEngineData IOData stackedplot(IOData) title('SI Engine Signals')```

The timetable contains over one hundred thousands observations of five variables measured at 10 Hz.

• Throttle position (degrees)

• Wastegate valve area (aperture percentage)

• Engine speed (RPM)

• Spark timing (degrees)

• Engine torque (N m)

Split the data into estimation (first 60000 data points) and validation (remaining data points) portions.

```eData = IOData(1:6e4,:); % portion used for estimation vData = IOData(6e4+1:end,:); % portion used for validation```

Downsample the training data by a factor of 10. This helps speed up the model training process and also limits the focus of the fit to a lower frequency region.

```% Downsample datasets 10 times eDataD = idresamp(eData,[1 10]); vDataD = idresamp(vData,[1 10]); eDataD.Properties.TimeStep```
```ans = duration 1 sec ```

The model training objective is to generate a dynamic system that has the engine torque as output and the other four variables as inputs.

Nonlinear ARX Model Identification

About nonlinear ARX Models

A nonlinear ARX model implements a two-stage transformation of the observed data to predict the future values. In the first step, it transforms the input/output signals into a finite-dimensional regressors, which are features based on time-delayed values of the signals. These regressors are mapped to the predicted output using a nonlinear but static function. You have many choices for the regressor formulas and the nonlinear functions.

You can use linear, polynomial, or trigonometric transformations of the time-delayed variables. You can also use your physical insight to define your own custom regressors. The functions that map the regressors to the outputs usually employ a parallel connection combination of a linear function, a nonlinear function and an offset (bias) term. You can chose the nonlinear functions to be sigmoid or wavelet based networks, a Gaussian process (GP) function, support vector regression (SVM) functions, regression tree ensembles (bagging or boosting forms), or other basis function expansions expressed in terms of custom unit functions.

In this example, you identify nonlinear ARX models based on different nonlinear functions. Begin by designating the input and output signals from the list of variables in the `eData` timetable.

```Inputs = ["ThrottlePosition","WastegateValve","EngineSpeed","SparkTiming"]; Output = "EngineTorque";```

Sigmoid Network Model

A sigmoid network model is a nonlinear ARX model that employs a parallel combination of a linear function, an offset term (bias), and a nonlinear function represented by a sum of dilated and translated sigmoid functions. The sigmoid network is represented by the `idSigmoidNetwork` object.

The parameters of the nonlinear ARX model are composed of parameters of the linear function, the nonlinear function and the offset. The modeling approach used in this example is incremental. You first estimate a linear, four-input, one-output model for the torque dynamics. Then, you extend the linear block by adding a single-hidden-layer sigmoid network to it in a parallel configuration to create a nonlinear ARX model.

Initial Linear Model

Create an object to specify options for estimating state-space models, and specify the option to minimize the simulation error (instead of the prediction error) during estimation.

`ssOpt = ssestOptions(Focus="simulation");`

You can uncomment and run the following command to launch an interactive plot for selecting both the optimal order and the corresponding `N4Horizon` parameter value.

```% linsys = ssest(eDataD, 1:10, ... % 'Feedthrough',true(1,4), ... % 'Ts',1, ... % 'DisturbanceModel','none', ... % 'InputName',Inputs, ... % 'OutputName',Output, ... % ssOpt);```

When you click Apply, an identified state-space model with your chosen order (here, 5) is created. Use the following commands to directly create this model without using the plot.

```% Specify the horizon and the order ssOpt.N4Horizon = [15 0 25]; Order = 5; % Estimate the linear system linsys = ssest(eDataD, Order, ... 'Feedthrough',true(1,4), ... 'Ts',1,... 'DisturbanceModel','none',... 'InputName',Inputs,... 'OutputName',Output,... ssOpt)```
```linsys = Discrete-time identified state-space model: x(t+Ts) = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x4 x5 x1 0.685 -0.2061 0.2727 0.1345 0.1283 x2 0.3301 -0.01208 0.6969 0.2742 0.3221 x3 0.0977 0.3312 0.4462 -0.5481 0.07495 x4 -0.1477 -0.1594 0.9024 -0.2664 0.0361 x5 0.0135 0.606 -0.4236 -0.2822 0.6817 B = ThrottlePosi WastegateVal EngineSpeed SparkTiming x1 0.0004832 -0.01389 0.0007426 0.001756 x2 0.004744 -0.08037 0.004403 0.008711 x3 0.002218 -0.04952 0.002774 0.002881 x4 0.002496 -0.04896 0.002671 0.003919 x5 -0.002972 0.0553 -0.003039 -0.005904 C = x1 x2 x3 x4 x5 EngineTorque -557.9 175.6 -85.25 10.78 49.91 D = ThrottlePosi WastegateVal EngineSpeed SparkTiming EngineTorque 0.1519 0.1686 0.001073 0.05432 K = EngineTorque x1 0 x2 0 x3 0 x4 0 x5 0 Sample time: 1 seconds Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: yes Disturbance component: none Number of free coefficients: 54 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "eDataD". Fit to estimation data: 46.45% FPE: 173, MSE: 171.3 Model Properties ```

Compare the simulated response of linsys against the measured engine torque used for training.

```clf compare(eDataD, linsys)```

The linear model is able to emulate the measured output, but with about a normalized root mean squared error (NRMSE) of 50%. However, as shown next, this linear model provides an excellent starting point for building a nonlinear ARX model.

Creating the Nonlinear Model

In the Sigmoid network based nonlinear ARX model, you have the option of disabling the linear function, the offset term, or both. You can also pick how many sigmoid units the nonlinear function employs. In this example, you retain the linear component and the offset term but change the number of sigmoid units to 11 (default is 10). The number of coefficients in the linear term and their initial values are set by using in the linear model `linsys` in the call to the model training command `nlarx`.

Create the nonlinear sigmoid-based function.

```NumberOfUnits = 11; % determined by trial and error NonlinearFcn = idSigmoidNetwork(NumberOfUnits)```
```NonlinearFcn = Sigmoid Network Nonlinear Function: Sigmoid network with 11 units Linear Function: uninitialized Output Offset: uninitialized Inputs: {1×0 cell} Outputs: {1×0 cell} NonlinearFcn: 'Sigmoid units and their parameters' LinearFcn: 'Linear function parameters' Offset: 'Offset parameters' ```

For model training, use the Levenberg-Marquardt search method and set the focus of training to "simulation" (minimize the error between the measured output and the model's simulated response). Also, normalize the model features and output using the "zscore" method.

```% Configure training options nlarxopt = nlarxOptions(Focus="simulation"); nlarxopt.SearchMethod = 'lm'; nlarxopt.SearchOptions.MaxIterations = 100; % Normalize regressors and output during training nlarxopt.NormalizationOptions.NormalizationMethod = 'zscore'; nlarxopt.Display = 'on';```

Use the data (`eDataD`), the linear model (`linsys`), and the nonlinear function (`NonlinearFcn`) in the `nlarx` command to train the sigmoid-network-based nonlinear ARX model.

```tic SigmoidNet = nlarx(eDataD,linsys,NonlinearFcn,nlarxopt)```
```SigmoidNet = Nonlinear ARX model with 1 output and 4 inputs Inputs: ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming Outputs: EngineTorque Regressors: Linear regressors in variables EngineTorque, ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming List of all regressors Output function: Sigmoid network with 11 units Sample time: 1 seconds Status: Termination condition: No improvement along the search direction with line search.. Number of iterations: 6, Number of function evaluations: 65 Estimated using NLARX on time domain data "eDataD". Fit to estimation data: 91.44% (simulation focus) FPE: 0.2113, MSE: 4.376 More information in model's "Report" property. Model Properties ```
`toc`
```Elapsed time is 149.138082 seconds. ```

The model `SigmoidNet` creates a set of linear regressors derived from the linear model linsys. It also uses the linear model coefficients to initialize its linear function coefficients.

`getreg(SigmoidNet)`
```ans = 29×1 cell {'EngineTorque(t-1)' } {'EngineTorque(t-2)' } {'EngineTorque(t-3)' } {'EngineTorque(t-4)' } {'EngineTorque(t-5)' } {'ThrottlePosition(t)' } {'ThrottlePosition(t-1)'} {'ThrottlePosition(t-2)'} {'ThrottlePosition(t-3)'} {'ThrottlePosition(t-4)'} {'ThrottlePosition(t-5)'} {'WastegateValve(t)' } {'WastegateValve(t-1)' } {'WastegateValve(t-2)' } {'WastegateValve(t-3)' } {'WastegateValve(t-4)' } {'WastegateValve(t-5)' } {'EngineSpeed(t)' } {'EngineSpeed(t-1)' } {'EngineSpeed(t-2)' } {'EngineSpeed(t-3)' } {'EngineSpeed(t-4)' } {'EngineSpeed(t-5)' } {'SparkTiming(t)' } {'SparkTiming(t-1)' } {'SparkTiming(t-2)' } {'SparkTiming(t-3)' } {'SparkTiming(t-4)' } {'SparkTiming(t-5)' } ```

Do a quick validation of the model quality by comparing its simulated response against the measured responses in the estimation and the validation datasets. Note that you check only against the decimated data here.

```subplot(211) compare(eDataD,SigmoidNet) % compare against estimation data title('Sigmoid Network Model: Comparison to Estimation Data') subplot(212) compare(vDataD,SigmoidNet) % compare against validation data title('Comparison to Validation Data')```

Gaussian Process (GP) Model

You can similarly create GP representation of the nonlinear dynamics.

Use an `idGaussianProcess` object to represent the nonlinear function. Create a GP representation based on the ARD Matern32 kernel. For this model, you do not use the linear model to initialize the model structure. Instead, create the regressors directly using the `linearRegressor` command.

```Kernel = 'ARDMatern32'; GP = idGaussianProcess(Kernel); MaxLag = 4; % maximum number of delays in each regressor OutputLags = 1:MaxLag; InputLags = 0:MaxLag; Reg = linearRegressor([Output, Inputs],{OutputLags,InputLags,InputLags,InputLags,InputLags})```
```Reg = Linear regressors in variables EngineTorque, ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming Variables: {'EngineTorque' 'ThrottlePosition' 'WastegateValve' 'EngineSpeed' 'SparkTiming'} Lags: {[1 2 3 4] [0 1 2 3 4] [0 1 2 3 4] [0 1 2 3 4] [0 1 2 3 4]} UseAbsolute: [0 0 0 0 0] TimeVariable: 't' Regressors described by this set ```

Use the estimation data (`eDataD`), the regressor set (`Reg`), and the nonlinear function (`GP`) in the `nlarx` command to train the GP-based nonlinear ARX model.

```nlarxopt.SearchOptions.MaxIterations = 20; tic GPModel = nlarx(eDataD,Reg,GP,nlarxopt)```
```|=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 0 | -4.31381e+03 | 6.390e+02 | 1.032e+00 | CONV | 2.065e+00 | 5.099e+00 | YES | | 1 | -4.31381e+03 | 6.390e+02 | 5.099e+00 | EXACT | -2.375e+00 | 5.099e+00 | NO | | 2 | -5.32039e+03 | 8.618e+02 | 2.550e+00 | EXACT | 2.711e-01 | 2.550e+00 | YES | | 3 | -5.34697e+03 | 8.088e+02 | 2.942e-02 | CONV | 1.947e+00 | 2.550e+00 | YES | | 4 | -5.48306e+03 | 7.906e+02 | 5.351e-01 | CONV | 5.787e-01 | 2.550e+00 | YES | | 5 | -5.58226e+03 | 2.694e+02 | 2.963e-01 | CONV | 7.755e-01 | 2.550e+00 | YES | | 6 | -6.01284e+03 | 1.156e+02 | 1.758e+00 | CONV | 1.627e+00 | 2.550e+00 | YES | | 7 | -6.57309e+03 | 8.562e+01 | 2.550e+00 | BNDRY | 1.412e+00 | 2.550e+00 | YES | | 8 | -7.52416e+03 | 4.205e+02 | 5.099e+00 | BNDRY | 1.280e+00 | 5.099e+00 | YES | | 9 | -7.52416e+03 | 4.205e+02 | 1.020e+01 | EXACT | -2.768e-01 | 1.020e+01 | NO | | 10 | -7.94111e+03 | 1.253e+02 | 5.099e+00 | BNDRY | 7.350e-01 | 5.099e+00 | YES | | 11 | -7.95215e+03 | 5.025e+02 | 1.100e+00 | CONV | 2.483e-01 | 5.099e+00 | YES | | 12 | -7.98119e+03 | 2.517e+02 | 7.176e-02 | CONV | 1.508e+00 | 5.099e+00 | YES | | 13 | -7.99465e+03 | 5.390e+01 | 8.582e-02 | CONV | 1.111e+00 | 5.099e+00 | YES | | 14 | -8.06705e+03 | 6.143e+01 | 3.915e+00 | CONV | 9.341e-01 | 5.099e+00 | YES | | 15 | -8.06705e+03 | 6.143e+01 | 5.099e+00 | EXACT | -1.172e-01 | 5.099e+00 | NO | | 16 | -8.10385e+03 | 2.499e+02 | 2.550e+00 | BNDRY | 3.154e-01 | 2.550e+00 | YES | | 17 | -8.15516e+03 | 7.925e+01 | 3.375e-01 | CONV | 1.513e+00 | 2.550e+00 | YES | | 18 | -8.15516e+03 | 7.925e+01 | 2.550e+00 | EXACT | -1.238e-01 | 2.550e+00 | NO | | 19 | -8.18832e+03 | 4.825e+01 | 7.055e-01 | CONV | 9.162e-01 | 1.275e+00 | YES | |=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 20 | -8.20116e+03 | 2.206e+01 | 7.438e-01 | CONV | 1.322e+00 | 1.275e+00 | YES | | 21 | -8.21145e+03 | 3.409e+01 | 1.037e+00 | CONV | 1.281e+00 | 1.275e+00 | YES | | 22 | -8.21145e+03 | 3.409e+01 | 2.550e+00 | BNDRY | -3.629e+00 | 2.550e+00 | NO | | 23 | -8.21145e+03 | 3.409e+01 | 1.126e+00 | CONV | -5.050e-01 | 1.275e+00 | NO | | 24 | -8.22018e+03 | 3.436e+01 | 5.504e-01 | CONV | 8.059e-01 | 6.374e-01 | YES | | 25 | -8.22136e+03 | 2.038e+01 | 3.781e-02 | CONV | 1.832e+00 | 1.275e+00 | YES | | 26 | -8.22597e+03 | 5.087e+01 | 3.347e-01 | CONV | 1.068e+00 | 1.275e+00 | YES | | 27 | -8.22597e+03 | 5.087e+01 | 1.275e+00 | EXACT | -5.391e-01 | 1.275e+00 | NO | | 28 | -8.22626e+03 | 3.608e+01 | 5.767e-03 | CONV | 1.734e+00 | 6.374e-01 | YES | | 29 | -8.22662e+03 | 7.018e+00 | 1.710e-02 | CONV | 9.861e-01 | 6.374e-01 | YES | | 30 | -8.23051e+03 | 2.485e+01 | 6.374e-01 | EXACT | 1.768e-01 | 6.374e-01 | YES | | 31 | -8.23420e+03 | 9.600e+00 | 6.374e-01 | BNDRY | 8.754e-01 | 6.374e-01 | YES | | 32 | -8.23669e+03 | 1.879e+01 | 7.030e-01 | CONV | 1.156e+00 | 1.275e+00 | YES | | 33 | -8.23669e+03 | 1.879e+01 | 1.275e+00 | EXACT | -1.072e+00 | 1.275e+00 | NO | | 34 | -8.23828e+03 | 1.612e+01 | 2.212e-01 | CONV | 9.035e-01 | 6.374e-01 | YES | | 35 | -8.23846e+03 | 1.247e+01 | 1.462e-02 | CONV | 1.865e+00 | 6.374e-01 | YES | | 36 | -8.24005e+03 | 8.546e+00 | 3.812e-01 | CONV | 1.064e+00 | 6.374e-01 | YES | | 37 | -8.24005e+03 | 8.546e+00 | 6.374e-01 | BNDRY | -9.651e-01 | 6.374e-01 | NO | | 38 | -8.24105e+03 | 5.501e+00 | 3.187e-01 | BNDRY | 1.035e+00 | 3.187e-01 | YES | | 39 | -8.24105e+03 | 5.501e+00 | 6.374e-01 | EXACT | -1.583e-01 | 6.374e-01 | NO | |=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 40 | -8.24211e+03 | 4.683e+00 | 3.187e-01 | BNDRY | 1.083e+00 | 3.187e-01 | YES | | 41 | -8.24332e+03 | 4.315e+00 | 6.374e-01 | BNDRY | 1.146e+00 | 6.374e-01 | YES | | 42 | -8.24399e+03 | 8.473e+00 | 6.155e-01 | CONV | 1.278e+00 | 1.275e+00 | YES | | 43 | -8.24418e+03 | 3.204e+00 | 2.804e-02 | CONV | 1.659e+00 | 1.275e+00 | YES | | 44 | -8.24438e+03 | 1.932e+01 | 3.526e-01 | CONV | 3.166e-01 | 1.275e+00 | YES | | 45 | -8.24449e+03 | 6.460e+00 | 6.603e-03 | CONV | 1.349e+00 | 1.275e+00 | YES | | 46 | -8.24482e+03 | 2.239e+00 | 8.697e-02 | CONV | 1.248e+00 | 1.275e+00 | YES | | 47 | -8.24528e+03 | 2.281e+00 | 6.766e-01 | CONV | 9.778e-01 | 1.275e+00 | YES | | 48 | -8.24530e+03 | 2.263e+00 | 2.098e-03 | CONV | 1.988e+00 | 1.275e+00 | YES | | 49 | -8.24672e+03 | 1.415e+01 | 1.275e+00 | BNDRY | 7.882e-01 | 1.275e+00 | YES | | 50 | -8.24694e+03 | 4.008e+00 | 2.013e-02 | CONV | 1.620e+00 | 2.550e+00 | YES | | 51 | -8.24779e+03 | 2.804e+00 | 6.365e-01 | CONV | 9.380e-01 | 2.550e+00 | YES | | 52 | -8.24791e+03 | 1.727e+00 | 1.250e-01 | CONV | 1.233e+00 | 2.550e+00 | YES | | 53 | -8.24791e+03 | 1.727e+00 | 1.104e+00 | CONV | -8.897e-01 | 2.550e+00 | NO | | 54 | -8.24813e+03 | 1.891e+00 | 4.673e-01 | CONV | 1.240e+00 | 1.275e+00 | YES | | 55 | -8.24821e+03 | 1.444e+00 | 6.598e-02 | CONV | 1.400e+00 | 1.275e+00 | YES | | 56 | -8.24827e+03 | 2.165e+00 | 1.506e-01 | CONV | 6.184e-01 | 1.275e+00 | YES | | 57 | -8.24830e+03 | 1.307e+00 | 2.383e-02 | CONV | 1.401e+00 | 1.275e+00 | YES | | 58 | -8.24830e+03 | 1.307e+00 | 1.275e+00 | BNDRY | -1.176e+00 | 1.275e+00 | NO | | 59 | -8.24872e+03 | 1.744e+00 | 5.747e-01 | CONV | 1.052e+00 | 6.374e-01 | YES | |=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 60 | -8.24872e+03 | 1.744e+00 | 2.782e-01 | CONV | -7.533e-01 | 1.275e+00 | NO | | 61 | -8.24881e+03 | 2.095e+00 | 8.867e-02 | CONV | 1.184e+00 | 6.374e-01 | YES | | 62 | -8.24898e+03 | 3.838e+00 | 2.840e-01 | CONV | 1.297e+00 | 6.374e-01 | YES | | 63 | -8.24900e+03 | 1.089e+00 | 8.141e-03 | CONV | 1.529e+00 | 6.374e-01 | YES | | 64 | -8.24900e+03 | 1.089e+00 | 6.374e-01 | EXACT | -1.407e+00 | 6.374e-01 | NO | | 65 | -8.24915e+03 | 1.202e+00 | 3.187e-01 | BNDRY | 1.168e+00 | 3.187e-01 | YES | | 66 | -8.24932e+03 | 1.885e+00 | 6.374e-01 | BNDRY | 1.156e+00 | 6.374e-01 | YES | | 67 | -8.24936e+03 | 4.521e-01 | 7.310e-02 | CONV | 1.016e+00 | 1.275e+00 | YES | | 68 | -8.24936e+03 | 4.521e-01 | 1.275e+00 | EXACT | -5.207e-01 | 1.275e+00 | NO | | 69 | -8.24947e+03 | 1.304e+00 | 6.374e-01 | BNDRY | 1.118e+00 | 6.374e-01 | YES | | 70 | -8.24950e+03 | 1.101e+00 | 9.168e-02 | CONV | 8.865e-01 | 1.275e+00 | YES | | 71 | -8.24950e+03 | 4.084e-01 | 2.281e-03 | CONV | 1.818e+00 | 1.275e+00 | YES | | 72 | -8.24966e+03 | 1.578e+00 | 1.148e+00 | CONV | 1.264e+00 | 1.275e+00 | YES | | 73 | -8.24966e+03 | 5.001e-01 | 2.343e-03 | CONV | 1.613e+00 | 2.550e+00 | YES | | 74 | -8.24966e+03 | 5.001e-01 | 2.550e+00 | EXACT | -3.572e+00 | 2.550e+00 | NO | | 75 | -8.24969e+03 | 4.874e-01 | 1.959e-01 | CONV | 1.030e+00 | 1.275e+00 | YES | | 76 | -8.24969e+03 | 4.874e-01 | 1.275e+00 | EXACT | -1.435e-01 | 1.275e+00 | NO | | 77 | -8.24970e+03 | 3.437e+00 | 1.852e-01 | CONV | 2.884e-01 | 6.374e-01 | YES | | 78 | -8.24970e+03 | 1.016e+00 | 1.305e-03 | CONV | 1.172e+00 | 6.374e-01 | YES | | 79 | -8.24971e+03 | 3.651e-01 | 1.948e-02 | CONV | 1.042e+00 | 6.374e-01 | YES | |=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 80 | -8.24971e+03 | 3.651e-01 | 6.374e-01 | EXACT | -9.870e-02 | 6.374e-01 | NO | | 81 | -8.24973e+03 | 7.005e-01 | 1.584e-01 | CONV | 9.709e-01 | 3.187e-01 | YES | | 82 | -8.24973e+03 | 5.043e-01 | 1.043e-03 | CONV | 1.891e+00 | 3.187e-01 | YES | | 83 | -8.24974e+03 | 3.541e-01 | 3.466e-02 | CONV | 1.100e+00 | 3.187e-01 | YES | | 84 | -8.24974e+03 | 3.541e-01 | 1.026e-01 | CONV | -1.469e+00 | 3.187e-01 | NO | | 85 | -8.24975e+03 | 1.976e-01 | 5.659e-02 | CONV | 1.017e+00 | 1.593e-01 | YES | | 86 | -8.24975e+03 | 5.387e-01 | 1.593e-01 | BNDRY | 3.360e-01 | 1.593e-01 | YES | | 87 | -8.24975e+03 | 1.334e-01 | 2.824e-02 | CONV | 9.966e-01 | 1.593e-01 | YES | | 88 | -8.24976e+03 | 1.070e-01 | 1.462e-01 | CONV | 1.009e+00 | 1.593e-01 | YES | | 89 | -8.24976e+03 | 1.070e-01 | 3.187e-01 | EXACT | -3.911e-01 | 3.187e-01 | NO | | 90 | -8.24976e+03 | 4.695e-01 | 1.593e-01 | BNDRY | 1.011e+00 | 1.593e-01 | YES | | 91 | -8.24976e+03 | 4.695e-01 | 3.187e-01 | EXACT | -4.858e-01 | 3.187e-01 | NO | | 92 | -8.24976e+03 | 2.509e-01 | 2.394e-03 | CONV | 1.396e+00 | 1.593e-01 | YES | | 93 | -8.24976e+03 | 2.509e-01 | 1.593e-01 | BNDRY | -5.821e+00 | 1.593e-01 | NO | | 94 | -8.24977e+03 | 5.567e-01 | 7.967e-02 | BNDRY | 1.247e+00 | 7.967e-02 | YES | | 95 | -8.24977e+03 | 1.315e-01 | 2.626e-04 | CONV | 9.851e-01 | 1.593e-01 | YES | | 96 | -8.24977e+03 | 1.184e-01 | 1.593e-01 | BNDRY | 1.087e+00 | 1.593e-01 | YES | | 97 | -8.24979e+03 | 1.389e-01 | 3.187e-01 | BNDRY | 1.053e+00 | 3.187e-01 | YES | | 98 | -8.24981e+03 | 2.413e-01 | 6.374e-01 | BNDRY | 1.011e+00 | 6.374e-01 | YES | | 99 | -8.24982e+03 | 2.294e-01 | 1.534e-01 | CONV | 6.901e-01 | 1.275e+00 | YES | |=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 100 | -8.24982e+03 | 7.302e-02 | 2.036e-02 | CONV | 1.013e+00 | 1.275e+00 | YES | | 101 | -8.24984e+03 | 2.291e-01 | 1.275e+00 | BNDRY | 1.148e+00 | 1.275e+00 | YES | | 102 | -8.24985e+03 | 4.335e-02 | 3.434e-02 | CONV | 1.012e+00 | 2.550e+00 | YES | | 103 | -8.24985e+03 | 4.335e-02 | 2.550e+00 | EXACT | -5.461e-01 | 2.550e+00 | NO | | 104 | -8.24985e+03 | 4.335e-02 | 1.275e+00 | EXACT | -7.068e-01 | 1.275e+00 | NO | | 105 | -8.24985e+03 | 1.342e-01 | 6.374e-01 | BNDRY | 1.006e+00 | 6.374e-01 | YES | | 106 | -8.24985e+03 | 1.342e-01 | 1.275e+00 | EXACT | -6.492e-01 | 1.275e+00 | NO | | 107 | -8.24986e+03 | 5.677e-02 | 6.096e-02 | CONV | 9.763e-01 | 6.374e-01 | YES | | 108 | -8.24986e+03 | 6.427e-01 | 6.374e-01 | BNDRY | 1.091e+00 | 6.374e-01 | YES | | 109 | -8.24986e+03 | 7.887e-02 | 2.711e-04 | CONV | 9.314e-01 | 1.275e+00 | YES | | 110 | -8.24986e+03 | 6.857e-02 | 1.257e-01 | CONV | 1.158e+00 | 1.275e+00 | YES | | 111 | -8.24987e+03 | 5.957e-01 | 1.275e+00 | BNDRY | 1.026e+00 | 1.275e+00 | YES | | 112 | -8.24987e+03 | 4.201e-02 | 3.418e-03 | CONV | 1.006e+00 | 2.550e+00 | YES | | 113 | -8.24987e+03 | 4.201e-02 | 2.550e+00 | EXACT | -3.250e-01 | 2.550e+00 | NO | | 114 | -8.24987e+03 | 4.201e-02 | 1.275e+00 | BNDRY | -1.119e+02 | 1.275e+00 | NO | | 115 | -8.24987e+03 | 1.784e-01 | 1.273e-01 | CONV | 1.053e+00 | 6.374e-01 | YES | | 116 | -8.24987e+03 | 1.013e-01 | 2.121e-03 | CONV | 4.568e-01 | 6.374e-01 | YES | | 117 | -8.24987e+03 | 1.082e-02 | 6.319e-04 | CONV | 1.032e+00 | 6.374e-01 | YES | | 118 | -8.24987e+03 | 1.082e-02 | 6.374e-01 | EXACT | -2.854e+00 | 6.374e-01 | NO | | 119 | -8.24987e+03 | 8.461e-02 | 3.187e-01 | BNDRY | 1.126e+00 | 3.187e-01 | YES | |=====================================================================================================| | ITER | FUN VALUE | NORM GRAD | NORM STEP | CG TERM | RHO | TRUST RAD | ACCEPT | |=====================================================================================================| | 120 | -8.24987e+03 | 3.413e-02 | 8.189e-03 | CONV | 8.917e-01 | 6.374e-01 | YES | | 121 | -8.24987e+03 | 1.322e-01 | 6.374e-01 | BNDRY | 1.097e+00 | 6.374e-01 | YES | | 122 | -8.24987e+03 | 2.415e-02 | 5.269e-05 | CONV | 1.056e+00 | 1.275e+00 | YES | | 123 | -8.24987e+03 | 1.411e-01 | 1.275e+00 | BNDRY | 1.263e+00 | 1.275e+00 | YES | | 124 | -8.24987e+03 | 1.192e-02 | 1.134e-03 | CONV | 1.021e+00 | 2.550e+00 | YES | | 125 | -8.24988e+03 | 6.432e-03 | 1.317e+00 | CONV | 1.442e+00 | 2.550e+00 | YES | | 126 | -8.24988e+03 | 3.648e-03 | 1.336e+00 | CONV | 1.448e+00 | 2.550e+00 | YES | | 127 | -8.24988e+03 | 2.015e-03 | 1.366e+00 | CONV | 1.447e+00 | 2.550e+00 | YES | | 128 | -8.24988e+03 | 1.008e-03 | 1.392e+00 | CONV | 1.447e+00 | 2.550e+00 | YES | | 129 | -8.24988e+03 | 5.493e-04 | 1.414e+00 | CONV | 1.446e+00 | 2.550e+00 | YES | | 130 | -8.24988e+03 | 2.889e-04 | 1.434e+00 | CONV | 1.446e+00 | 2.550e+00 | YES | Infinity norm of the final gradient = 2.889e-04 Two norm of the final step = 1.434e+00, TolX = 1.000e-12 Relative infinity norm of the final gradient = 5.753e-07, TolFun = 1.000e-06 EXIT: Local minimum found. GPModel = Nonlinear ARX model with 1 output and 4 inputs Inputs: ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming Outputs: EngineTorque Regressors: Linear regressors in variables EngineTorque, ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming List of all regressors Output function: Gaussian process function using a ARDMatern32 kernel Sample time: 1 seconds Status: Termination condition: Divergent gradient calculation.. Number of iterations: 1, Number of function evaluations: 27 Estimated using NLARX on time domain data "eDataD". Fit to estimation data: 99.97% (simulation focus) FPE: 2.673e-05, MSE: 5.029e-05 More information in model's "Report" property. Model Properties ```
`toc`
```Elapsed time is 920.782575 seconds. ```

Validate the model.

```subplot(211) compare(eDataD,GPModel) % compare against estimation data title('GP Model: Comparison to Estimation Data') subplot(212) compare(vDataD,GPModel) % compare against validation data title('Comparison to Validation Data')```

Support Vector Machine (SVM) Model

As a third option, build an SVM-based nonlinear ARX model.

Unlike the GP- and sigmoid-network-based representations, SVM models do not provide a parallel connection of separate linear and nonlinear functions. Use an `idSupportVectorMachine` object to represent the function.

```% Nonlinear function configuration Kernel = 'Gaussian'; EpsilonMargin = 3.7204e-5; BoxConstraint = 5.2530; SVMFcn = idSupportVectorMachine(Kernel); SVMFcn.BoxConstraint = BoxConstraint; SVMFcn.EpsilonMargin = EpsilonMargin; % Regressor configuration MaxLag = 1; OutputLags = 1:MaxLag; InputLags = 0:MaxLag; Reg2 = linearRegressor([Output, Inputs],{OutputLags,InputLags,InputLags,InputLags,InputLags});```

Use the data (`eDataD`), the regressor set (`Reg2`), and the nonlinear function (SVMFcn) in the `nlarx` command to train the SVM based nonlinear ARX model. Note that unlike the GP- and sigmoid-network-based models, SVM-based models do not support a simulation focus. The model must be trained with a prediction focus. The validation is still based on comparing the simulated response of the trained model.

```nlarxopt.Focus = 'prediction'; nlarxopt.SearchOptions.MaxIterations = 100; tic SVMModel = nlarx(eDataD,Reg2,SVMFcn,nlarxopt)```
```SVMModel = Nonlinear ARX model with 1 output and 4 inputs Inputs: ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming Outputs: EngineTorque Regressors: Linear regressors in variables EngineTorque, ThrottlePosition, WastegateValve, EngineSpeed, SparkTiming List of all regressors Output function: Support Vector Machine function using a Gaussian kernel Sample time: 1 seconds Status: Estimated using NLARX on time domain data "eDataD". Fit to estimation data: 98.44% (prediction focus) FPE: 4.301, MSE: 0.1445 More information in model's "Report" property. Model Properties ```
`toc`
```Elapsed time is 28.535420 seconds. ```

Validate the model.

```subplot(211) compare(eDataD,SVMModel) % compare against estimation data title('SVM Model: Comparison to Estimation Data') subplot(212) compare(vDataD,SVMModel) % compare against validation data title('Comparison to Validation Data (Downsampled)')```

Result Validation

In this example, you create three different types of nonlinear ARX models based on the nature of the function that maps the regressors to the output. All three models show comparably good fits to the validation data.

```clf compare(vDataD,SigmoidNet,GPModel,SVMModel)```

However, note that the training and validation has been performed on the downsampled datasets. If your application requires running the nonlinear ARX model at the original rate of 10 Hz, you can do so in Simulink® by using the rate transition blocks at the input/output ports of the nonlinear ARX Model block.

```% Create an iddata representation of the data in vData to feed the signals to the Simulink % model using an IDDATA Source block. zsim = iddata(vData,InputName=Inputs,OutputName=Output); zsim.Tstart = 0; % Open the simulation model and simulate it using zsim as the source of input data. mdl = 'nlarx_simulator'; open_system(mdl)```

`sim(mdl);`

After the initial transients owing to mismatching initial conditions, all the three models are able to emulate the measured engine torque with high accuracy. The nonlinear ARX models thus created can be used as a data-based proxy for the original high-fidelity system for use in resource-constrained environments.