Hammerstein-Wiener Model of SI Engine Torque Dynamics

This example uses:

This example describes modeling the nonlinear torque dynamics of a spark-ignition (SI) engine as a Hammerstein-Wiener 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 Hammerstein-Wiener models, see Hammerstein-Wiener 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.

Hammerstein-Wiener Model Identification

About Hammerstein-Wiener Models

Hammerstein-Wiener models describe dynamic systems using one or two static nonlinear blocks in series with a linear block. The linear block is a discrete transfer function that represents the dynamic component of the model.

The input and output nonlinearities can be set to flexible basis function expansions such as sigmoid or wavelet networks, or piecewise linear forms. They can also be physics-inspired forms such as saturation or deadzone. You do not have to use nonlinearities at both ends of the linear model. When the model contains only the input nonlinearity, it is called a Hammerstein model. When it contains only the output nonlinearity, it is called a Wiener model.

In this example, you identify Hammerstein-Wiener 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";

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 nonlinear elements to it in a series configuration to create a Hammerstein-Wiener model.

Initial Linear Model

Estimate a discrete-time state-space model. The ssest command allows for searching for the optimal value in a given range. However, for this example, you use a fixed order of seven, which is larger than what the ssest order selector dialog suggests but leads to a better final nonlinear 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. Also specify the horizon, the maximum number of iterations, and the order.

ssOpt = ssestOptions(Focus="simulation");
ssOpt.N4Horizon = [15 0 26];
ssOpt.SearchOptions.MaxIterations = 200;
Order = 7;

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        x6        x7
   x1    0.6307   -0.3216    0.2104   -0.0621   0.08057   0.06688   -0.4296
   x2    0.2057   -0.3095    0.4587  -0.07681    0.1713    0.5015    -1.665
   x3   0.02398   -0.1231    0.3411    0.7639    0.6176   0.07823     -1.12
   x4   0.03332   0.04544   -0.9544   -0.2606    0.3623   -0.5218    0.1301
   x5  -0.06289   -0.3276    0.2222    -0.419    0.2827    0.4606    0.1373
   x6   -0.1268    0.2479   -0.1423   -0.1647    -0.169    0.2188      1.64
   x7  -0.04426  -0.05109  -0.08225  0.005721    0.1281   -0.1702     1.033
 
  B = 
       ThrottlePosi  WastegateVal   EngineSpeed   SparkTiming
   x1     -0.002453      -0.01188     0.0003596      0.003926
   x2      -0.02644       -0.1191      0.003484       0.03948
   x3      -0.02189       -0.1015      0.002956        0.0345
   x4       0.01706        0.0783     -0.002274      -0.02633
   x5       -0.0164      -0.07734      0.002247       0.02601
   x6       0.01384       0.06275     -0.001798        -0.021
   x7      0.001481      0.006872     -0.000192     -0.002306
 
  C = 
                     x1      x2      x3      x4      x5      x6      x7
   EngineTorque   558.3  -174.4   98.83   17.78   61.55   31.32  -375.9
 
  D = 
                 ThrottlePosi  WastegateVal   EngineSpeed   SparkTiming
   EngineTorque        0.1375        0.1187     0.0008919       0.06466
 
  K = 
       EngineTorque
   x1             0
   x2             0
   x3             0
   x4             0
   x5             0
   x6             0
   x7             0
 
Sample time: 1 seconds

Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: yes
   Disturbance component: none
   Number of free coefficients: 88
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                            
Estimated using SSEST on time domain data "eDataD".
Fit to estimation data: 46.39%                     
FPE: 173.9, MSE: 171.7                             
 
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 to only about 50% on the Normalized Root Mean Squared Error (NRMSE) measure. However, as shown next, it provides an excellent starting point for the creation of the Hammerstein-Wiener model.

Hammerstein Model

Identify a model that uses a sigmoid network as input nonlinearity for each of its four inputs. A sigmoid network is a feedforward, single-hidden-layer network represented by the idSigmoidNetwork object.

This nonlinear function is composed of parameters corresponding to the linear term, the offset, and a nonlinear term expressed as a weighted sum of dilated and shifted sigmoid units. 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. For this example, exclude both the linear component and the offset term, and set the number of sigmoid units to 3.

Create the sigmoid-based nonlinear function.

NumberOfUnits = 3; % determined by trial and error
NonlinearFcn = idSigmoidNetwork(NumberOfUnits,false,false)

NonlinearFcn = 
Sigmoid Network

 Nonlinear Function: Sigmoid network with 3 units
 Linear Function: not in use
 Output Offset: not in use

          Inputs: {1×0 cell}
         Outputs: {1×0 cell}
    NonlinearFcn: 'Sigmoid units and their parameters'
       LinearFcn: 'Linear function parameters'
          Offset: 'Offset parameters'

Create a template model (using the idnlhw model constructor) whose linear component is set to linsys and the input nonlinearity to NonlinearFcn.

HammersteinModel = idnlhw(linsys, NonlinearFcn, [])

HammersteinModel =

Hammerstein-Wiener model with 1 output and 4 inputs

Linear transfer function matrix corresponding to the orders:
   nb = [8 8 8 8]
   nf = [7 7 7 7]
   nk = [0 0 0 0]

Input nonlinearities:
  Input 1:  Sigmoid network with 3 units
  Input 2:  Sigmoid network with 3 units
  Input 3:  Sigmoid network with 3 units
  Input 4:  Sigmoid network with 3 units

Output nonlinearity: None
Sample time: 1 seconds

Status:                                                         
Created by direct construction or transformation. Not estimated.
More information in model's "Report" property.

Model Properties

Inspect the linear component of the model.

HammersteinModel.B % numerator of the linear transfer function

ans=1×4 cell array
    {[0.1375 -0.0179 -0.0559 -0.0440 -0.0480 0.0161 0.0401 -0.0222]}    {[0.1187 -0.1175 0.0391 -0.1343 0.1394 -0.1045 0.1484 -0.0888]}    {[8.9186e-04 -0.0025 0.0025 3.1157e-04 -0.0026 0.0042 -0.0053 0.0027]}    {[0.0647 -0.0655 0.0766 -0.0850 0.0083 0.0122 -0.0368 0.0279]}

HammersteinModel.F % denominator of the linear transfer function

ans=1×4 cell array
    {[1 -1.9356 1.8174 -1.7491 1.2668 -0.4497 0.0571 9.3696e-04]}    {[1 -1.9356 1.8174 -1.7491 1.2668 -0.4497 0.0571 9.3696e-04]}    {[1 -1.9356 1.8174 -1.7491 1.2668 -0.4497 0.0571 9.3696e-04]}    {[1 -1.9356 1.8174 -1.7491 1.2668 -0.4497 0.0571 9.3696e-04]}

Furthermore, fix the poles of the model by freezing the coefficients of the denominator ( $F (q)$ ) of the linear model.

HammersteinModel.Ffree = false; % fix the denominator to its current value

For model training, we use the lsqnonlin search method and do not normalize the I/O signals.

% Configure training options
nlhwopt = nlhwOptions('Display','on');
nlhwopt.SearchMethod = 'lsqnonlin';
nlhwopt.SearchOptions.MaxIterations = 300;

% choose to NOT normalize the I/O signals
nlhwopt.Normalize = false;

Estimate the parameters of the template model (HammersteinModel) to fit the data (eDataD) using the nlhw command.

tic
HammersteinModel = nlhw(eDataD,HammersteinModel,nlhwopt)

HammersteinModel =

Hammerstein-Wiener model with 1 output and 4 inputs

Linear transfer function matrix corresponding to the orders:
   nb = [8 8 8 8]
   nf = [7 7 7 7]
   nk = [0 0 0 0]

Input nonlinearities:
  Input 1:  Sigmoid network with 3 units
  Input 2:  Sigmoid network with 3 units
  Input 3:  Sigmoid network with 3 units
  Input 4:  Sigmoid network with 3 units

Output nonlinearity: None
Sample time: 1 seconds

Status:                                                                                        
Termination condition: Maximum number of iterations or number of function evaluations reached..
Number of iterations: 301, Number of function evaluations: 302                                 
                                                                                               
Estimated using NLHW on time domain data "eDataD".                                             
Fit to estimation data: 72.39%                                                                 
FPE: 46.51, MSE: 45.53                                                                         
More information in model's "Report" property.

Model Properties

toc

Elapsed time is 28.918469 seconds.

Perform a 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,HammersteinModel)  % compare against estimation data
title('Hammerstein Model: Comparison to Estimation Data')

subplot(212)
compare(vDataD,HammersteinModel)  % compare against validation data
title('Comparison to Validation Data (Downsampled)')

Hammerstein-Wiener Model

Identify a model that uses piecewise linear functions for both input and output nonlinearities. You have a linear model that you want to extend by adding nonlinear elements. To ensure that the identified model is no worse than the linear one, you freeze the linear block parameters (the $B (q)$ and $F (q)$ polynomials) completely.

PWLin = idPiecewiseLinear(12);  % input nonlinearities
PWLout = idPiecewiseLinear(3); % output nonlinearity
HWModel = idnlhw(linsys, PWLin, PWLout); 

% Fix the linear component completely
HWModel.Bfree = false;
HWModel.Ffree = false;

% Configure estimation options
nlhwopt = nlhwOptions('Display','on');
nlhwopt.SearchMethod = 'lm';
nlhwopt.SearchOptions.MaxIterations = 100;

% Estimate the breakpoint locations of PWLin, PWLout
HWModel = nlhw(eDataD,HWModel,nlhwopt)

HWModel =

Hammerstein-Wiener model with 1 output and 4 inputs

Linear transfer function matrix corresponding to the orders:
   nb = [8 8 8 8]
   nf = [7 7 7 7]
   nk = [0 0 0 0]

Input nonlinearities:
  Input 1:  Piecewise linear with 12 break-points
  Input 2:  Piecewise linear with 12 break-points
  Input 3:  Piecewise linear with 12 break-points
  Input 4:  Piecewise linear with 12 break-points

Output nonlinearity: Piecewise linear with 3 break-points
Sample time: 1 seconds

Status:                                                       
Termination condition: Maximum number of iterations reached.. 
Number of iterations: 100, Number of function evaluations: 527
                                                              
Estimated using NLHW on time domain data "eDataD".            
Fit to estimation data: 72.61%                                
FPE: 46.52, MSE: 44.82                                        
More information in model's "Report" property.

Model Properties

% Validate the model
subplot(211)
compare(eDataD,HWModel)  % compare against estimation data
title('Hammerstein-Wiener Model: Comparison to Estimation Data')

subplot(212)
compare(vDataD,HWModel)  % compare against validation data
title('Comparison to Validation Data')

Result Validation

In this example, you create two models based on different choices of the input and output nonlinearities. The two models show comparably good fits to the validation data.

clf
compare(vDataD,HammersteinModel,HWModel)

However, note that the training and validation has been performed on the downsampled datasets. If your application requires running the 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 Hammerstein-Wiener 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 = 'nlhw_simulator';
open_system(mdl)

sim(mdl);

Both the models are able to emulate the measured engine torque with slightly better results for the HammersteinModel. The models thus created can be used as a data-based proxy for the original high-fidelity system for use in resource-constrained environments.