---
abstract: 'Hysteresis-controlled devices are widely used in industrial applications. For example, cooling devices usually contain a two-point controller, resulting in a nonlinear hybrid system with two discrete states. Dynamic models of systems are essential for optimizing such industrial supply technology. However, conventional system identification approaches cannot handle hysteresis-controlled devices. Thus, the new identification method Sparse Identification of Nonlinear Dynamics (SINDy) is extended to handle hybrid systems. In this new method (SINDyHybrid), tailored basis functions in form of relay hysterons are added to the library which is used by SINDy. Experiments with a hysteresis controlled water basin show that this approach correctly identifies state transitions of hybrid systems and also succeeds in modeling the dynamics of the discrete system states. A novel proximity hysteron achieves the robustness of this method. The impacts of the sampling rate and the signal noise ration of the measurement data are examined accordingly.'
author:
- 
- 
bibliography:
- 'IEEEabrv.bib'
- 'literature.bib'
title: 'System identification of a hysteresis-controlled chiller plant using SINDy [^1] '
---

SINDyHybrid, system identification, nonlinear systems, hysteresis, hybrid dynamical systems, sparse regression

Introduction
============

Nowadays, data-driven approaches for modeling dynamics without a fixed model structure gain relevance. However, classical grey-box modeling is still popular because it allows to include domain-specific knowledge, constraints, and to achieve transparency. The presented work provides an example of combining explicit knowledge with a flexible data-driven approach in order to estimate hysteresis behavior.

A traditional approach of obtaining a system model combines physical laws and technical parameters to so called white-box-models. Alternatively, grey-box-models can be used, which have a fixed model structure but chosen parameters are tuned using measurement data. Methods from machine learning partially abandon fixed model structures, like in the case of neural networks. These approaches do not provide a clear transparency, which implies a need for empirical validation. Furthermore, pure data-driven methods require for extensive data for all considered situations and operation points. However, due to high costs, efforts or dangers it is not always feasible to conduct real experiments.

This article proposes an extension of the Sparse Identification of Nonlinear Dynamics (SINDy) framework developed and presented by S. Brunton and J.N. Kutz in their prominent publication [@BruntonSINDy.2016]. Adding the robust identification of hysteresis behavior in hybrid systems holds the potential to identify all discrete states of a hysteresis controlled system and build a nonlinear model for the dynamics on all states in one step. The approach aims to handle diverse systems, to include system knowledge and to result in an easily interpretable model.

Section \[sec:stateOfTheArt\] will cover basic theory regarding system identification and current approaches. A new concept is presented in Section \[sec:Concept\]. Part \[sec:practicalEvaluation\] contains the practical experiments and their results while Part \[sec:discussion\] gives a deepening discussion. Section \[sec:conclusion\] summarizes the findings and gives a brief outlook to further applications.

State of the art {#sec:stateOfTheArt}
================

Whenever a certain degree of precision is required for tasks like optimization or anomaly detection, simple dynamic models may deviate too much from the actual system, especially in the critical regions of hysteretic systems, where the state transitions happen. In order to pursue model-based approaches for these tasks, sophisticated models are required.

System Identification for nonlinear hybrid systems
--------------------------------------------------

### Hybrid systems.

The considered type of dynamic systems contains both discrete and continuous parts. Depending on a discrete state, a different behavior is observed. A famous example is the temperature characteristics of water: the behavior changes with reaching the threshold of 100 $^\circ C$ due to the phase transition. Therefore, it makes sense to create two distinct models: one for the discrete state “temperature below 100 $^\circ C$" and one for “temperature above 100 $^\circ C$".

### PWARX models. {#PWARX}

Piece-Wise-Affine (PWA) methods divide the state-space into discrete regions, commonly shaped as polyhedronscite[Ferrari-Trecate.2000]{}. For every region a continuous model is set up, e.g. the well known linear *autoregressive model with exogenous input* (ARX-model). This procedure is called PWARX[@pwarx_paoletti]. If the transition from one state to another is affected by hysteresis behavior, PWA models are not sufficient anymore, because they assume fixed transition borders from one discrete state to another. They are unable to include memory of the past trajectory of the system’s state. [@Fang.2015] present a data-driven approach to model sticky valves using the Preisach model for hysteresis, which originated from describing natural hysteresis in ferromagnetic materials. They succeed in modeling the state transitions, but they separate the identification of the transition model from the identification of dynamic models for the discrete states. Also, the number of discrete states must be known beforehand. This will be overcome by using SINDyHybrid.

### SARX models. {#SARX}

**S**witched **A**ffine Modelle (*SARX*) models can be used if transitions of the state are triggered by external events. It is also possible to use **p**iece**w**ise continuously differentiable **n**onlinear functions (*PWNARX*, *SNARX*). An overview of related techniques is given in[@Hutchison.2008]. Furthermore, the authors present a support-vector regression to separate the domains with hyperplanes of the form $$\label{eq:pwa-svr}
    h(x) = \sum_{i=1}^{N}(\beta_i k_c (x,x_i) + b_c) = 0 ,$$ with $N$ samples of p-dimensional data $x = [x_1, ..., x_N]$, kernel function $k_c(.,.)$ and tunable parameters $\beta_i, b_c$. The sign of $h$ gives the classification.

### Black-box-techniques.

An alternative approach would be to use neural networks to learn the system dynamics, but the resulting models are usually difficult to check for plausibility and stability, which poses further problems for industrial applications[@Ljung_BlackBox_2001].

Hysteresis Modeling {#hysteresisModeling}
-------------------

Whenever a system’s behavior is not only depending on its current internal states, but also on its past trajectory, the system is called hysteretic. This nonlinear effect occurs in nature phenomena like in magnetization of ferromagnetic materials, but also in artificial systems like two-point controllers. The Hungarian physicist Ferenc Preisach presented the first suitable model to describe hysteretic behavior of magnetism in 1935[@Preisach.1935; @PreisachInMemorial]. His hysteron-operator allows for approximation of natural hysteresis using weighted stairs. The relay operator $R_{\alpha, \beta}$, also called relay hysteron, is expressed by Eq. and illustrated in Fig.\[fig:preisachRelay\]. $$\label{eq:relayHysteron}
    y(x) = \begin{cases}
    1, & \text{if $x \geq \beta$.} \\
    0, & \text{if $x \leq \alpha$.} \\
    y_p, & \text{if $ \alpha < x < \beta $.}
    \end{cases}$$

![Relay hysteron [@RelayHyst1][]{data-label="fig:preisachRelay"}](Preisach_RelaySVG){width="0.82\linewidth"}

The value $y_p$ always corresponds to the prior value of $y(t)$. The complex shape of a magnetic hysteresis curve can be achieved by superposition of many elementary relay hysterons.

Modern approaches of hysteresis modeling are discussed in[@Hassani.2014]. The Krasnosel’skii–Pokrovskii (KP) model[@Krasnoselskii.1989] formalizes Preisach’s model in a mathematical way and allows for different operators to be included. Next to operator-based approaches, hysteresis can also be described by differential equations, e.g. using the Bouc–Wen model[@wen1976method] or the extended Bouc-Wen model, which requires 13 parameters to approximate the complex shape of the hysteresis[@Zhang.2001]. The here presented work utilizes the concept of hysterons which is the basis of Preisach’s model.

System Identification
---------------------

The SINDy framework is based upon the concepts of the Koopman operator and of sparsity. This section will give a brief introduction to these concepts, as well as present an overview of the SINDy methodology and its applications. In the following Section \[sec:Concept\], a concept is shown which implements the hysteresis model into the given framework.

### Koopman operator.

The Koopman operator $\mathcal{K}$[@Koopman.1931] is a linear but typically infinite-dimensional operator, which is capable of describing the full dynamics of an underlying system. For this paper we focus on the practical implementation of Koopman operator theory. An in-depth overview about theoretical concepts and applications can be found in [@Budisic.2012].

Each dynamic system of the form with state $x \in \mathbb{R}^n$ and time $t$ has a linear representation of the form $$\frac{d}{dt}g(x) = \mathcal{K} g(x),$$

where $g = [h_1(x), ...]^T$ is a function of possibly infinite observables $h_i(x)$ of the state. The Koopman operator $\mathcal{K}$ evolves observable functions $h_i(x)$ instead of direct traces of the states, which can be arbitrary functions. An example is shown in Eq. . $$\label{eq:observables_example}
g(x) = [h_1(x), h_2(x)...]^T = [1, x, x^2, x^3, \sin(x), \sqrt{x}, e^x, ...]^T$$

These observable functions span a subspace of the Hilbert space[@BruntonKoopmanInvariant.2016b]. In general, the resulting system is infinite dimensional but linear. This enables tools like spectral analysis without any inaccuracy caused by linearization[@Budisic.2012]. In 2005, Igor Mezic presented the application of the Koopman operator for spectral analysis of high-dimensional, nonlinear systems[@Mezic.2005]. The Dynamic Mode Decomposition (DMD) extracts eigenvalues and modes of the Koopman operator from measurement data[@Schmid2010], i.e. it finds the most probable system matrix $A$ for The extended-DMD (e-DMD) additionally approximates eigenfunctions of the Koopman operator[@WilliamsKernel; @WilliamsEDMD.2015].

### Sparsity.

The concept of sparsity is useful to find solutions of underdetermined systems of linear equations. A solution $x \in \mathbb{C}^n$ is considered *sparse*, when most of the entries of $x$ vanish,[@Cands2014MathematicsOS].

Next to offering numerical advantages, a solution with a limited number of terms facilitates the manual interpretation of the results. The problem of finding sparse solutions can be connected to model reduction techniques, as discussed in [@loiseau_brunton_2018].

[@Donoho.2006] showed that an underdetermined system of equations can be converted to a convex problem using a sparsity promoting condition. Consider solving $$y = \theta s^{\prime}$$ for $s^{\prime}$, with $y \in \mathbb{R}^n$, $\theta \in \mathbb{R}^{n \times m}$, $s^{\prime}   \in \mathbb{R}^m$ and $m>n$. The parameter vector $ s^{\prime}$ needs to be determined in order to fit the measurement vector $y$ using the product of the observation matrix $\theta$ and $ s^{\prime}$.

The common $\text{L}_2$-regularized regression promotes a high number of involved terms. Contrariwise, by using the $\text{L}_1$-norm as regularization term, terms of minor impact are neglected from the solution. The $\text{L}_1$-regularized regression (see Eq. ) with weight $\lambda$ can be solved e.g. by the LASSO described by [@TibshiraniLasso]. The solution $s$ is optimal in $\text{L}_1$-sense and sparse within the space spanned by the columns of $\theta$, i.e. the basis functions. $$\label{eq:l1regression}
s = \operatorname*{arg\,min}_{s^{\prime}} \| \theta s^{\prime} - y \|_2 + \lambda \| s^{\prime} \|_1.$$

### Sparse Identification of Nonlinear Dynamics (SINDy).

Since 2015, Nathan Kutz, Steven Brunton and Joshua Proctor have developed the Sparse Identification of Nonlinear Dynamics method[@BruntonSINDy.2016], which is the basis for our data-driven approach to identify hysteretic behavior. First, a library $\Theta(X)$ with possible basis functions is built. SINDy aims to find the minimal amount of those functions, with which a given signal can be approximated. In contrast to DMD, which results in a linear model, SINDy can include arbitrary nonlinearities. As a further development, SINDYc can include feedback and control[@BRUNTON_SINDYc]. A successful application of SINDy in model predictive control is reported in [@Kaiser.].

The main principle of SINDy is the approximation of signals using only few basis functions, i.e. *sparse* solutions. The vector $x(t) = \begin{bmatrix}
x_1(t) &    x_2(t)  & \dots & x_n(t)
\end{bmatrix}^T \in \mathbb{R}^n $ represents all states of the system at time $t$. Measurement data of these states $x(t)$ and their derivatives $\dot{x}(t)$ at the time-steps $t_1, t_2, \dots , t_m$ form two data matrices $X$ and $\dot{X}$, see Eq. , . $$\begin{aligned}
\label{eq:dataMatricesX}
    X &= \begin{bmatrix}
        x^T(t_1) \\ 
        x^T(t_2) \\
        \vdots \\
        x^T(t_m)
    \end{bmatrix}
    =   \newcommand{\mymatrix}[1]{\ensuremath{\overset{\xrightarrow[\hphantom{#1}]{\text{\scriptsize States}}}{#1}\left\downarrow\vphantom{#1}\right.}}
    \mymatrix{\begin{bmatrix}
        x_1(t_1)    &   x_2(t_1)    & \dots &   x_n(t_1) \\
        x_1(t_2)    &   x_2(t_2)    & \dots &   x_n(t_2) \\
        \vdots      &   \vdots      & \ddots&   \vdots   \\
        x_1(t_m)    &   x_2(t_m)    & \dots &   x_n(t_m) \\
        \end{bmatrix}}
        \text{\scriptsize Time}
\\
\dot{X} &=
\begin{bmatrix}
\dot{x}^T (t_1) \\ 
\dot{x}^T (t_2) \\
\vdots \\
\dot{x}^T (t_m)
\end{bmatrix}
= \begin{bmatrix}
    \dot{x}_1(t_1)  &   \dot{x}_2(t_1)  & \dots &   \dot{x}_n(t_1) \\
    \dot{x}_1(t_2)  &   \dot{x}_2(t_2)  & \dots &   \dot{x}_n(t_2) \\
    \vdots      &   \vdots      & \ddots&   \vdots   \\
    \dot{x}_1(t_m)  &   \dot{x}_2(t_m)  & \dots &   \dot{x}_n(t_m) 
    \end{bmatrix}
    \label{eq:dataMatricesXdot}\end{aligned}$$

As a next step, a library $\Theta(X)$ is built from candidate functions, in which $X$ is evaluated. Typical candidates are polynomials, trigonometric functions and constants, but arbitrary functions can be included as well. An example is given in Eq. . $$\begin{aligned}
\Theta(X)\!= \!
\begin{bmatrix}
\vline  &\vline & \vline    &\vline     &       &   \vline  & \vline  &         & \\
1       &   X   & X^{P_2}   & X^{P_3}   & \!\dots\! & \sin(X)   & \cos(X) & \!\dots\!   & e^X  \\
\vline  &\vline & \vline    &\vline     &       &   \vline  & \vline  &         &
\end{bmatrix}
\label{eq:bibliothekStruktur}\end{aligned}$$

The terms $ X^{P_n} $ are polynomials of degree $n$. This includes cross terms of the single states. For the case of quadratic polynomials, this yields to Eq. . $$\delimitershortfall=0pt
\begin{footnotesize}
    X^{P_2} \!= \!
\begin{bmatrix}
x_1^2(t_1) \!&\! x_1(t_1)x_2(t_1) \!& \!\dots\! &\! x_2^2(t_1) &\! x_2(t_1)x_3(t_1) & \!\dots\! &\! x_n^2(t_1) \\
x_1^2(t_2) \!&\! x_1(t_2)x_2(t_2) \!& \!\dots\! &\! x_2^2(t_2) &\! x_2(t_2)x_3(t_2) & \!\dots\! &\! x_n^2(t_2) \\
\vdots     \!&\! \vdots           \!& \!\ddots\!&\! \vdots     &\!  \vdots        & \!\ddots\!&\!\vdots     \\
x_1^2(t_m) \!&\! x_1(t_m)x_2(t_m) \!& \!\dots\! &\! x_2^2(t_m) &\! x_2(t_m)x_3(t_m) & \!\dots\! &\! x_n^2(t_m)
    \end{bmatrix}\label{eq:polynomialLib}
\end{footnotesize}$$

So, the dynamic can be expressed by $\dot{X} = \Theta(X) \Xi$, where $\Theta(X)$ is the evaluated library and $\Xi = \begin{bmatrix}
\xi_1 & \xi_2 & \dots & \xi_n
\end{bmatrix}$ is the coefficient matrix. The coefficients can be determined using sparse regression techniques, e.g. sequential least-squares regression, as has been done in [@BruntonSINDy.2016].

For time-discrete systems, the equation $$\label{eq:sindyDiscrete}
X_{k+1} = \Theta(X_k) \Xi.$$ can be solved accordingly. SINDy allows for easy identification of few dominant functions out of a basis library. Choosing an appropriate library, which contains sufficient functions to describe the dynamics, is essential for a successful identification.

### Criteria for model quality.

There are several methods to compare different models and to decide which one to use. Parameter tuning often uses cross-validation[@kohavi1995study], which bases on the separation of data-sets for training and for validation. The fewer parameters a model has, the less flexible it is, and the less likely it is that overfitting will take place. Therefore the *Akaike Information Criterion* (AIC)[@akaike1974new] was defined to include the number of parameters into the quality of the model: $$\label{eq:aic_definition}
    AIC = 2k - 2\ln(RSS).$$ The integer $k$ gives the number of parameters and the residual sum of squares (RSS) quantifies the prediction error[@ManganModelSelection.2017]. Therefore, the AIC is a trade-off between prediction and model complexity.\
In order to be able to compare errors between different time series, the Mean Absolute Scaled Error *MASE*[@Hyndman.2006] can be used, which relates the error to the overall change of the time series. For a small number of samples, there is a special correction term for the AIC. This so called AICc is defined in Eq. following[@Burnham.2002 Chap. 6.8]. $$\label{eq:aicc}
    AICc = AIC + \frac{2k^2 + 2k}{n-k-1}.$$ The best model achieves the smallest AIC, or AICc respectively. The term $\Delta_j = AIC_i - AIC_{min}$ is a measurement for the difference between models. For $\Delta_j < 2$, models are similar, whereas for $\Delta_j > 10$, the models are significantly different[@KennethBurnham2004]. In context of SINDy, the AIC was applied successfully for model selection[@ManganModelSelection.2017].

Concept for identification of hysteresis-controlled systems {#sec:Concept}
===========================================================

This section provides a technical concept for system identification of hysteretic systems using SINDy, named SINDyHybrid. Based upon theoretical considerations of SINDy and hybrid systems, we developed an approach to incorporate tailored basis functions into the function library. To assess the results, we define criteria for model quality. SINDy has been used by [@BruntonSINDy.2016] to model the dynamics of systems with a bifurcation parameter $\mu$, which changes the qualitative behavior of the system’s dynamics. In order to incorporate $\mu$ into the observables, the equation $\dot{x} = f(x)$ was extended by $\mu$, such that $$\begin{aligned}
    \dot{x} &= f(x, \mu), \label{eq:sindyBifurcation}\\
    \dot{\mu} &= 0. \label{eq:MuDerivative}\end{aligned}$$ The function library $\Theta$ then consists of functions of $x$ and $\mu$. The value of $\mu$ however remains unchanged, as indicated by Eq. . A more general formulation has been used in [@BRUNTON_SINDYc], where feedback control signals were included by choosing $$\dot{\mu} = g(x). \label{eq:MuDerivativeSindyC}$$

Expanding the theory for hybrid systems
---------------------------------------

We now further generalize the approach of Eq. by choosing $$\dot{\mu} = g(x, \mu), \label{eq:MuDerivativeGeneral}$$ which connects changes of $\mu$ to current values of $\mu$. This represents the behavior of a hybrid system. Special cases of Eq. can be linked to the previously presented methods of modeling hybrid systems. Choosing $$\dot{\mu} = g(t)$$ results in the formulation of SARX models (see \[SARX\]). In cases where the change of the qualitative behavior of the system depends on only current values of the state $x$, the change of $\mu$ can be represented by $$\dot{\mu} = g(x(t)), \label{eq:MuDerivativePWARX}$$ which is the premise of PieceWise-Affine (PWA) models (see \[PWARX\]). For a linear function $g$, the regions for each state have a linear border. For a nonlinear $g$, complex borders can be modeled between regions, as is the case for NPWA models. If the function $f$ in Eq. is linear, the formulation is equivalent to PWARX models, otherwise it results in PWNARX models.\
Hysteresis models not only depend on the current state of the system, but also on past values up to a horizon $q$. Let $$\dot{\mu} = g(x(t-q), x(t-q+1), ..., x(t-1), x(t), \mu) = g(x_q, \mu).$$ This way, $\mu(t)$ serves as a storage for information of past values of the system. The general formulation Eq. is equivalent to $$\dot{\mu} = g(x(t), \mu(t)). \label{eq:MuDerivativeSimple}$$

Therefore it is sufficient to model the change of $\mu$ by current values of $x$ and $\mu$, even for hysteretic systems. The question of how to choose a proper function $g$ will be addressed next.

Modeling discrete states {#modelingDiscreteStates}
------------------------

In order to find a suited function $g$ from Eq. to model changes of discrete systems, one has to consider the way the change of states takes place. The transition can for example be either abrupt or smooth, following a certain shape. In the case of two-point controlled systems the state transition is abrupt: As soon as a certain threshold gets passed, the system behavior changes more or less instantaneously. This assumption holds for many hybrid systems.

In \[hysteresisModeling\] we discussed the Preisach model for identification of hysteresis, which uses the relay operator $R_{\alpha, \beta}$. Based upon the idea that state transitions in many hybrid systems happen abruptly, the relay hysteron seems to be a suited candidate for the function $g$.

Preparing the data
------------------

The first step of data processing is affine scaling of all data to a sensible and comparable range, e.g. from -1 to 1. This promotes an equal treatment of all time series, not depending on the absolute values of the data.

Before identification with SINDy begins, the possible state indicators (relay hysterons, Section \[modelingDiscreteStates\]) have to be computed from the data. This is done in three steps: Building indicator functions, pairwise matching of indicator functions and evaluating the relay hysterons.

We assume a change of the system behavior as soon as a value $x$ rises above $\beta \in \mathbb{R} $ or falls below $\alpha \in \mathbb{R} $. We also assume $\alpha$ and $\beta$ are included in our data. The first step is to build differences between all time series in our data. If the units of the time series are known, it is helpful only to build differences between commensurable data, e.g. time series which share the same unit, as exemplarily done in Eq. . $$\begin{aligned}
        \label{eq:differences}
        \tilde{x}_1 &= x-\alpha \\
        \tilde{x}_2 &= x-\beta \\
        \tilde{x}_3 &= \alpha - \beta
    \end{aligned}$$

We then create indicator functions masking these differences depending on the sign of the value. $I(\tilde{x}_2)$ masks ranges, where $x \geq \beta$. $\bar{I}(\tilde{x}_1)$ on the other hand masks ranges, where $x < \alpha$. $$\begin{aligned}
\begin{split}
I(x) &= \begin{cases}
1,  & \text{if $x \geq 0 $} \\
0,  & \text{if $x < 0 $} \\
\end{cases} \label{eq:indicator1} \\
\bar{I}(x) &= \begin{cases} 
1,  & \text{if $x < 0 $} \\
0,  & \text{if $x \geq 0 $} \\
\end{cases}
\end{split}\end{aligned}$$

A hysteron can be formed out of two indicator functions, which are not both true at the same time. For all matching indicator pairs $I_\alpha, I_\beta$, the hysteron $H_{\alpha, \beta}$ is built. This way, $I(\tilde{x}_2)$ and $\bar{I}(\tilde{x}_1)$ would also be combined. For each hysteron, a complementary hysteron $\bar{H}_{\alpha, \beta}$ will also be built. The evaluation is done according to

$$\label{eq:hysteron}
H(k)_{\alpha, \beta} = \begin{cases}
1, & \text{if $I_\beta = 1$.} \\
0, & \text{if $I_\alpha = 1$.} \\
H(k-1), & \text{for $I_\beta = I_\alpha = 0$.}
\end{cases}$$

for each step $k>1$.

This is equivalent to the formulation of the relay hysteron in Eq. . Presuming alternating switches between the states, $H(1)_{\alpha, \beta}$ can be initialized with the opposite state of the first occuring switch. This presumption has to be reviewed depending on the use case.

Robust handling of transitions – the proximity hysteron {#sec:proximity}
-------------------------------------------------------

A problem in the design of hysterons is the fact, that the critical threshold values $\alpha$ and $\beta$ from Eq. are only reached for a short moment. Due to a slow sampling frequency or noisy data it may happen that these critical points are not even included in the data at all. In order to improve robustness of the relay hysterons, we developed a “proximity-hysteron", which utilizes an $\epsilon$-range around the actual threshold value for switches between states. The evaluation is done according to $$\label{eq:proxHysteron}
H_{\epsilon}(k)_{\alpha, \beta} = \begin{cases}
1, & \text{if $I_{\epsilon, \beta}  = 1$.} \\
0, & \text{if $I_{\epsilon, \alpha} = 1$.} \\
H_{\epsilon}(k-1), & \text{for $I_{\epsilon, \beta} = I_{\epsilon, \alpha} = 0$.}
\end{cases}$$

with the proximity indicator functions $$\begin{aligned}
\begin{split}
I_{\epsilon_\beta}(x) &= \begin{cases}
1,  & \text{if $x = \min\{\max( Z_{x_\beta}(t) ), \beta\} $.} \\
0,  & \text{otherwise.} 
\end{cases} \\
I_{\epsilon_\alpha}(x) &= \begin{cases}
1,  & \text{if $x = \max\{\min( Z_{x_\alpha}(t) ), \alpha\} $.} \\
0,  & \text{otherwise.} 
\end{cases}
\end{split}\end{aligned}$$

and the connected sets and\
$Z_{x_\alpha}(t) = \{x \mid x < \alpha + \epsilon_\alpha\}$, with $\epsilon_\beta,\epsilon_\alpha \in \mathbb{R} > 0$, which include $x(t)$.

The proximity indicator functions define an $\epsilon$-range, which moves the threshold for the transition in the corresponding direction (see Fig.\[fig:proximity\_hysteron\]). However, the transition between states does not happen as soon as the shifted threshold is reached, but only when the extremum within the connected set of points within the $\epsilon$-range is reached. An exception occurs, when the actual threshold $\alpha$ or $\beta$ is reached. In this case, the switch is done regardless of whether this is the extremum within this set.

![Switches of the proximity hysteron[]{data-label="fig:proximity_hysteron"}](proximityHysteronSVG){width="\linewidth"}

Due to the introduction of the proximity hysteron, state transitions can even be localized under harsh conditions – if the threshold points do not occur in the measurement data. The connection between signal, threshold values, indicator functions and hysterons is visualized in Fig.\[fig:buildingOfHysterons\].

![Interaction between signal, threshold values, indicator functions and hysterons[]{data-label="fig:buildingOfHysterons"}](simSignal "fig:"){width="\linewidth"}\
![Interaction between signal, threshold values, indicator functions and hysterons[]{data-label="fig:buildingOfHysterons"}](simIndicatorsSVG "fig:"){width="\linewidth"}\
![Interaction between signal, threshold values, indicator functions and hysterons[]{data-label="fig:buildingOfHysterons"}](simHysteronSVG.pdf "fig:"){width="\linewidth"}

System Identification with SINDy
--------------------------------

Once the hysterons are built and other preprocessing measures have been taken, the SINDy algorithm can start. The schematics are shown in \[fig:sindyFlow\]. Creating the data matrices includes choosing the right state vector $X$. The past values $H(k-1)_i$ for $i = 1, ..., m$ (see Eq. ) of all $m$ hysterons have to be included, holding information about the system’s current discrete state. Past measurement values $x(k-1), ..., x(k-q) $ up to a defined horizon $q$ can also be considered, allowing a representation resembling an FIR-filter. The state at time $k$ can be expressed by

$$\begin{aligned}
\begin{split}
    X(k) =& [x_1(k), ... , x_n(k), H_1(k-1), \bar{H}_1(k-1), ...
    \\& \bar{H}_m(k-1), x_1(k-1), ... , \bar{H}_m(k-q)]
\end{split}\end{aligned}$$

![The SINDy algorithm with tailored basis functions for identification of hysteretic systems[]{data-label="fig:sindyFlow"}](sindyFlowSVG){width="0.6\linewidth"}

As for the design of the library $\Theta$, polynomial basis functions are useful and easy to include. The function $\varphi_{poly}(x_1, x_2, ..., x_n)$ evaluates all polynomials between $x_1$ to $x_n$ up to a predefined polynomial degree. Hysterons are evolved by $\varphi_{relay}(x_1, x_2, ..., x_n, H_i)$. After evaluating all basis functions $\phi$, all basis functions that are unaffected by hysterons are multiplied with the updated hysterons. The finished library $\Theta$ consists of three parts: One without hysterons, the second with cross terms of hysterons and non hysteretic terms and the third of only the updated hysterons. Once the library is evaluated, a sequential least squares regression with a tuning parameter $\lambda$ is used to solve $ X(k+1) = \Theta(X(k)) \Xi$ for the sparse coefficients $\Xi$, as proposed in[@BruntonSINDy.2016].

For a library consisting of only polynomial basis functions and tailored basis functions for the propagation of hysterons, the evaluation and calculation of the next step is shown in Fig. \[fig:sindyLib\].

![Building and evaluation of the library[]{data-label="fig:sindyLib"}](bibGrafikSVG){width="\linewidth"}

Practical evaluation and results {#sec:practicalEvaluation}
================================

In order to evaluate the presented concept, we used different model parameters and differing data quality to identify a simple system. Our aim was to determine under which conditions the proposed identification procedure succeeds in finding a suited model, which parameters influence model quality, and where the limitations lay. The influence of the size of the library is studied in the first experiment by varying the degree of polynomials used. The robustness to sample rate and measurement noise are further aspects under test.

Example: Hysteresis-controlled tank system
------------------------------------------

The level of the basin in Fig. \[fig:Wassertank\] is controlled to a certain height by a *two-point controller* with an upper set point $h_{max}$ and a lower set point $h_{min}$. A leak causes a constant drain $Q_{out} \in \mathbb{R} < 0$, which neglects the influence of the current water level. The controller with control signal $u(t)$ turns the pump on whenever $h_{min}$ is reached, resulting in a steady inflow $Q_{in} \in \mathbb{R} > 0$. The pump is turned off when the level raises over $h_{max}$.

\[tank\_model\]

The length of the supply pipe $l_P$ can be used to analyze the effect of delay.

The simulation is implemented in Simulink^^. For this example, 20 different settings were processed, whereof 16 served as training- and four as validation set. The dynamics of the tank system is given in Eq. . $$\begin{aligned}
\begin{split}
\label{eq:tank_dynamics}
    h(t+1) = h(t) + u(t-l_p) + Q_{out} \text{, with } \\
    u(t) =  \begin{cases}
    Q_{in}, & \text{if pump on.} \\
    0, & \text{if pump off.}
    \end{cases}
\end{split}\end{aligned}$$ As a result, the pump works in rectangular pulses and the water level changes in a saw-tooth curve.

Influence of polynomial degree
------------------------------

The library prepared for SINDy consists of polynomials and hysterons. In a first experiment, we vary the degree of the polynomials. Figure \[fig:exp2\_sim\_h\] contains the simulated prediction of the filling level. Apparently, polynomials of low degree lead to the best results. For polynomials of degree one, the identified model is shown in Eq. . It consists of 4 out of 55 functions of the library, meaning a high degree of sparsity. Due to the hybrid form of hysteron $H_1$, it has the same structure as the correct model in Eq. . The models of higher degree manage to reproduce state transitions at the correct critical points, but they induce an inaccuracy for the slope.\
$$\begin{aligned}
        \tilde{h}_{poly1}(k+1) &= h(k)  + 1.56 \times 10^{-4}(4.98 Q_{in}  (1 - H_1) + Q_{out} ) \label{eq:poly1_model}, \nonumber \\
        H_{1}(k) &= \begin{cases}
        1, & \text{if } h > h_{max}. \\
        0, & \text{if } h < h_{min}. \\
        H_{1}(k-1), & \text{else}.
        \end{cases} \end{aligned}$$

![image](polydegSimHSVG){width="\linewidth"} \[fig:exp2\_sim\_h\]

Influence of sample rate
------------------------

As a further objective, we analyzed the effect of the sampling rate. As explained in Part \[sec:proximity\], this method cannot localize transitions if the critical points are not included in the given data, which may happen for slow sampling rates. Figure \[fig:exp3\_normal\_sim\_h\] shows, that state transitions are not detected using standard hysterons for slower sampling than $0.02$ s.

![image](samplingrateNormalHystSVG){width="\linewidth"} \[fig:exp3\_normal\_sim\_h\]

As a remedy, we presented the proximity hysteron in Section \[sec:proximity\]. The effect can be seen in Fig. \[fig:exp3\_proxi\_sim\_h\]. With this extension, the algorithm achieves good results for all sampling rates. The best model is shown in Eq. and was built from a library of 23 basis functions. It is equivalent to the model in Eq. , as $1 - H_1 = \bar{H}_1$.\
$$\begin{aligned}
    \tilde{h}_{prox, 0.01}(k+1) &= h(k)  + 1.56 \times 10^{-4}(4.98 Q_{in} \ \bar{H}_1 + Q_{out} ) \label{eq:snr_model}, \nonumber \\
    H_{1}(k) &= \begin{cases}
        1, & \text{if } h > h_{max}. \\
        0, & \text{if } h < h_{min}. \\
        H_{1}(k-1), & \text{else}.
        \end{cases} \end{aligned}$$

![image](samplingrateProxHystSVG){width="\linewidth"} \[fig:exp3\_proxi\_sim\_h\]

Effect of noisy data
--------------------

Beside an insufficient sampling rate, noise is a typical problem for data-driven approaches. In this experiment we investigate its influence with the SNR values 1000, 100, 50 and 10. The results are depicted in Fig.\[fig:exp4\_sim\_h\].

![image](snrSimHSVG){width="\linewidth"} \[fig:exp4\_sim\_h\]

It can be seen that noise influences the identifcation quality very harshly. The approach works properly for an SNR of 1000. For an SNR of 100, the slope of the dynamics is not modeled accurately, resulting in a phase shift. The model for SNR 50 apparently identified two discrete states, but the state transitions as well as the dynamics within the states are faulty. For SNR 10, the measurements are too noisy to identify discrete states. As such, the model tries to satisfy all measurements by converging to the average.

Discussion {#sec:discussion}
==========

The results show, that SINDy with taylored basis functions in form of hysterons (SINDyHybrid) is a suitable method for the identification of hysteresis-controlled systems. The continuous system dynamics and the discrete states do not have do be modeled separately. With given signals and guesses for the thresholds, the hysteresis relation is found automatically. Adequate models were identified under varying sample rate and signal noise ratio. Designing a suited library remains the main challenge for the user. It was shown that systemic knowledge can be included in the function library, and that specific adaptations, e.g. with regard to poor data quality, can be integrated. The library should not be larger than necessary, as this complicates the regression and hinders identification of the correct basis functions. The resulting models are sparse and thus easily interpretable. This enables testing for plausibility, which is required for many industrial applications. However, this method may show problems for systems with many discrete states. For each state, a hysteron will be identified, which is multiplied with functions describing the dynamics for that particular state. This results in a solution with a high number of overall terms, which is not according to the idea behind SINDy that the optimal solution is sparse.

It must be noted, that creating the hysteron requires knowledge about the actual threshold value of the state transition. For two-point controlled systems, those threshold values can usually be read from data. The $\epsilon$-range around the threshold values used by the proximity hysteron can compensate inaccuracies of the threshold values, and this way even compensate aging effects.

Conclusion {#sec:conclusion}
==========

The presented work proposes a novel data-driven approach for identification of hybrid systems, SINDyHybrid. This method is exemplarily applied to a hysteresis-controlled system. The flexibility of the SINDy framework enables the integration of tailored basis functions, which can incorporate domain-specific knowledge. Based upon the relay-operator of the Preisach model, we developed the proximity hysteron, which allows for robust identification of state transitions and dynamics.

[^1]: The presented work was partially performed in the project *Reinforcement Learning for complex automation technology*, supported by Investitionsbank Berlin (IBB), co-financed by the European Regional Development Fund.
