Wiener–Hammerstein System¶

Note: The example shown in this notebook is taken from the companion book Nonlinear System Identification and Forecasting: Theory and Practice with SysIdentPy.

The description content primarily derives from the benchmark website - Nonlinear Benchmark and associated paper - Wiener-Hammerstein benchmark with process noise. For a detailed description, readers are referred to the linked references.

The nonlinear benchmark website stands as a significant contribution to the system identification and machine learning community. The users are encouraged to explore all the papers referenced on the site.

This benchmark focuses on a Wiener-Hammerstein electronic circuit where process noise plays a significant role in distorting the output signal.

The Wiener-Hammerstein structure is a well-known block-oriented system which contains a static nonlinearity sandwiched between two Linear Time-Invariant (LTI) blocks (Figure 2). This arrangement presents a challenging identification problem due to the presence of these LTI blocks.

Figure 2: the Wiener-Hammerstein system

In Figure 2, the Wiener-Hammerstein system is illustrated with process noise \(e_x(t)\) entering before the static nonlinearity \(f(x)\), sandwiched between LTI blocks represented by \(R(s)\) and \(S(s)\) at the input and output, respectively. Additionally, small, negligible noise sources \(e_u(t)\) and \(e_y(t)\) affect the measurement channels. The measured input and output signals are denoted as \(u_m(t)\) and \(y_m(t)\).

The first LTI block \(R(s)\) is effectively modeled as a third-order lowpass filter. The second LTI subsystem \(S(s)\) is configured as an inverse Chebyshev filter with a stop-band attenuation of \(40 dB\) and a cutoff frequency of \(5 kHz\). Notably, \(S(s)\) includes a transmission zero within the operational frequency range, complicating its inversion.

The static nonlinearity \(f(x)\) is implemented using a diode-resistor network, resulting in saturation nonlinearity. Process noise \(e_x(t)\) is introduced as filtered white Gaussian noise, generated from a discrete-time third-order lowpass Butterworth filter followed by zero-order hold and analog low-pass reconstruction filtering with a cutoff of \(20 kHz\).

Measurement noise sources \(e_u(t)\) and \(e_y(t)\) are minimal compared to \(e_x(t)\). The system's inputs and process noise are generated using an Arbitrary Waveform Generator (AWG), specifically the Agilent/HP E1445A, sampling at \(78125 Hz\), synchronized with an acquisition system (Agilent/HP E1430A) to ensure phase coherence and prevent leakage errors. Buffering between the acquisition cards and the system's inputs and outputs minimizes measurement equipment distortion.

The benchmark provides two standard test signals through the benchmarking website: a random phase multi sine and a sine-sweep signal. Both signals have a \(rms\) value of \(0.71 Vrms\) and cover frequencies from DC to \(15 kHz\) (excluding DC). The sine-sweep spans this frequency range at a rate of \(4.29 MHz/min\). These test sets serve as targets for evaluating the model's performance, emphasizing accurate representation under varied conditions.

The Wiener-Hammerstein benchmark highlights three primary nonlinear system identification challenges:

1. Process Noise: Significant in the system, influencing output fidelity.
2. Static Nonlinearity: Indirectly accessible from measured data, posing identification challenges.
3. Output Dynamics: Complex inversion due to transmission zero presence in \(S(s)\).

The goal of this benchmark is to develop and validate robust models using separate estimation data, ensuring accurate characterization of the Wiener-Hammerstein system's behavior.

Required Packages and Versions¶

This case study was verified with SysIdentPy 0.9.0 on Python 3.12.12 and nonlinear-benchmarks==1.0.1. Install the repository checkout and the official benchmark loader explicitly:

python -m pip install -e .
python -m pip install nonlinear-benchmarks==1.0.1

Use a virtual environment to isolate the optional loader. Numerical results should be recomputed if the environment or model configuration changes.

SysIdentPy configuration¶

In this section, we will demonstrate the application of SysIdentPy to the Wiener-Hammerstein system dataset. The following code will guide you through the process of loading the dataset, configuring the SysIdentPy parameters, and building a model for Wiener-Hammerstein system.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sysidentpy.model_structure_selection import FROLS, AOLS, MetaMSS
from sysidentpy.basis_function import Polynomial, Fourier
from sysidentpy.utils.display_results import results
from sysidentpy.parameter_estimation import (
    LeastSquares,
    BoundedVariableLeastSquares,
    NonNegativeLeastSquares,
    LeastSquaresMinimalResidual,
)

from sysidentpy.metrics import root_mean_squared_error
from sysidentpy.utils.plotting import plot_results

import nonlinear_benchmarks

train_val, test = nonlinear_benchmarks.WienerHammerBenchMark(atleast_2d=True)
x_train, y_train = train_val
x_test, y_test = test

We used the nonlinear_benchmarks package to load the data. The user is referred to the package documentation to check the details of how to use it.

The following plot detail the training and testing data of the experiment.

plot_n = 800

plt.figure(figsize=(15, 4))
plt.plot(x_train[:plot_n])
plt.plot(y_train[:plot_n])
plt.title("Experiment: training data")
plt.legend(["x_train", "y_train"])
plt.show()

plt.figure(figsize=(15, 4))
plt.plot(x_test[:plot_n])
plt.plot(y_test[:plot_n])
plt.title("Experiment: testing data")
plt.legend(["x_test", "y_test"])
plt.show()

The external benchmark provides useful context for the experiment. A direct ranking, however, requires the same split, initialization window and normalization; those conditions are made explicit below before any comparison is drawn.

State of the art results presented in the benchmarking paper. In this section we are only working with the Wiener-Hammerstein results, which are presented in the \(W-H\) column.

Results¶

We will start with a basic configuration of FROLS using a polynomial basis function with degree equal 2. The xlag and ylag are set to \(7\) in this first example. Because the dataset is considerably large, we will start with n_info_values=50. This means the FROLS algorithm will not include all regressors when calculating the information criteria used to determine the model order. While this approach might result in a suboptimal model, it is a reasonable starting point for our first attempt.

n = test.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=7,
    ylag=7,
    basis_function=basis_function,
    estimator=LeastSquares(unbiased=False),
    n_info_values=50,
)

model.fit(X=x_train, y=y_train)
if model.max_lag > n:
    raise ValueError("The model lag exceeds the benchmark initialization window.")
start = n - model.max_lag
yhat = model.predict(X=x_test[start:], y=y_test[start:n])
yhat = yhat[model.max_lag :]
rmse = root_mean_squared_error(y_test[n:], yhat)
nrmse = rmse / np.std(y_test[n:])
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=1000,
    title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(nrmse, 4)}",
)

The first configuration gives RMSE \(0.020007\) and NRMSE \(0.082029\). We started with xlag=ylag=7 to establish a compact baseline. The benchmarking paper uses longer memories in some models, so the next configuration sets xlag=ylag=10.

x_train, y_train = train_val
x_test, y_test = test

n = test.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=10,
    ylag=10,
    basis_function=basis_function,
    estimator=LeastSquares(unbiased=False),
    n_info_values=50,
)

model.fit(X=x_train, y=y_train)
if model.max_lag > n:
    raise ValueError("The model lag exceeds the benchmark initialization window.")
start = n - model.max_lag
yhat = model.predict(X=x_test[start:], y=y_test[start:n])
yhat = yhat[model.max_lag :]
rmse = root_mean_squared_error(y_test[n:], yhat)
nrmse = rmse / np.std(y_test[n:])
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=1000,
    title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(nrmse, 4)}",
)

The 10-lag configuration improves the result to RMSE \(0.015202\) and NRMSE \(0.062328\). For now, we are not optimizing model complexity. The information- criterion trace nevertheless shows that the 50-regressor model changes little after many of the later additions.

plt.plot(model.info_values)

So, what happens if we set a model with half of the regressors?

x_train, y_train = train_val
x_test, y_test = test

n = test.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=10,
    ylag=10,
    basis_function=basis_function,
    estimator=LeastSquares(unbiased=False),
    n_info_values=50,
    n_terms=25,
    order_selection=False,
)

model.fit(X=x_train, y=y_train)
if model.max_lag > n:
    raise ValueError("The model lag exceeds the benchmark initialization window.")
start = n - model.max_lag
yhat = model.predict(X=x_test[start:], y=y_test[start:n])
yhat = yhat[model.max_lag :]
rmse = root_mean_squared_error(y_test[n:], yhat)
nrmse = rmse / np.std(y_test[n:])
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=1000,
    title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(nrmse, 4)}",
)

The fixed 25-term model gives RMSE \(0.018809\) and NRMSE \(0.077117\). It is more compact than the automatically selected 50-term model, at the cost of a larger error. The three current results are summarized below.

The historical figures stored four decimal places, and the reproduced values remain the same at that precision. Published tables may use a different normalization or split, so they are not treated as a direct ranking. Users who want to investigate the deep state-space alternatives can explore the deepsysid package.

This basic configuration can serve as a starting point for users to develop even better models using SysIdentPy. Give it a try!