Silver Box System¶

Reproducibility¶

This tutorial was verified with SysIdentPy 0.9.0 on Python 3.12.12 and additionally requires nonlinear-benchmarks==1.0.1. Numerical results must be recomputed when the environment or model configuration changes.

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 mainly derives (copy and paste) from the associated paper - Three free data sets for development and benchmarking in nonlinear system identification. For a detailed description, readers are referred to the linked reference.

The Silverbox system can be seen as an electronic implementation of the Duffing oscillator. It is build as a 2^nd order linear time-invariant system with a 3^rd degree polynomial static nonlinearity around it in feedback. This type of dynamics are, for instance, often encountered in mechanical systems Nonlinear Benchmark - Silverbox.

In this case study, we will create a NARX model for the Silver box benchmark. The Silver box represents a simplified version of mechanical oscillating processes, which are a critical category of nonlinear dynamic systems. Examples include vehicle suspensions, where shock absorbers and progressive springs play vital roles. The data generated by the Silver box provides a simplified representation of such combined components. The electrical circuit generating this data closely approximates, but does not perfectly match, the idealized models described below.

As described in the original paper, the system was excited using a general waveform generator (HPE1445A). The input signal begins as a discrete-time signal \(r(k)\), which is converted to an analog signal \(r_c(t)\) using zero-order-hold reconstruction. The actual excitation signal \(u_0(t)\) is then obtained by passing \(r_c(t)\) through an analog low-pass filter \(G(p)\) to eliminate high-frequency content around multiples of the sampling frequency. Here, \(p\) denotes the differentiation operator. Thus, the input is given by:

The input and output signals were measured using HP1430A data acquisition cards, with synchronized clocks for the acquisition and generator cards. The sampling frequency was:

The silver box uses analog electrical circuitry to generate data representing a nonlinear mechanical resonating system with a moving mass \(m\), viscous damping \(d\), and a nonlinear spring \(k(y)\). The electrical circuit is designed to relate the displacement \(y(t)\) (the output) to the force \(u(t)\) (the input) by the following differential equation:

The nonlinear progressive spring is described by a static, position-dependent stiffness:

The signal-to-noise ratio is sufficiently high to model the system without accounting for measurement noise. However, measurement noise can be included by replacing \(y(t)\) with the artificial variable \(x(t)\) in the equation above, and introducing disturbances \(w(t)\) and \(e(t)\) as follows:

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 Silver box dataset. The following code will guide you through the process of loading the dataset, configuring the SysIdentPy parameters, and building a model for mentioned system.

import numpy as np
import matplotlib.pyplot as plt

from sysidentpy.model_structure_selection import FROLS
from sysidentpy.basis_function import Polynomial, Fourier
from sysidentpy.parameter_estimation import LeastSquares
from sysidentpy.metrics import root_mean_squared_error
from sysidentpy.utils.plotting import plot_results

import nonlinear_benchmarks

train_val, test = nonlinear_benchmarks.Silverbox(atleast_2d=True)

x_train, y_train = train_val.u, train_val.y
test_multisine, test_arrow_full, test_arrow_no_extrapolation = test
x_test, y_test = test_multisine.u, test_multisine.y

n = test_multisine.state_initialization_window_length

We used the nonlinear_benchmarks package to load the data. The user is referred to the package documentation - GerbenBeintema/nonlinear_benchmarks: The official data load for http://www.nonlinearbenchmark.org/ (github.com) to check the details of how to use it.

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

plt.plot(x_train)
plt.plot(y_train, alpha=0.3)
plt.title("Experiment 1: training data")
plt.show()

plt.plot(x_test)
plt.plot(y_test, alpha=0.3)
plt.title("Experiment 1: testing data")
plt.show()

plt.plot(test_arrow_full.u)
plt.plot(test_arrow_full.y, alpha=0.3)
plt.title("Experiment 2: training data")
plt.show()

plt.plot(test_arrow_no_extrapolation.u)
plt.plot(test_arrow_no_extrapolation.y, alpha=0.2)
plt.title("Experiment 2: testing data")
plt.show()

Important Note

The goal of this benchmark is to develop a model that outperforms the state-of-the-art (SOTA) model presented in the benchmarking paper. However, the results in the paper differ from those provided in the GitHub repository.

It appears that the values shown in the paper actually represent the training time, not the error metrics. I will contact the authors to confirm this information. According to the Nonlinear Benchmark website, the information is as follows:

where the values in the "Training time" column matches the ones presented as error metrics in the paper.

While we await confirmation of the correct values for this benchmark, we will demonstrate the performance of SysIdentPy. However, we will refrain from making any comparisons or attempting to improve the model at this stage.

Results¶

We will start (as we did in every other case study) with a basic configuration of FROLS using a polynomial basis function with degree equal to 2. The xlag and ylag are set to \(7\) in this first example. Because the dataset is large, we use err_tol and set n_terms=40; the search stops when the error-reduction tolerance is reached or 40 regressors have been tested. While this approach may result in a suboptimal model, it is a reasonable starting point. There are three tests: multisine, arrow (full), and arrow (without extrapolation).

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=7,
    ylag=7,
    basis_function=basis_function,
    estimator=LeastSquares(),
    err_tol=0.999,
    n_terms=40,
    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:])
rmse_mv = 1000 * rmse
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")
plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=30000,
    figsize=(15, 4),
    title=f"Multisine. Model -> RMSE (x1000) mv: {round(rmse_mv, 4)}",
)

plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=300,
    figsize=(15, 4),
    title=f"Multisine. Model -> RMSE (x1000) mv: {round(rmse_mv, 4)}",
)

x_train, y_train = train_val.u, train_val.y
test_multisine, test_arrow_full, test_arrow_no_extrapolation = test
x_test, y_test = test_arrow_full.u, test_arrow_full.y

n = test_arrow_full.state_initialization_window_length

basis_function = Polynomial(degree=3)
model = FROLS(
    xlag=14,
    ylag=14,
    basis_function=basis_function,
    estimator=LeastSquares(),
    err_tol=0.9999,
    n_terms=80,
    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:])
rmse_mv = 1000 * rmse

print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")

plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=30000,
    figsize=(15, 4),
    title=f"Arrow (full). Model -> RMSE (x1000) mv: {round(rmse_mv, 4)}",
)

plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=300,
    figsize=(15, 4),
    title=f"Arrow (full). Model -> RMSE (x1000) mv: {round(rmse_mv, 4)}",
)

x_train, y_train = train_val.u, train_val.y
test_multisine, test_arrow_full, test_arrow_no_extrapolation = test
x_test, y_test = test_arrow_no_extrapolation.u, test_arrow_no_extrapolation.y

n = test_arrow_no_extrapolation.state_initialization_window_length

basis_function = Polynomial(degree=3)
model = FROLS(
    xlag=14,
    ylag=14,
    basis_function=basis_function,
    estimator=LeastSquares(),
    err_tol=0.9999,
    n_terms=40,
    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:])
rmse_mv = 1000 * rmse
print(f"RMSE: {rmse:.6f}; NRMSE: {nrmse:.6f}")

plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=30000,
    figsize=(15, 4),
    title=f"Arrow (no extrapolation). Model -> RMSE (x1000) mv: {round(rmse_mv, 4)}",
)

plot_results(
    y=y_test[n:],
    yhat=yhat,
    n=300,
    figsize=(15, 4),
    title=f"Free Run simulation. Model -> RMSE (x1000) mv: {round(rmse_mv, 4)}",
)

The current free-run results, evaluated after the 50-sample initialization window supplied by the loader, are:

The complete arrow test is harder than the version without extrapolation. The metric excludes the initialization window once; it does not remove model.max_lag a second time. The external deep state-space results remain valuable context, but are not ranked directly against this table without first matching the split, initialization and normalization.