10. Case Studies

M4 Dataset¶

The M4 dataset is a well known resource for time series forecasting, offering a wide range of data series used to test and improve forecasting methods. Created for the M4 competition organized by Spyros Makridakis, this dataset has driven many advancements in forecasting techniques.

The M4 dataset includes 100,000 time series from various fields such as demographics, finance, industry, macroeconomics, and microeconomics, which were selected randomly from the ForeDeCk database. The series come in different frequencies (yearly, quarterly, monthly, weekly, daily, and hourly), making it a comprehensive collection for testing forecasting methods.

In this case study, we will focus on the hourly subset of the M4 dataset. This subset consists of time series data recorded hourly, providing a detailed and high-frequency look at changes over time. Hourly data presents unique challenges due to its granularity and the potential for capturing short-term fluctuations and patterns.

The M4 dataset provides a standard benchmark to compare different forecasting methods, allowing researchers and practitioners to evaluate their models consistently. With series from various domains and frequencies, the M4 dataset represents real-world forecasting challenges, making it valuable for developing robust forecasting techniques. The competition and the dataset itself have led to the creation of new algorithms and methods, significantly improving forecasting accuracy and reliability.

We will present a end to end walkthrough using the M4 hourly dataset to demonstrate the capabilities of SysIdentPy. SysIdentPy offers a range of tools and techniques designed to effectively handle the complexities of time series data, but we will focus on fast and easy setup for this case. We will cover model selection and evaluation metrics specific to the hourly dataset.

By the end of this case study, you will have a solid understanding of how to use SysIdentPy for forecasting with the M4 hourly dataset, preparing you to tackle similar forecasting challenges in real-world scenarios.

Required Packages and Versions¶

To ensure that you can replicate this case study, it is essential to use specific versions of the required packages. Below is a list of the packages along with their respective versions needed for running the case studies effectively.

To install all the required packages, you can create a requirements.txt file with the following content:

sysidentpy==0.4.0
datasetsforecast==0.0.8
pandas==2.2.2
numpy==1.26.0
matplotlib==3.8.4
s3fs==2024.6.1

Then, install the packages using:

pip install -r requirements.txt

• Ensure that you use a virtual environment to avoid conflicts between package versions.
• Versions specified are based on compatibility with the code examples provided. If you are using different versions, some adjustments in the code might be necessary.

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 warnings
import numpy as np
import pandas as pd
from pandas.errors import SettingWithCopyWarning
import matplotlib.pyplot as plt

from sysidentpy.model_structure_selection import FROLS, AOLS
from sysidentpy.basis_function._basis_function import Polynomial
from sysidentpy.parameter_estimation import LeastSquares
from sysidentpy.metrics import root_relative_squared_error, symmetric_mean_absolute_percentage_error
from sysidentpy.utils.plotting import plot_results

from datasetsforecast.m4 import M4, M4Evaluation

warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=UserWarning)
warnings.simplefilter(action='ignore', category=SettingWithCopyWarning)

train = pd.read_csv('https://auto-arima-results.s3.amazonaws.com/M4-Hourly.csv')
test = pd.read_csv('https://auto-arima-results.s3.amazonaws.com/M4-Hourly-test.csv').rename(columns={'y': 'y_test'})

The following plots provide a visualization of the training data for a small subset of the time series. The plot shows the raw data, giving you an insight into the patterns and behaviors inherent in each series.

By observing the data, you can get a sense of the variety and complexity of the time series we are working with. The plots can reveal important characteristics such as trends, seasonal patterns, and potential anomalies within the time series. Understanding these elements is crucial for the development of accurate forecasting models.

However, when dealing with a large number of different time series, it is common to start with broad assumptions rather than detailed individual analysis. In this context, we will adopt a similar approach. Instead of going into the specifics of each dataset, we will make some general assumptions and see how SysIdentPy handles them.

This approach provides a practical starting point, demonstrating how SysIdentPy can manage different types of time series data without too much work. As you become more familiar with the tool, you can refine your models with more detailed insights. For now, let's focus on using SysIdentPy to create the forecasts based on these initial assumptions.

Our first assumption is that there is a 24-hour seasonal pattern in the series. By examining the plots below, this seems reasonable. Therefore, we'll begin building our models with ylag=24.

ax = train[train["unique_id"]=="H10"].reset_index(drop=True)["y"].plot(figsize=(15, 2), title="H10")
xcoords = [a for a in range(24, 24*30, 24)]

for xc in xcoords:
    plt.axvline(x=xc, color='red', linestyle='--', alpha=0.5)

Lets check build a model for the H20 group before we extrapolate the settings for every group. Because there are no input features, we will be using a NAR model type in SysIdentPy. To keep things simple and fast, we will start with Polynomial basis function with degree \(1\).

unique_id = "H20"
y_id = train[train["unique_id"]==unique_id]["y"].values.reshape(-1, 1)
y_val = test[test["unique_id"]==unique_id]["y_test"].values.reshape(-1, 1)

basis_function = Polynomial(degree=1)
model = FROLS(
    order_selection=True,
    ylag=24,
    estimator=LeastSquares(),
    basis_function=basis_function,
    model_type="NAR",
)

model.fit(y=y_id)
y_val = np.concatenate([y_id[-model.max_lag :], y_val])
y_hat = model.predict(y=y_val, forecast_horizon=48)
smape = symmetric_mean_absolute_percentage_error(y_val[model.max_lag::], y_hat[model.max_lag::])

plot_results(y=y_val[model.max_lag :], yhat=y_hat[model.max_lag :], n=30000, figsize=(15, 4), title=f"Group: {unique_id} - SMAPE {round(smape, 4)}")

Probably, the result are not optimal and will not work for every group. However, let's check how this setting performs against the winner model  M4 time series competition: the Exponential Smoothing with Recurrent Neural Networks (ESRNN).

esrnn_url = 'https://github.com/Nixtla/m4-forecasts/raw/master/forecasts/submission-118.zip'
esrnn_forecasts = M4Evaluation.load_benchmark('data', 'Hourly', esrnn_url)
esrnn_evaluation = M4Evaluation.evaluate('data', 'Hourly', esrnn_forecasts)

esrnn_evaluation

The following code took only 49 seconds to run on my machine (AMD Ryzen 5 5600x processor, 32GB RAM at 3600MHz). Because of its efficiency, I didn't create a parallel version. By the end of this use case, you will see how SysIdentPy can be both fast and effective, delivering good results without too much optimization.

r = []
ds_test = list(range(701, 749))
for u_id, data in train.groupby(by=["unique_id"], observed=True):
    y_id = data["y"].values.reshape(-1, 1)
    basis_function = Polynomial(degree=1)
    model = FROLS(
        ylag=24,
        estimator=LeastSquares(),
        basis_function=basis_function,
        model_type="NAR",
        n_info_values=25,
    )
    try:
        model.fit(y=y_id)
        y_val = y_id[-model.max_lag :].reshape(-1, 1)
        y_hat = model.predict(y=y_val, forecast_horizon=48)
        r.append(
            [
                u_id*len(y_hat[model.max_lag::]),
                ds_test,
                y_hat[model.max_lag::].ravel()
            ]
        )
    except Exception:
        print(f"Problem with {u_id}")

results_1 = pd.DataFrame(r, columns=["unique_id", "ds", "NARMAX_1"]).explode(['unique_id', 'ds', 'NARMAX_1'])
results_1["NARMAX_1"] = results_1["NARMAX_1"].astype(float)#.clip(lower=10)
pivot_df = results_1.pivot(index='unique_id', columns='ds', values='NARMAX_1')
results = pivot_df.to_numpy()

M4Evaluation.evaluate('data', 'Hourly', results)

The initial results are reasonable, but they don't quite match the performance of ESRNN. These results are based solely on our first assumption. To better understand the performance, let’s examine the groups with the worst results.

The following plot illustrates two such groups, H147 and H136. Both exhibit a 24-hour seasonal pattern.

However, a closer look reveals an additional insight: in addition to the daily pattern, these series also show a weekly pattern. Observe how the data looks like when we split the series into weekly segments.

xcoords = list(range(0, 168*5, 168))
filtered_train = train[train["unique_id"] == "H147"].reset_index(drop=True)

fig, ax = plt.subplots(figsize=(10, 1.5 * len(xcoords[1:])))
for i, start in enumerate(xcoords[:-1]):
    end = xcoords[i + 1]
    ax = fig.add_subplot(len(xcoords[1:]), 1, i + 1)
    filtered_train["y"].iloc[start:end].plot(ax=ax)
    ax.set_title(f'H147 -> Slice {i+1}: Hour {start} to {end-1}')

plt.tight_layout()
plt.show()

Therefore, we will build models setting ylag=168.

Note that this is a very high number for lags, so be careful if you want to try it with higher polynomial degrees because the time to run the models can increase significantly. I tried some configurations with polynomial degree equal to 2 and only took \(6\) minutes to run (even less, using AOLS), without making the code run in parallel. As you can see, SysIdentPy can be very fast and you can make it faster by applying parallelization.

# this took 2min to run on my computer.
r = []
ds_test = list(range(701, 749))
for u_id, data in train.groupby(by=["unique_id"], observed=True):
    y_id = data["y"].values.reshape(-1, 1)
    basis_function = Polynomial(degree=1)
    model = FROLS(
            ylag=168,
            estimator=LeastSquares(),
            basis_function=basis_function,
            model_type="NAR",
        )
    try:
        model.fit(y=y_id)
        y_val = y_id[-model.max_lag :].reshape(-1, 1)
        y_hat = model.predict(y=y_val, forecast_horizon=48)
        r.append(
            [
                u_id*len(y_hat[model.max_lag::]),
                ds_test,
                y_hat[model.max_lag::].ravel()
            ]
        )
    except Exception:
        print(f"Problema no id {u_id}")

results_1 = pd.DataFrame(r, columns=["unique_id", "ds", "NARMAX_1"]).explode(['unique_id', 'ds', 'NARMAX_1'])
results_1["NARMAX_1"] = results_1["NARMAX_1"].astype(float)#.clip(lower=10)
pivot_df = results_1.pivot(index='unique_id', columns='ds', values='NARMAX_1')
results = pivot_df.to_numpy()
M4Evaluation.evaluate('data', 'Hourly', results)

Now, the results are much closer to those of the ESRNN model! While the Symmetric Mean Absolute Percentage Error (SMAPE) is slightly worse, the Mean Absolute Scaled Error (MASE) is better when comparing against ESRNN, leading to a very similar Overall Weighted Average (OWA) metric. Remarkably, these results are achieved using only simple AR models. Next, let's see if the AOLS method can provide even better results.

r = []
ds_test = list(range(701, 749))
for u_id, data in train.groupby(by=["unique_id"], observed=True):
    y_id = data["y"].values.reshape(-1, 1)
    basis_function = Polynomial(degree=1)
    model = AOLS(
        ylag=168,
        basis_function=basis_function,
        model_type="NAR",
        # due to high lag settings, k was increased to 6 as an initial guess
        k=6,
    )
    try:
        model.fit(y=y_id)
        y_val = y_id[-model.max_lag :].reshape(-1, 1)
        y_hat = model.predict(y=y_val, forecast_horizon=48)
        r.append(
            [
                u_id*len(y_hat[model.max_lag::]),
                ds_test,
                y_hat[model.max_lag::].ravel()
            ]
        )
    except Exception:
        print(f"Problema no id {u_id}")

results_1 = pd.DataFrame(r, columns=["unique_id", "ds", "NARMAX_1"]).explode(['unique_id', 'ds', 'NARMAX_1'])
results_1["NARMAX_1"] = results_1["NARMAX_1"].astype(float)#.clip(lower=10)
pivot_df = results_1.pivot(index='unique_id', columns='ds', values='NARMAX_1')
results = pivot_df.to_numpy()
M4Evaluation.evaluate('data', 'Hourly', results)

The Overall Weighted Average (OWA) is even better than that of the ESRNN model! Additionally, the AOLS method was incredibly efficient, taking only 6 seconds to run. This combination of high performance and rapid execution makes AOLS a compelling alternative for time series forecasting in cases with multiple series.

Before we finish, let's verify how the performance of the H147 model has improved with the ylag=168 setting.

Based on the M4 benchmark paper, we could also clip the predictions lower than 10 to 10 and the results would be slightly better. But this is left to the user.

We could achieve even better performance with some fine-tuning of the model configuration. However, I’ll leave exploring these alternative adjustments as an exercise for the user. However, keep in mind that experimenting with different settings does not always guarantee improved results. A deeper theoretical knowledge can often lead you to better configurations and, hence, better results.

Coupled Eletric Device¶

The CE8 coupled electric drives dataset presents a compelling use case for demonstrating the performance of SysIdentPy. This system involves two electric motors driving a pulley with a flexible belt, creating a dynamic environment ideal for testing system identification tools.

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.

System Overview¶

The CE8 system, illustrated in Figure 1, features: - Two Electric Motors: These motors independently control the tension and speed of the belt, providing symmetrical control around zero. This enables both clockwise and counterclockwise movements. - Pulley Mechanism: The pulley is supported by a spring, introducing a lightly damped dynamic mode that adds complexity to the system. - Speed Control Focus: The primary focus is on the speed control system. The pulley’s angular speed is measured using a pulse counter, which is insensitive to the direction of the velocity.

Figure 1. CE8 system design.

Sensor and Filtering¶

The measurement process involves: - Pulse Counter: This sensor measures the angular speed of the pulley without regard to the direction. - Analogue Low Pass Filtering: This reduces high-frequency noise, followed by anti-aliasing filtering to prepare the signal for digital processing. The dynamic effects are mainly influenced by the electric drive time constants and the spring, with the low pass filtering having a minimal impact on the output.

SOTA Results¶

SysIdentPy can be used to build robust models for identifying and modeling the complex dynamics of the CE8 system. The performance will be compared against a benchmark provided by Max D. Champneys, Gerben I. Beintema, Roland Tóth, Maarten Schoukens, and Timothy J. Rogers - Baselines for Nonlinear Benchmarks, Workshop on Nonlinear System Identification Benchmarks, 2024.

The benchmark evaluate the average metric between the two experiments. That's why the SOTA method do not have the better metric for test 1, but it is still the best overall. The goal of this case study is not only to showcase the robustness of SysIdentPy but also provides valuable insights into its practical applications in real-world dynamic systems.

Required Packages and Versions¶

To ensure that you can replicate this case study, it is essential to use specific versions of the required packages. Below is a list of the packages along with their respective versions needed for running the case studies effectively.

To install all the required packages, you can create a requirements.txt file with the following content:

sysidentpy==0.4.0
pandas==2.2.2
numpy==1.26.0
matplotlib==3.8.4
nonlinear_benchmarks==0.1.2

Then, install the packages using:

pip install -r requirements.txt

• Ensure that you use a virtual environment to avoid conflicts between package versions.
• Versions specified are based on compatibility with the code examples provided. If you are using different versions, some adjustments in the code might be necessary.

SysIdentPy configuration¶

In this section, we will demonstrate the application of SysIdentPy to the CE8 coupled electric drives dataset. This example showcases the robust performance of SysIdentPy in modeling and identifying complex dynamic systems. The following code will guide you through the process of loading the dataset, configuring the SysIdentPy parameters, and building a model for CE8 system.

This practical example will help users understand how to effectively utilize SysIdentPy for their own system identification tasks, leveraging its advanced features to handle the complexities of real-world dynamic systems. Let's dive into the code and explore the capabilities of SysIdentPy.

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

from sysidentpy.model_structure_selection import FROLS
from sysidentpy.basis_function import Polynomial, Fourier
from sysidentpy.utils.display_results import results
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.CED(atleast_2d=True)
data_train_1, data_train_2 = train_val
data_test_1, data_test_2   = 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 both experiments. Here we are trying to get two models, one for each experiment, that have a better performance than the mentioned baselines.

plt.plot(data_train_1.u)
plt.plot(data_train_1.y)
plt.title("Experiment 1: training data")
plt.show()

plt.plot(data_test_1.u)
plt.plot(data_test_1.y)
plt.title("Experiment 1: testing data")
plt.show()

plt.plot(data_train_2.u)
plt.plot(data_train_2.y)
plt.title("Experiment 2: training data")
plt.show()

plt.plot(data_test_2.u)
plt.plot(data_test_2.y)
plt.title("Experiment 2: testing data")
plt.show()

Results¶

First, we will set the exactly same configuration to built models for both experiments. We can have better models by optimizing the configurations individually, but we will start simple.

A basic configuration of FROLS using a polynomial basis function with degree equal 2 is defined. The information criteria will be the default one, the aic. The xlag and ylag are set to \(7\) in this first example.

Model for experiment 1:

y_train = data_train_1.y
y_test = data_test_1.y
x_train = data_train_1.u
x_test = data_test_1.u

n = data_test_1.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=7,
    ylag=7,
    basis_function=basis_function,
    estimator=LeastSquares(),
    info_criteria="aic",
    n_info_values=120
)

model.fit(X=x_train, y=y_train)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])
rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title=f"Free Run simulation. Model 1 -> RMSE: {round(rmse, 4)}")

Model for experiment 2:

y_train = data_train_2.y
y_test = data_test_2.y
x_train = data_train_2.u
x_test = data_test_2.u

n = data_test_2.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=7,
    ylag=7,
    basis_function=basis_function,
    estimator=LeastSquares(),
    info_criteria="aic",
    n_info_values=120
)

model.fit(X=x_train, y=y_train)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])
rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title=f"Free Run simulation. Model 2 -> RMSE: {round(rmse, 4)}")

The first configuration for experiment 1 is already better than the LTI ARX, LTI SS, GRU, LSTM, MLP NARX, MLP FIR, OLSTM, and the SOTA models shown in the benchmark table. Better than 8 out 11 models shown in the benchmark. For experiment 2, its better than LTI ARX, LTI SS, GRU, RNN, LSTM, OLSTM, and pNARX (7 out 11). It's a good start, but let's check if the performance improves if we set a higher lag for both xlag and ylag.

The average metric is \((0.1131 + 0.1059)/2 = 0.1095\), which is very good, but worse than the SOTA (\(0.0945\)). We will now increase the lags for x and y to check if we get a better model. Before increasing the lags, the information criteria is shown:

xaxis = np.arange(1, model.n_info_values + 1)
plt.plot(xaxis, model.info_values)
plt.xlabel("n_terms")
plt.ylabel("Information Criteria")

It can be observed that after 22 regressors, adding new regressors do not improve the model performance (considering the configuration defined for that model). Because we want to try models with higher lags and higher nonlinearity degree, the stopping criteria will be changed to err_tol instead of information criteria. This will made the algorithm runs considerably faster.

# experiment 1
y_train = data_train_1.y
y_test = data_test_1.y
x_train = data_train_1.u
x_test = data_test_1.u

n = data_test_1.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=14,
    ylag=14,
    basis_function=basis_function,
    estimator=LeastSquares(),
    err_tol=0.9996,
    n_terms=22,
    order_selection=False
)

model.fit(X=x_train, y=y_train)
print(model.final_model.shape, model.err.sum())
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])

rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])

plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title=f"Free Run simulation. Model 1 -> RMSE: {round(rmse, 4)}")

# experiment 2
y_train = data_train_2.y
y_test = data_test_2.y
x_train = data_train_2.u
x_test = data_test_2.u

n = data_test_2.state_initialization_window_length

basis_function = Polynomial(degree=2)
model = FROLS(
    xlag=14,
    ylag=14,
    basis_function=basis_function,
    estimator=LeastSquares(),
    info_criteria="aicc",
    err_tol=0.9996,
    n_terms=22,
    order_selection=False
)

model.fit(X=x_train, y=y_train)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])

rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])

plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title=f"Free Run simulation. Model 2 -> RMSE: {round(rmse, 4)}")

In the first experiment, the model showed a slight improvement, while the performance of the second experiment experienced a minor decline. Increasing the lag settings with these configurations did not result in significant changes. Therefore, let's set the polynomial degree to \(3\) and increase the number of terms to build the model to n_terms=40 if the err_tol is not reached. It's important to note that these values are chosen empirically. We could also adjust the parameter estimation technique, the err_tol, the model structure selection algorithm, and the basis function, among other factors. Users are encouraged to employ hyperparameter tuning techniques to find the optimal combinations of hyperparameters.

# experiment 1
y_train = data_train_1.y
y_test = data_test_1.y
x_train = data_train_1.u
x_test = data_test_1.u

n = data_test_1.state_initialization_window_length

basis_function = Polynomial(degree=3)
model = FROLS(
    xlag=14,
    ylag=14,
    basis_function=basis_function,
    estimator=LeastSquares(),
    err_tol=0.9996,
    n_terms=40,
    order_selection=False
)

model.fit(X=x_train, y=y_train)
print(model.final_model.shape, model.err.sum())
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])

rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])

plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title=f"Free Run simulation. Model 1 -> RMSE: {round(rmse, 4)}")

# experiment 2
y_train = data_train_2.y
y_test = data_test_2.y
x_train = data_train_2.u
x_test = data_test_2.u

n = data_test_2.state_initialization_window_length

basis_function = Polynomial(degree=3)
model = FROLS(
    xlag=14,
    ylag=14,
    basis_function=basis_function,
    estimator=LeastSquares(),
    info_criteria="aicc",
    err_tol=0.9996,
    n_terms=40,
    order_selection=False
)

model.fit(X=x_train, y=y_train)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])

rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])

plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title=f"Free Run simulation. Model 2 -> RMSE: {round(rmse, 4)}")

As shown in the plot, we have surpassed the state-of-the-art (SOTA) results with an average metric of \((0.0969 + 0.0731)/2 = 0.0849\). Additionally, the metric for the first experiment matches the best model in the benchmark, and the metric for the second experiment slightly exceeds the benchmark's best model. Using the same configuration for both models, we achieved the best overall results!

Wiener-Hammerstein¶

The description content primarily derives from the benchmark website and associated paper. 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 multisine and a sine-sweep signal. Both signals have an \(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¶

To ensure that you can replicate this case study, it is essential to use specific versions of the required packages. Below is a list of the packages along with their respective versions needed for running the case studies effectively.

To install all the required packages, you can create a requirements.txt file with the following content:

sysidentpy==0.4.0
pandas==2.2.2
numpy==1.26.0
matplotlib==3.8.4
nonlinear_benchmarks==0.1.2

Then, install the packages using:

pip install -r requirements.txt

• Ensure that you use a virtual environment to avoid conflicts between package versions.
• Versions specified are based on compatibility with the code examples provided. If you are using different versions, some adjustments in the code might be necessary.

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 goal of this benchmark it to get a model that have a better performance than the SOTA model provided in the benchmarking paper.

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 sub-optimal 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)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])
rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])
rmse_sota = rmse/y_test.std()
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=1000, title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(rmse_sota, 4)}")

The first configuration is already better than the SOTA models shown in the benchmark table! We started using xlag=ylag=7 to have an idea of how well SysIdentPy would handle this dataset, but the results are pretty good already! However, the benchmarking paper indicates that they used higher lags for their models. Let's check what happens if we set 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)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])
rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])
rmse_sota = rmse/y_test.std()
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=1000, title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(rmse_sota, 4)}")

The performance is even better now! For now, we are not worried about the model complexity (even in this case where we are comparing to a deep state neural network...). However, if we check the model order and the AIC plot, we see that the model have 50 regressors , but the AIC values do not change much after each added regression.

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)
y_test = np.concatenate([y_train[-model.max_lag:], y_test])
x_test = np.concatenate([x_train[-model.max_lag:], x_test])
yhat = model.predict(X=x_test, y=y_test[:model.max_lag, :])
rmse = root_mean_squared_error(y_test[model.max_lag + n :], yhat[model.max_lag + n:])
rmse_sota = rmse/y_test.std()
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=1000, title=f"SysIdentPy -> RMSE: {round(rmse, 4)}, NRMSE: {round(rmse_sota, 4)}")

As shown in the figure above, the results still outperform the SOTA models presented in the benchmarking paper. The SOTA results from the paper could likely be improved as well. Users are encouraged to explore the deepsysid package, which can be used to build deep state neural networks.

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

Air Passenger Demand Forecasting - A Benchmarking¶

In this case study, we explore the capabilities of SysIdentPy by applying it to the Air Passenger dataset, a classic time series dataset widely used for evaluating time series forecasting methods. The primary goal of this analysis is to demonstrate that SysIdentPy can serve as a strong alternative for time series modeling, rather than to assert that one library is superior to another.

Dataset Overview¶

The Air Passenger dataset consists of monthly totals of international airline passengers from 1949 to 1960. This dataset is characterized by its strong seasonal patterns, trend components, and variability, making it an ideal benchmark for evaluating various time series forecasting methods. Specifically, the dataset includes:

• Total Monthly Passengers: The number of passengers (in thousands) for each month.
• Time Period: From January 1949 to December 1960, providing 144 data points.

The dataset exhibits clear seasonal fluctuations and a trend, which poses a significant challenge for forecasting methods. It serves as a well-known benchmark for assessing the performance of different time series models due to its inherent complexity and well-documented behavior.

Comparison with Other Libraries¶

We will compare the performance of SysIdentPy with other popular time series modeling libraries, focusing on the following tools:

• sktime: An extensive library for time series analysis in Python, offering various modeling techniques. For this case study, we will use:
• AutoARIMA: Automatically selects the best ARIMA model based on the data.
• BATS (Bayesian Structural Time Series): A model that captures complex seasonal patterns and trends.
• TBATS (Trigonometric, Box-Cox, ARMA, Trend, and Seasonal): A model designed to handle multiple seasonal patterns.
• Exponential Smoothing: A method that applies weighted averages to forecast future values.
• Prophet: Developed by Facebook, it is particularly effective for capturing seasonality and holiday effects.

• AutoETS (Automatic Exponential Smoothing): Selects the best exponential smoothing model for the data.

• SysIdentPy: A library designed for system identification and time series modeling. We will focus on:

• MetaMSS (Meta-heuristic Model Structure Selection): Uses metaheuristic algorithms to select the best model structure.
• AOLS (Accelerated Orthogonal Least Squares): A method for selecting relevant regressors in a model.
• FROLS (Forward Regression with Orthogonal Least Squares, using polynomial base functions): A regression technique for model structure selection with polynomial terms.
• NARXNN (Nonlinear Auto-Regressive model with Exogenous Inputs using Neural Networks): A flexible method for modeling nonlinear time series with external inputs.

Objective¶

The objective of this case study is to evaluate and compare the performance of these methods on the Air Passenger dataset. We aim to assess how well each library handles the complex seasonal and trend components of the data and to showcase SysIdentPy as a viable option for time series forecasting.

Required Packages and Versions¶

To ensure that you can replicate this case study, it is essential to use specific versions of the required packages. Below is a list of the packages along with their respective versions needed for running the case studies effectively.

To install all the required packages, you can create a requirements.txt file with the following content:

sysidentpy==0.4.0
pystan==2.19.1.1
holidays==0.11.2
fbprophet==0.7.1
neuralprophet==0.2.7
pandas==1.3.2
numpy==1.23.3
matplotlib==3.8.4
pmdarima==1.8.3
scikit-learn==0.24.2
scipy==1.9.1
sktime==0.8.0
statsmodels==0.12.2
tbats==1.1.0
torch==1.12.1

Then, install the packages using:

pip install -r requirements.txt

• Ensure that you use a virtual environment to avoid conflicts between package versions. This practice isolates your project’s dependencies and prevents version conflicts with other projects or system-wide packages. Additionally, be aware that some packages, such as sktime and neuralprophet, may install several dependencies automatically during their installation. Setting up a virtual environment helps manage these dependencies more effectively and keeps your project environment clean and reproducible.
• Versions specified are based on compatibility with the code examples provided. If you are using different versions, some adjustments in the code might be necessary.

Let's begin by importing the necessary packages and setting up the environment for this analysis.

from warnings import simplefilter
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.signal.signaltools

def _centered(arr, newsize):
    # Return the center newsize portion of the array.
    # this is needed due a conflict error using the versions of the packages defined
    # for this example
    newsize = np.asarray(newsize)
    currsize = np.array(arr.shape)
    startind = (currsize - newsize) // 2
    endind = startind + newsize
    myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
    return arr[tuple(myslice)]


scipy.signal.signaltools._centered = _centered
from sysidentpy.model_structure_selection import FROLS
from sysidentpy.model_structure_selection import AOLS
from sysidentpy.model_structure_selection import MetaMSS
from sysidentpy.basis_function import Polynomial
from sysidentpy.utils.plotting import plot_results
from torch import nn
from sysidentpy.neural_network import NARXNN

from sktime.datasets import load_airline
from sktime.forecasting.ets import AutoETS
from sktime.forecasting.arima import ARIMA, AutoARIMA
from sktime.forecasting.base import ForecastingHorizon
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
from sktime.forecasting.fbprophet import Prophet
from sktime.forecasting.tbats import TBATS
from sktime.forecasting.bats import BATS
from sktime.forecasting.model_selection import temporal_train_test_split
from sktime.performance_metrics.forecasting import mean_squared_error
from sktime.utils.plotting import plot_series

from neuralprophet import NeuralProphet
from neuralprophet import set_random_seed

simplefilter("ignore", FutureWarning)
np.seterr(all="ignore")
%matplotlib inline
loss = mean_squared_error

We use the sktime method to load the data. Besides, 23 samples is used as test data, following the definitions in the sktime examples.

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)  # 23 samples for testing
plot_series(y_train, y_test, labels=["y_train", "y_test"])
fh = ForecastingHorizon(y_test.index, is_relative=False)
print(y_train.shape[0], y_test.shape[0])

The following image shows the data of the system to be modeled.

Results¶

Because we have several different models to test, the results are summarized in the following table. The user you will see that no hyperparameter tuning was made for SysIdentPy model. The idea here is to show how simple it can be to build good models in SysIdentPy.

SysIdentPy: FROLS¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
y_train = y_train.values.reshape(-1, 1)
y_test = y_test.values.reshape(-1, 1)
basis_function = Polynomial(degree=1)
sysidentpy = FROLS(
    order_selection=True,
    ylag=13,  # the lags for all models will be 13
    basis_function=basis_function,
    model_type="NAR",
)

sysidentpy.fit(y=y_train)
y_test = np.concatenate([y_train[-sysidentpy.max_lag :], y_test])
yhat = sysidentpy.predict(y=y_test, forecast_horizon=23)
frols_loss = loss(
    pd.Series(y_test.flatten()[sysidentpy.max_lag :]),
    pd.Series(yhat.flatten()[sysidentpy.max_lag :]),
)
print(frols_loss)
plot_results(y=y_test[sysidentpy.max_lag :], yhat=yhat[sysidentpy.max_lag :])
>>> 805.95

SysIdentPy: AOLS¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
y_train = y_train.values.reshape(-1, 1)
y_test = y_test.values.reshape(-1, 1)
df_train, df_test = temporal_train_test_split(y, test_size=23)
df_train = df_train.reset_index()
df_train.columns = ["ds", "y"]
df_train["ds"] = pd.to_datetime(df_train["ds"].astype(str))
df_test = df_test.reset_index()
df_test.columns = ["ds", "y"]
df_test["ds"] = pd.to_datetime(df_test["ds"].astype(str))

sysidentpy_AOLS = AOLS(
    ylag=13, k=2, L=1, model_type="NAR", basis_function=basis_function
)

sysidentpy_AOLS.fit(y=y_train)
y_test = np.concatenate([y_train[-sysidentpy_AOLS.max_lag :], y_test])
yhat = sysidentpy_AOLS.predict(y=y_test, steps_ahead=None, forecast_horizon=23)

aols_loss = loss(
    pd.Series(y_test.flatten()[sysidentpy_AOLS.max_lag :]),
    pd.Series(yhat.flatten()[sysidentpy_AOLS.max_lag :]),
)

print(aols_loss)
plot_results(y=y_test[sysidentpy_AOLS.max_lag :], yhat=yhat[sysidentpy_AOLS.max_lag :])
>>> 476.64

SysIdentPy: MetaMSS¶

set_random_seed(42)
y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
y_train = y_train.values.reshape(-1, 1)
y_test = y_test.values.reshape(-1, 1)

sysidentpy_metamss = MetaMSS(
    basis_function=basis_function, ylag=13, model_type="NAR", test_size=0.17
)

sysidentpy_metamss.fit(y=y_train)

y_test = np.concatenate([y_train[-sysidentpy_metamss.max_lag :], y_test])
yhat = sysidentpy_metamss.predict(y=y_test, steps_ahead=None, forecast_horizon=23)

metamss_loss = loss(
    pd.Series(y_test.flatten()[sysidentpy_metamss.max_lag :]),
    pd.Series(yhat.flatten()[sysidentpy_metamss.max_lag :]),
)

print(metamss_loss)
plot_results(
    y=y_test[sysidentpy_metamss.max_lag :], yhat=yhat[sysidentpy_metamss.max_lag :]
)

>>> 450.99

SysIdentPy: Neural NARX¶

The network architecture is just the same as the one used in to show how to build a Neural NARX model in SysIdentPy docs.

import torch
torch.manual_seed(42)

y = load_airline()
# the split here will use 36 as test size just because the network will use the first values as initial conditions. It could be done like the others methods by concatenating the values
y_train, y_test = temporal_train_test_split(y, test_size=36)
y_train = y_train.values.reshape(-1, 1)
y_test = y_test.values.reshape(-1, 1)
x_train = np.zeros_like(y_train)
x_test = np.zeros_like(y_test)

class NARX(nn.Module):
    def __init__(self):
        super().__init__()
        self.lin = nn.Linear(13, 20)
        self.lin2 = nn.Linear(20, 20)
        self.lin3 = nn.Linear(20, 20)
        self.lin4 = nn.Linear(20, 1)
        self.relu = nn.ReLU()


    def forward(self, xb):
        z = self.lin(xb)
        z = self.relu(z)
        z = self.lin2(z)
        z = self.relu(z)
        z = self.lin3(z)
        z = self.relu(z)
        z = self.lin4(z)
        return z

narx_net = NARXNN(
    net=NARX(),
    ylag=13,
    model_type="NAR",
    basis_function=Polynomial(degree=1),
    epochs=900,
    verbose=False,
    learning_rate=2.5e-02,
    optim_params={},  # optional parameters of the optimizer
)

narx_net.fit(y=y_train)
yhat = narx_net.predict(y=y_test, forecast_horizon=23)

narxnet_loss = loss(
    pd.Series(y_test.flatten()[narx_net.max_lag :]),
    pd.Series(yhat.flatten()[narx_net.max_lag :]),
)

print(narxnet_loss)
plot_results(y=y_test[narx_net.max_lag :], yhat=yhat[narx_net.max_lag :])

sktime models¶

The following models are the ones available in the sktime package.

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)  # 23 samples for testing
plot_series(y_train, y_test, labels=["y_train", "y_test"])
fh = ForecastingHorizon(y_test.index, is_relative=False)

sktime: Exponential Smoothing¶

es = ExponentialSmoothing(trend="add", seasonal="multiplicative", sp=12)
y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
es.fit(y_train)
y_pred_es = es.predict(fh)
plot_series(y_test, y_pred_es, labels=["y_test", "y_pred"])
es_loss = loss(y_test, y_pred_es)
es_loss
>>> 910.46

sktime: AutoETS¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
ets = AutoETS(auto=True, sp=12, n_jobs=-1)
ets.fit(y_train)
y_pred_ets = ets.predict(fh)
plot_series(y_test, y_pred_ets, labels=["y_test", "y_pred"])
ets_loss = loss(y_test, y_pred_ets)
ets_loss
>>> 1739.11

sktime: AutoArima¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)

auto_arima = AutoARIMA(sp=12, suppress_warnings=True)
auto_arima.fit(y_train)
y_pred_auto_arima = auto_arima.predict(fh)
plot_series(y_test, y_pred_auto_arima, labels=["y_test", "y_pred"])
autoarima_loss = loss(y_test, y_pred_auto_arima)
autoarima_loss
>>> 1714.47

sktime: Arima¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
manual_arima = ARIMA(
    order=(13, 1, 0), suppress_warnings=True
)  # seasonal_order=(0, 1, 0, 12)
manual_arima.fit(y_train)
y_pred_manual_arima = manual_arima.predict(fh)
plot_series(y_test, y_pred_manual_arima, labels=["y_test", "y_pred"])
manualarima_loss = loss(y_test, y_pred_manual_arima)
manualarima_loss
>>> 2085.42

sktime: BATS¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
bats = BATS(sp=12, use_trend=True, use_box_cox=False)
bats.fit(y_train)
y_pred_bats = bats.predict(fh)
plot_series(y_test, y_pred_bats, labels=["y_test", "y_pred"])
bats_loss = loss(y_test, y_pred_bats)
bats_loss
>>> 7286.64

sktime: TBATS¶

y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
tbats = TBATS(sp=12, use_trend=True, use_box_cox=False)
tbats.fit(y_train)
y_pred_tbats = tbats.predict(fh)
plot_series(y_test, y_pred_tbats, labels=["y_test", "y_pred"])
tbats_loss = loss(y_test, y_pred_tbats)
tbats_loss
>>> 7448.43

sktime: Prophet¶

set_random_seed(42)
y = load_airline()
y_train, y_test = temporal_train_test_split(y, test_size=23)
z = y.copy()
z = z.to_timestamp(freq="M")
z_train, z_test = temporal_train_test_split(z, test_size=23)
prophet = Prophet(
    seasonality_mode="multiplicative",
    n_changepoints=int(len(y_train) / 12),
    add_country_holidays={"country_name": "Germany"},
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
)
prophet.fit(z_train)
y_pred_prophet = prophet.predict(fh.to_relative(cutoff=y_train.index[-1]))
y_pred_prophet.index = y_test.index
plot_series(y_test, y_pred_prophet, labels=["y_test", "y_pred"])
prophet_loss = loss(y_test, y_pred_prophet)
prophet_loss
>>> 1186.00

Neural Prophet¶

set_random_seed(42)
df = pd.read_csv(r".\datasets\air_passengers.csv")
m = NeuralProphet(seasonality_mode="multiplicative")
df_train = df.iloc[:-23, :].copy()
df_test = df.iloc[-23:, :].copy()

m = NeuralProphet(seasonality_mode="multiplicative")
metrics = m.fit(df_train, freq="MS")
future = m.make_future_dataframe(
    df_train, periods=23, n_historic_predictions=len(df_train)
)
forecast = m.predict(future)
plt.plot(forecast["yhat1"].values[-23:])
plt.plot(df_test["y"].values)
neuralprophet_loss = loss(forecast["yhat1"].values[-23:], df_test["y"].values)
neuralprophet_loss
>>> 501.24

The final results can be summarized as follows, resulting in the table presented in the beginning of this case study:

results = {
    "Exponential Smoothing": es_loss,
    "ETS": ets_loss,
    "AutoArima": autoarima_loss,
    "Manual Arima": manualarima_loss,
    "BATS": bats_loss,
    "TBATS": tbats_loss,
    "Prophet": prophet_loss,
    "SysIdentPy (Polynomial Model)": frols_loss,
    "SysIdentPy (Neural Model)": narxnet_loss,
    "SysIdentPy (AOLS)": aols_loss,
    "SysIdentPy (MetaMSS)": metamss_loss,
    "NeuralProphet": neuralprophet_loss,
}
sorted(results.items(), key=lambda result: result[1])

System With Hysteresis - Modeling a Magneto-rheological Damper Device¶

The memory effects between quasi-static input and output make the modeling of hysteretic systems very difficult. Physics-based models are often used to describe the hysteresis loops, but these models usually lack the simplicity and efficiency required in practical applications involving system characterization, identification, and control. As detailed in Martins, S. A. M. and Aguirre, L. A., NARX models have proven to be a feasible choice to describe the hysteresis loops. See Chapter 8 for a detailed background. However, even considering the sufficient conditions for rate independent hysteresis representation, classical structure selection algorithms fails to return a model with decent performance and the user needs to set a multi-valued function to ensure the occurrence of the bounding structure \(\mathcal{H}\) (Martins, S. A. M. and Aguirre, L. A.).

Even though some progress has been made, previous work has been limited to models with a single equilibrium point. The present case study aims to present new prospects in the model structure selection of hysteretic systems regarding the cases where the models have multiple inputs and it is not restricted concerning the number of equilibrium points. For that, the MetaMSS algorithm will be used to build a model for a magneto-rheological damper (MRD) considering the mentioned sufficient conditions.

A Brief description of the Bouc-Wen model of magneto-rheological damper device¶

The data used in this study-case is the Bouc-Wen model (Bouc, R, Wen, Y. X.) of an MRD whose schematic diagram is shown in the figure below.

The model for a magneto-rheological damper proposed by Spencer, B. F. and Sain, M. K..

The general form of the Bouc-Wen model can be described as (Spencer, B. F. and Sain, M. K.):

where \(z\) is the hysteretic model output, \(x\) the input and \(g[\cdot]\) a nonlinear function of \(x\), \(z\) and \(sign (dx/dt)\). Spencer, B. F. and Sain, M. K. proposed the following phenomenological model for the aforementioned device:

where \(f\) is the damping force, \(c_1\) and \(c_0\) represent the viscous coefficients, \(E\) is the input voltage, \(x\) is the displacement and \(\dot{x}\) is the velocity of the model. The parameters of the system (see table below) were taken from Leva, A. and Piroddi, L..

For this particular study, both displacement and voltage inputs, \(x\) and \(E\), respectively, were generated by filtering a white Gaussian noise sequence using a Blackman-Harris FIR filter with \(6\)Hz cutoff frequency. The integration step-size was set to \(h = 0.002\), following the procedures described in Martins, S. A. M. and Aguirre, L. A.. These procedures are for identification purposes only since the inputs of a MRD could have several different characteristics.

The data used in this example is provided by the Professor Samir Angelo Milani Martins.

The challenges are:

• it possesses a nonlinearity featuring memory, i.e. a dynamic nonlinearity;
• the nonlinearity is governed by an internal variable z(t), which is not measurable;
• the nonlinear functional form in the Bouc Wen equation is nonlinear in the parameter;
• the nonlinear functional form in the Bouc Wen equation does not admit a finite Taylor series expansion because of the presence of absolute values

Required Packages and Versions¶

To ensure that you can replicate this case study, it is essential to use specific versions of the required packages. Below is a list of the packages along with their respective versions needed for running the case studies effectively.

To install all the required packages, you can create a requirements.txt file with the following content:

sysidentpy==0.4.0
pandas==2.2.2
numpy==1.26.0
matplotlib==3.8.4
scikit-learn==1.4.2

Then, install the packages using:

pip install -r requirements.txt

• Ensure that you use a virtual environment to avoid conflicts between package versions.
• Versions specified are based on compatibility with the code examples provided. If you are using different versions, some adjustments in the code might be necessary.

SysIdentPy Configuration¶

import numpy as np
from sklearn.preprocessing import MaxAbsScaler, MinMaxScaler
import pandas as pd
import matplotlib.pyplot as plt

from sysidentpy.model_structure_selection import FROLS
from sysidentpy.basis_function import Polynomial
from sysidentpy.utils.display_results import results
from sysidentpy.parameter_estimation import LeastSquares
from sysidentpy.metrics import root_relative_squared_error
from sysidentpy.utils.plotting import plot_results

df = pd.read_csv("boucwen_histeretic_system.csv")
scaler_x = MaxAbsScaler()
scaler_y = MaxAbsScaler()

init = 400
x_train = df[["E", "v"]].iloc[init:df.shape[0]//2, :]
x_train["sign_v"] = np.sign(df["v"])
x_train = scaler_x.fit_transform(x_train)

x_test = df[["E", "v"]].iloc[df.shape[0]//2 + 1:df.shape[0] - init, :]
x_test["sign_v"] = np.sign(df["v"])
x_test = scaler_x.transform(x_test)

y_train = df[["f"]].iloc[init:df.shape[0]//2, :].values.reshape(-1, 1)
y_train = scaler_y.fit_transform(y_train)

y_test = df[["f"]].iloc[df.shape[0]//2 + 1:df.shape[0] - init, :].values.reshape(-1, 1)
y_test = scaler_y.transform(y_test)

# Plotting the data
plt.figure(figsize=(10, 8))
plt.suptitle('Identification (training) data', fontsize=16)

plt.subplot(221)
plt.plot(y_train, 'k')
plt.ylabel('Force - Output')
plt.xlabel('Samples')
plt.title('y')
plt.grid()
plt.axis([0, 1500, -1.5, 1.5])

plt.subplot(222)
plt.plot(x_train[:, 0], 'k')
plt.ylabel('Control Voltage')
plt.xlabel('Samples')
plt.title('x_1')
plt.grid()
plt.axis([0, 1500, 0, 1])

plt.subplot(223)
plt.plot(x_train[:, 1], 'k')
plt.ylabel('Velocity')
plt.xlabel('Samples')
plt.title('x_2')
plt.grid()
plt.axis([0, 1500, -1.5, 1.5])

plt.subplot(224)
plt.plot(x_train[:, 2], 'k')
plt.ylabel('sign(Velocity)')
plt.xlabel('Samples')
plt.title('x_3')
plt.grid()
plt.axis([0, 1500, -1.5, 1.5])

plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()

Let's check how is the histeretic behavior considering each input:

plt.plot(x_train[:, 0], y_train)
plt.xlabel("x1 - Voltage")
plt.ylabel("y - Force")

plt.plot(x_train[:, 1], y_train)
plt.xlabel("x2 - Velocity")
plt.ylabel("y - Force")

plt.plot(x_train[:, 2], y_train)
plt.xlabel("u3 - sign(Velocity)")
plt.ylabel("y - Force")

Now, we can just build a NARX model:

basis_function = Polynomial(degree=3)
model = FROLS(
    xlag=[[1], [1], [1]],
    ylag=1,
    basis_function=basis_function,
    estimator=LeastSquares(),
    info_criteria="aic",
)

model.fit(X=x_train, y=y_train)
yhat = model.predict(X=x_test, y=y_test[:model.max_lag :, :])
rrse = root_relative_squared_error(y_test[model.max_lag :], yhat[model.max_lag :])
print(rrse)
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title="FROLS: sign(v) and MaxAbsScaler")
>>> 0.0450

If we remove the sign(v) input and try to build a NARX model using the same configuration, the model diverge, as can be seen in the following figure:

basis_function = Polynomial(degree=3)
model = FROLS(
    xlag=[[1], [1]],
    ylag=1,
    basis_function=basis_function,
    estimator=LeastSquares(),
    info_criteria="aic",
)

model.fit(X=x_train[:, :2], y=y_train)
yhat = model.predict(X=x_test[:, :2], y=y_test[:model.max_lag :, :])
rrse = root_relative_squared_error(y_test[model.max_lag :], yhat[model.max_lag :])
print(rrse)
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title="FROLS: MaxAbsScaler, discarding sign(v)")
>>> nan

If we use the MetaMSS algorithm instead, the results are better.

from sysidentpy.model_structure_selection import MetaMSS

basis_function = Polynomial(degree=3)
model = MetaMSS(
    xlag=[[1], [1]],
    ylag=1,
    basis_function=basis_function,
    estimator=LeastSquares(),
    random_state=42,
)

model.fit(X=x_train[:, :2], y=y_train)
yhat = model.predict(X=x_test[:, :2], y=y_test[:model.max_lag :, :])
rrse = root_relative_squared_error(y_test[model.max_lag :], yhat[model.max_lag :])
print(rrse)
plot_results(y=y_test[model.max_lag :], yhat=yhat[model.max_lag :], n=10000, title="MetaMSS: MaxAbsScaler, discarding sign(v)")
>>> 0.24