2.1.1. Why do we benchmark?

A first advantage of benchmarking is the performance assessment of a developed approach. Furthermore, if the chosen benchmark datasets cover relevant challenges in the field, obtaining good results on these benchmarks establishes the credibility of the proposed approach. Finally, when using widespread benchmarks, the competitiveness of the proposed approach can be established and validated at a low cost for the authors, as many other results on the same data will be available in the literature. This bypasses the time-consuming task of implementing and running other approaches on the examples one wants to compare. Finally, a well-constructed benchmark can also raise the attention of open challenges in the research field.

2.1.2. What constitutes a good benchmark?

A good benchmark problem consists of many different components. A clear overview of the main elements involved is provided by Reddi (Reference Reddi2024, Chapter 11): standardized datasets, predefined tasks, evaluation metrics, baseline results, reproducibility, and results interpretation guidelines. In the next paragraphs, we revisit these aspects and briefly discuss them from a system identification point of view.

2.1.3. Standardized datasets

Data forms the starting point of most identification algorithms. To be able to compare multiple identification methods, the same data should be used. Preferably, these datasets should match the wide variety of data encountered in engineering practice ranging from low to high signal-to-noise ratios, from short (expensive) to long (cheap) experiments, and from highly informative to less informative experiments. However, starting from a fixed dataset ignores the important input/experiment design aspect of the system identification cycle. Hence, we should also accommodate benchmark systems where custom datasets can be generated.

2.1.4. Predefined tasks and evaluation metrics

The end goal of a model is often not the model itself, but, for example, the design of a controller, a structural health monitoring, or a predictive maintenance tool. Nevertheless, if we want to provide a comparison of identification algorithms, a clearly defined task needs to be specified. In system identification tasks, this is often minimizing the (Root) Mean Square Error (RMSE), or a normalized variation of it (e.g. the best-fit rate), between the actual measured outputs and the predicted or simulated output. In some cases, some distance between the actual system parameters (e.g., resistor and capacitor values in an electric network) and the estimated model parameters are considered, but this significantly restricts the model structures that can be used as it should be closely tied to the considered nature of the parameters.

At the same time, the RMSE is often used as the main quantitative evaluation metric. However, focusing on such a metric alone would be damaging to the general development of the field as it ignores a wide array of qualities that can be incorporated by a method or by the identified model and the algorithm used to estimate it: e.g., the model and algorithmic complexity, the robustness of the algorithm, uncertainty quantification, or the computational demand for the estimation and the use of the model. Clearly, one needs to leave room for nuance and detailed explanation of the reported results to account for a wide range of user demands due to various end uses of an estimated model. Examples of such considerations include, for instance, the recent awareness regarding “green AI” and “green-box” system identification: e.g., a model that has a 10% higher error but that can be obtained with only 1% of the computational cost can be highly relevant for a wide range of applications (Batselier, Reference Batselier2022; Memmel et al., Reference Memmel, Menzen, Schuurmans, Wesel and Batselier2024).

2.1.5. Result interpretation

As argued earlier, comparing identification approaches on a limited number of figures of merit only (e.g., RMSE), would not be effective. The end use of the model often adds additional constraints. Hence, when interpreting results, the qualitative description of the benefits and limitations of an approach is paramount. Furthermore, benchmark problems only cover a limited set of the full problem space. Focusing only on benchmark results would bias the field too much toward a (limited) set of (nonlinear) system identification problems.

2.1.6. Baseline results

Including baseline results and their implementations provides a key idea of what minimum performance should be expected, and what the classical tools within the field can obtain. An obvious baseline result for nonlinear system identification could, for instance, be the result of a linear model. Any nonlinear model performing worse than this is most likely of little interest to the community. At the same time, these baselines limit the danger of cherry-picking the approaches one compares with when they aim to publish their results. As the baselines are published together with the benchmark data or are very visible, reviewers will also have easy access to them.

3.1.1. Wiener–Hammerstein (WH)

The Wiener–Hammerstein system class is a well-known block-oriented structure (Schoukens and Tiels, Reference Schoukens and Tiels2017). It contains a static nonlinearity that is sandwiched in between two LTI blocks. As the signal that governs the nonlinear behavior is not directly measured (it is hidden by the two LTI blocks), this results in a problem that is harder to identify, e.g., input–output representations such as NAR(MA)X models.

The Wiener–Hammerstein benchmark (Schoukens and Ljung, Reference Schoukens and Ljung2009) was the subject of an invited session during the IFAC Symposium on System Identification in St. Malo, France in 2009, and has attracted quite some attention since. The system is an electronic circuit with a Wiener–Hammerstein structure built by Gerd Vandersteen (Reference Vandersteen1997), where the nonlinearity is realized as a diode resistor network. It is excited by a highly informative input signal and features a low-noise, large 188,000 samples long measured dataset.

3.1.2. Parallel Wiener–Hammerstein (ParWH)

A Parallel Wiener–Hammerstein system is obtained by connecting multiple Wiener–Hammerstein systems in parallel. The system considered in this benchmark is an electronic circuit with a two-branch parallel Wiener–Hammerstein structure (Schoukens et al., Reference Schoukens, Marconato, Pintelon, Vandersteen and Rolain2015). The front and back LTI blocks of each branch are third order continuous-time IIR filters. The static nonlinearity of each branch is realized with a diode-resistor network. The measured training dataset is over 3 million data points large, containing measurements at 5 different amplitude levels, with multiple realizations of a random phase multisine input signal at each amplitude level (Schoukens et al., Reference Schoukens, Marconato, Pintelon, Vandersteen and Rolain2015; Schoukens, Reference Schoukens2020). Twenty realizations at each amplitude level can be used for training and validation, and one realization per amplitude level is available for testing. The aforementioned challenges in identifying a Wiener–Hammerstein model are further complicated due to its parallel branch nature and the increased order of the system dynamics (12th order).

3.1.3. Silverbox (SB)

The Silverbox system is an electronic implementation of a mass-spring-damper system with a cubic spring (Wigren and Schoukens, Reference Wigren and Schoukens2013). This is also sometimes known as a forced Duffing oscillator. However, the current benchmark is constructed such that it never enters a bifurcating regime: for the applied input signals, a fading memory behavior is observed. Similar to the Wiener–Hammerstein case, the Silverbox system has been widely studied in the literature, and this benchmark also offers highly informative data, in a low-noise setting, and a long data record. It can be considered as one of the “entry-level” benchmark problems that are considered here. The training dataset contains concatenated multisine realizations, while the testing dataset contains a multisine excitation and a random Gaussian noise input with growing amplitude experiment.

3.1.4. Cascaded Tanks (CTs)

The cascaded tanks system consists of two water tanks with free outlets fed by a pump (Schoukens et al., Reference Schoukens, Mattsson, Wigren and Noël2016). The input controls a water pump that delivers the water from a reservoir into the upper water tank. The water of the upper tank flows through a small opening into the lower tank, and finally through a small opening from the lower tank back into the reservoir. The water level of the lower tank is considered as the system output (measured in V using a capacitive sensor). This benchmark combines soft and hard (due to overflow) nonlinearities to be identified based on relatively short data records.

3.1.5. Electro-Mechanical Positioning System (EMPS)

The EMPS is a standard configuration of a drive system for prismatic joints of robots or machine tools (Janot et al., Reference Janot, Gautier and Brunot2019). The main source of nonlinearity is the friction effects that are present in the setup. Due to the presence of a pure integrator in the system, the measurements are obtained in a closed-loop setting.

3.1.6. Coupled Electric Drives (CEDs)

The CED system consists of two electric motors that drive a pulley using a flexible belt (Wigren and Schoukens, Reference Wigren and Schoukens2017). The pulley is held by a spring, resulting in a (lightly) damped dynamic mode. The drive control is symmetric around zero, hence both clockwise and counter-clockwise movements are possible. The angular speed of the pulley is measured as an output with a pulse counter, and this sensor is insensitive to the sign of the velocity. The available datasets are short, which presents a challenge when identifying a model.

3.1.7. F16 Ground Vibration Test (F16)

The F-16 Ground Vibration Test benchmark (Noël and Schoukens, Reference Noël and Schoukens2020) features a high-order system with clearance and friction nonlinearities at the mounting interface of the payloads. The experimental data made available were acquired on a full-scale F-16 aircraft on the occasion of the Siemens LMS Ground Vibration Testing Master Class, held in September 2014. During the test campaign, two dummy payloads were mounted at the wing tips to simulate the mass and inertia properties of real devices typically equipping an F-16 in flight. One shaker was attached underneath the right wing to apply input signals. The aircraft structure was instrumented with accelerometers. A preliminary investigation showed that the back connection of the right-wing-to-payload interface was the predominant source of nonlinear distortions in the aircraft dynamics, and is therefore the focus of this benchmark study. Measurements were performed using multisine and sinesweep experiments at 7 different amplitude levels each. Amplitude levels 1, 3, 5, and 7 are used for training and validation while the others are used for testing. Note that the “standardized dataset” used here only comprises the random phase multisine and sinesweep measurements; the odd multisine measurements are omitted as they cover a much broader frequency range. The standardized dataset also only considers the first of the three offered outputs in the dataset.

4.1.1. LTI SS

Discrete-time LTI state-space models have been commonly used for decades in the system identification literature. Here, LTI state-space models with an output error noise model are estimated:

(4.1) $$ {\displaystyle \begin{array}{c}{x}_{t+1}={Ax}_t+{Bu}_t\\ {}{\hat{y}}_t={Cx}_t+{Du}_t\end{array}} $$

where $ {u}_t $ is the input, $ {x}_t $ is the state, and $ {\hat{y}}_t $ is the modeled output at time $ t $ . The matrices $ A $ , $ B $ , $ C $ , and $ D $ are paramerized in a “free” form, that is, any element in the matrix is freely adjustable. The state-space models are estimated using the prediction error framework (ssest command in Matlab). The resulting nonlinear optimization problem is initialized using a subspace estimate (Ljung, Reference Ljung1999).

The means of the measured input and output training data are subtracted before estimating the LTI state-space models. Hence, minimizing the prediction error results in a Best Linear Approximation model of the systems considered, following the lines of Enqvist and Ljung (Reference Enqvist and Ljung2005) and Pintelon and Schoukens (Reference Pintelon and Schoukens2012). The resulting mean values of the input and output training data are part of the resulting model as they are respectively subtracted from and added to any validation or test input provided to the model and resulting model output.

The number of states considered by the LTI SS approach for each benchmark is reported in Table 1. The state order $ {n}_x $ is selected based on the provided benchmark system descriptions, except for the EMPS benchmark, where a grid search with validation dataset illustrated that $ {n}_x=4 $ resulted in significantly better results than $ {n}_x=2 $ or $ {n}_x=5 $ , and slightly better than $ {n}_x=3 $ . The F16 LTI state-space model is estimated using data that was 8x downsampled (using the Matlab downsample command). To report the results, the model data is upsampled back to the original sampling rate (using the Matlab interp command).

Table 1.

Number of states used by the LTI SS models considered in the present paper

Table 1. Long description

The table consists of two rows and eight columns.

Row 1 (Headers): The first cell is empty. The subsequent columns are labeled S B, E M P S, W H, C T, C E D, Par W H, and F 16.

Row 2 (Data): The first cell contains the variable n sub x. The corresponding values for the models are:

* S B: 2

* E M P S: 4

* W H: 6

* C T: 2

* C E D: 3

* Par W H: 12

* F 16: 20

4.1.2. LTI ARX

This model is described in more detail in the auto-regressive models subsection.

4.2.1. LTI ARX

The first auto-regressive model considered here is the LTI auto-regressive with exogenous inputs (ARX) model. This is a linear formulation with,

(4.4) $$ \Phi (H)= H\alpha $$

where $ \alpha $ is an unknown vector of linear weights that can be evaluated by the solution of the linear equation $ \mathbf{y}= H\alpha $ . Many solvers are appropriate for this task. The results in the present paper use a QR-decomposition-based approach.

4.2.2. pNARX

The first nonlinear ARX (NARX) model considered in this study is the classical polynomial NARX model (Billings, Reference Billings2013). Here, the nonlinearity is defined as a linear combination of polynomial terms $ {P}_n $ computed from $ H $ :

(4.5) $$ \Phi (H)={P}_n(H)\hskip0.1em \alpha $$

To promote stability, only univariate polynomial terms are retained, that is, each basis function depends on a single regressor only, and cross-terms such as $ {y}_{t-1}^2{u}_{t-1} $ are excluded. A Legendre orthogonal polynomial basis up to order $ n $ is applied to each regressor independently to improve numerical conditioning. Concretely, for a single row $ {h}^{\top }=\left[{y}_{t-1},\dots, {y}_{t-{n}_y},{u}_t,\dots, {u}_{t-{n}_u+1}\right]\in {\mathrm{\mathbb{R}}}^r $ of $ H $ (with $ r={n}_y+{n}_u $ ), the feature map is

(4.6) $$ \phi \left({h}^{\top}\right)=\left[{\mathrm{\ell}}_1\left({y}_{t-1}\right),\dots, {\mathrm{\ell}}_n\left({y}_{t-1}\right),\dots, {\mathrm{\ell}}_1\left({u}_{t-{n}_u+1}\right),\dots, {\mathrm{\ell}}_n\left({u}_{t-{n}_u+1}\right)\right]\in {\mathrm{\mathbb{R}}}^{1\times rn} $$

Where $ {\mathrm{\ell}}_k $ is the $ k $ th term in the Legendre basis, so that $ \Phi (H)\in {\mathrm{\mathbb{R}}}^{N\times rn} $ , and $ \alpha \in {\mathrm{\mathbb{R}}}^{rn} $ . The order of the polynomial is selected based on AIC cross-validation scores on the hold-out validation set, up to a maximum order of 10. The identified model orders are reproduced in Table 3 for the convenience of the reader. It is expected that the inclusion of multinomial terms may improve the identification on some benchmarks. However, a full multinomial basis leads to a very high number of basis terms for only a modest number of lags and so is not considered here.

Table 3.

Polynomial orders used by the PNARX models considered in the present paper

Table 3. Long description

The table consists of two rows and eight columns.

Row 1 (Header): The first cell is empty. The subsequent headers from left to right are S B, E M P S, W H, C T, C E D, Par W H, and F 16.

Row 2 (Polynomial order):

- S B: 7

- E M P S: 3

- W H: 5

- C T: 3

- C E D: 2

- Par W H: dash (no value)

- F 16: 2

4.2.3. GP NARX

Another common choice for the nonlinearity in NARX models is to use a Gaussian process (GP) model. GPs offer flexible kernel regression and quantification of uncertainty with few hyperparameters (Rasmussen and Williams, Reference Rasmussen and Williams2006). The nonlinearity is therefore,

(4.7) $$ \Phi (H)\sim \mathcal{GP}\left(m(K),K\Big(H,H\Big)\right) $$

The GP is entirely parametrized by the kernel function $ K $ and the mean function $ m $ . In this study, GP NARX models are constructed with a squared-exponential kernel and zero-mean function. Training GPs requires the specification of several kernel hyperparameters $ {\theta}_K=\{{\sigma}_n^2,{\sigma}_f^2,\ell \} $ . By convention, these are optimized by minimization of the negative-log marginal-likelihood with respect to the training data (Rasmussen and Williams, Reference Rasmussen and Williams2006). For a zero-mean GP, this is given by,

(4.8) $$ -\log p\left({\theta}_K|\mathbf{y}\right)\propto \frac{1}{2}{\mathbf{y}}^T{\Sigma}^{-1}\mathbf{y}+\frac{1}{2}\log \mid \Sigma \mid $$

with $ \Sigma =K\left(H,H\right)+{\sigma}_n^2I $ , for the proposed noise variance term $ {\sigma}_n^2 $ .

Here, this is achieved with the well-known stochastic gradient descent algorithm ADAM (Kingma and Ba, Reference Kingma and Ba2014). For each benchmark, the hyperparameters are optimized from random initializations for a total of 1000 steps with a learning rate of $ {10}^{-2} $ . The hyperparameters with the best score on the hold-out validation set are then used for evaluation of the test set.

Due to the computational cost ( $ \mathcal{O}\left({n}^3\right) $ for $ n $ datapoints) of training GPs, for benchmark sets with more than 1000 examples (EMPS, WH, SB, ParWH, F16), only the first 1000 points of the training set are used to train the model. For training on larger sequences, sparse GP approaches (such as those in Titsias, Reference Titsias2009; Hensman et al., Reference Hensman, Fusi and Lawrence2013) could be applied; however, these are not considered in this work following our philosophy of producing only “baseline” results to which more complex approaches can be compared.

4.2.4. MLP NARX

The final nonlinearity considered for NARX models in this study is to use a multilayer perceptron (MLP) for $ \Phi $ . Here, an MLP with a single hidden layer and hyperbolic tangent activation functions is considered. The number of nodes in the hidden layer is selected by cross-validation in the hold-out validation set with sizes of 2, 4, 8, 16, and 32 considered. Lags are selected by cross-validation of an ARX model by the procedure described earlier. Each network is trained for $ {10}^4 $ iterations of the ADAM (Kingma and Ba, Reference Kingma and Ba2014) optimizer, with 5 random initializations except for the two smallest benchmark datasets (CED and CT), which are optimized for $ 2\times {10}^4 $ iterations. The overall best network is selected as the one with the best score on the hold-out validation set. This network is then used for evaluation on the unseen testing set.

4.3.1. MLP FIR

The first structure considered in this study is a (nonlinear) finite impulse response (FIR) model. This model structure is equivalent to removing the hidden ‘memory” state from the recurrence relation given earlier.

(4.10) $$ {\mathbf{y}}_{t+1}=F\left({u}_{t+1},\dots {u}_{t-{n}_u+1}\right) $$

Where $ F $ is an MLP with a single hidden layer. The MLP FIR is essentially equivalent to the MLP NARX model with no lags on the output. Although this architecture is not recurrent, it will serve to highlight the extent to which a recurrent formulation improves the identification for each of the benchmark systems.

4.3.2. RNN

The canonical formulation of an RNN (Schuster and Paliwal, Reference Schuster and Paliwal1997) constructs the hidden cell as a single dense layer in an MLP, with a hyperbolic tangent activation function – other choices are clearly available. This architecture is employed here with the size of the hidden layer decided by cross-validation as described earlier.

4.3.3. GRU

A known issue with the RNN is that as the gradient is propagated backward in time, successive applications of the chain rule over activation functions can cause the numerical values of the gradients to become arbitrarily small or large. This issue is known as vanishing (resp. exploding) gradients (Pascanu et al., Reference Pascanu, Mikolov and Bengio2013). Several variants of the basic RNN cell have been proposed to try and combat this issue, by applying a “gating” function (often a sigmoid or similar) to the value of the memory state in each iteration. The gated recurrent network (GRU) (Chung et al., Reference Chung, Gulcehre, Cho and Bengio2014) is one such architecture that employs this technique.

4.3.4. LSTM

The original extension to the RNN is the famed long-short term memory (LSTM) unit (Hochreiter and Schmidhuber, Reference Hochreiter and Schmidhuber1997). This architecture propagates an extended hidden state consisting of both “long-term” and “short-term” memory. The optimised LSTM (OLSTM) is a variant of the LSTM cell that concatenates the hidden states prior to the activation functions for efficient computation. Results from both the LSTM and OLSTM networks are reported in this study.