Computationally efficient identification of continuous-time Lur'e-type systems with stability guarantees

In this paper, we propose a parametric system identification approach for a class of continuous-time Lur’e-type systems. Using the Mixed-Time-Frequency (MTF) algorithm, we show that the steady-state model response and the gradient of the model response with respect to its parameters can be computed in a numerically fast and efficient way, allowing efficient use of global and local optimization methods to solve the identification problem. Furthermore, by enforcing the identified model to be inside the set of convergent models, we certify a stability property of the identified model, which allows for reliable generalized usage of the model also for other excitation signals than those used to identify the model. The effectiveness and benefits of the proposed approach are demonstrated in a simulation case study. Furthermore, we have experimentally shown that the proposed approach provides fast identification of both medical equipment and patient parameters in mechanical ventilation and, thereby, enables improved patient treatment.


Introduction
Accurate dynamical models of complex systems are required for model-based controller design, analysis of the dynamic behavior of the system under study and improved system design.System identification deals with the construction of such dynamic models from observed input and output data.For the class of linear systems, many identification techniques exist (Ljung, 1987), together with user-friendly software, see, e.g., Ljung (1995).However, many physical systems are essentially nonlinear and, depending on the application, one can choose to model it in a linear framework by neglecting (non-significant) nonlinear phenomena, or can choose to also model the nonlinearities and, therewith, increase the model accuracy at the cost of model complexity.
A major challenge in system identification is to enforce a form of stability on the identified model (Ljung, 2010;Schoukens & Ljung, 2019;Tobenkin, Manchester, & Megretski, 2017).Even if it is known that the true system under study enjoys some stability property, it is not guaranteed that the identified model preserves this stability property due to the finite length of the data set, noisy data sets or the presence of unmodeled dynamics.Several methods that enforce exponential stability are proposed in the literature for the identification of linear time-invariant (LTI) models (Hakvoort & Hof, 1994;Lacy & Bernstein, 2003;Tobenkin et al., 2017;Umenberger & Manchester, 2018;Van Gestel, Suykens, Van Dooren, & De Moor, 2001).For nonlinear systems, however, the literature on identification of nonlinear models with some guaranteed form of stability is scarce and the consequences of not enforcing stability are more pronounced and more complex in the nonlinear context (Decuyper, De Troyer, Runacres, Tiels, & Schoukens, 2018;Decuyper, Runacres, Schoukens, & Tiels, 2020;Schoukens & Ljung, 2019).Firstly, it is well-known that nonlinear models can exhibit multiple stable solutions being attractive for different sets of initial conditions (Khalil, 1996).Consequently, the response of the identified model might be close to the measured system response, but might as well be far off or even become unbounded, depending on the initial condition.Secondly, nonlinear models can have a large sensitivity to the excitation signal, and, therefore, their response can be significantly different, or even become unbounded, even for the slightest change in the excitation signal (Khalil, 1996).Such models have poor generalization capabilities to other excitation signals than the ones used during identification.With the instability issue of the identified model in mind, Tobenkin et al. (2017) and Umenberger and Manchester (2018) developed a method that imposes global incremental stability for a class of identified black-box nonlinear state-space models, where the incremental stability property avoids the above problems of high sensitivity to input changes and initial conditions.In that approach, the data set includes the state sequence of the underlying nonlinear system, which is non-trivial and rather case specific to obtain in practice.Other types of stability properties, e.g., asymptotic stability of the origin for zero inputs and input-to-state stability, are enforced in Doyle, Pearson, and Ogunnaike (2002), Milanese and Novara (2005), Suykens, de Moor, and Vandewalle (2000), Suykens, Vandewalle, and De Moor (1997) and Umlauft, Lederer, and Hirche (2017), mostly for autonomous black-box nonlinear models.Note, however, that models that enjoy these other types of stability properties can still exhibit a large sensitivity to input changes and initial conditions.
To address the challenge of identifying nonlinear models with a form of incremental stability, we focus on a practically relevant class of nonlinear systems, namely Lur'e-type systems (Khalil, 1996).These systems are block structured in the sense that all the LTI dynamics are captured in an LTI block and separated from the nonlinearities that are captured in a static nonlinear block in the feedback loop.The class of Lur'e-type systems encompasses so-called Wiener, Hammerstein and nonlinear feedback systems and, therefore, captures a large class of nonlinear systems.In the system identification literature, Lur'e-type systems are also referred to as nonlinear feedback systems (Giri & Bai, 2010), and NonLinear Linear Fractional Representation (NL-LFR) systems (Schoukens & Tóth, 2020;Vanbeylen, 2013).
Different approaches for the identification of Lur'e-type systems have been proposed in literature (Hu & Ding, 2013;Mulders, Vanbeylen, & Usevich, 2014;Paduart, Horvath, & Schoukens, 2004;Schoukens, Nemeth, Crama, Rolain, & Pintelon, 2003;Schoukens & Tiels, 2017;Schoukens & Tóth, 2020;Van Mulders, Schoukens, & Vanbeylen, 2013;Van Pelt & Bernstein, 2001;Vanbeylen, 2013;Wang, Shen, Wu, & Ji, 2016).What most of these approaches have in common is that, as a final step in the identification procedure, a non-convex cost function has to be minimized using gradient-based optimization routines.Furthermore, all these approaches are formulated in a discrete-time setting.Although continuous-time counterparts of most of these methods can be formulated readily, it is an open problem how to perform the gradient-based minimization of the cost function in the continuous-time setting in a numerically efficient way.More specifically, any gradient-based optimization routine requires the computation of model responses to evaluate the cost function and requires the computation of the gradient of the cost function with respect to the model parameters.While model simulation of discrete-time nonlinear models boils down to a series of computationally cheap algebraic operations, the standard means of model simulation for continuous-time nonlinear models is numerical forward integration, which is computationally extremely expensive.
In this paper, we present a parametric identification approach for a class of continuous-time single-input single-output (SISO) Lur'e-type systems that: (i) guarantees the identified model to be exponentially convergent and, (ii) uses numerically efficient tools to compute model responses and gradient information to be used for gradient-based optimization.Point (i) guarantees that the identified model is globally exponentially convergent, which is a property of (non)linear models that guarantees the uniqueness, boundedness and global exponential stability of the steady-state solution (Pavlov, van de Wouw, & Nijmeijer, 2006).As a consequence, for any bounded excitation, the identified model 'forgets' its initial condition and, therefore, exhibits a uniquely defined, bounded steady-state solution, which is globally exponentially stable and depends solely on the applied excitation.Furthermore, the identified model enjoys the property that a small variation of the excitation signal results in a small variation of the steadystate output, which adds to the robustness of the identification result.Point (ii) enables system identification of continuous-time nonlinear models in a computationally efficient way.Hereto, the so-called mixed-time-frequency (MTF) algorithm is employed to enable fast computation of steady-state model responses under periodic excitations (Pavlov, Hunnekens, Wouw, & Nijmeijer, 2013).The MTF algorithm facilitates usage of global optimization techniques (Locatelli & Schoen, 2013) and also enables efficient computation of the gradient of the underlying cost function with respect to the model parameters, which can be effectively used in any gradient-based optimization routine.
We demonstrate the performance and the computational benefits of the proposed identification approach in two case studies.The first case study is a simulation example, where we identify a black-box nonlinear model.We show that the proposed identification strategy is significantly faster than conventional methods as the computation time is reduced from hours to only minutes.The second case study is an experimental study on mechanical ventilation, which is used to regulate the breathing of patients in respiratory distress during nursery of intensive-care and acutely ill surgical patients.Particularly, during the COVID-19 crisis, mechanical ventilation is extensively used to support breathing of hospitalized patients with severe lung damage.In our longstanding relation with Demcon Macawi Respiratory Systems, Best, The Netherlands, we mainly focused on the pressure control aspect (Hunnekens, Kamps, & van de Wouw, 2018;Reinders, Hunnekens, Heck, Oomen, & van de Wouw, 2020).In this paper, we address the identification of the parameters of a first-principle model.The model includes the patient's lung parameters, which reveal important patient health information that is used in the medical decision-making process (Amato et al., 2015;Lachmann, 1992) and for pressure control purposes (Hunnekens et al., 2018;Tehrani, Rogers, Lo, Malinowski, Afuwape, Lum, Grundl, & Terry, 2004).We show that our identification method significantly reduces the computation time, which is crucial as it enables faster patient treatment with the aim to avoid negative consequences for the patient's lungs (Gajic, Dara, Mendez, Adesanya, Festic, Caples, et al., 2004;Slutsky & Ranieri, 2013) and, saves valuable time of the medical practitioner.
In summary, the main contribution of this paper is a computationally efficient approach to the identification of continuoustime convergent Lur'e-type systems that guarantees that the identified model is also globally exponentially convergent.We have demonstrated theoretically that this approach is statistically consistent and we have demonstrated numerically that this approach is effective and computationally efficient.Furthermore, we have shown that the developed approach is an enabler for innovation in mechanical ventilation as it allows to identify physical parameters in a fast way, which can subsequently be used for improved patient treatment.
A preliminary part of this work is presented in Shakib, Pogromsky, Pavlov, and van de Wouw (2019).Compared to Shakib et al. (2019), the current paper considers an extended model class and a more generic model parametrization.Furthermore, this paper includes a statistical consistency analysis of the estimator, a detailed overview of methods to obtain an initial convergent model, a novel simulation case study and an experimental validation of the proposed identification approach.
The remainder of this paper is structured as follows.Section 2 introduces the considered class of Lur'e-type systems and recalls sufficient conditions for convergence.Section 3 formally poses the identification problem as a constrained optimization problem and also shows the statistical consistency of the estimator.Section 4 provides an overview of methods to minimize the cost function, constrained to the set of convergent models.Section 5 introduces the MTF algorithm and also shows how the gradient of the cost function is computed accurately and efficiently.Section 6 presents the performance of the identification approach in a simulation study.Section 7 describes the mechanical ventilation case study.Section 8 closes with concluding remarks.

Convergent Lur'e-type systems
We consider SISO Lur'e-type systems described by where x(t) ∈ R n is the state vector, y(t) ∈ R is the unmeasured feedback signal, z 0 (t) ∈ R is the output, w(t) ∈ R is the user-defined external input, φ(y(t), w(t)) : R × R → R is a static nonlinear function.We assume that φ(0, w) = 0 ∀ w ∈ R. Fig. 1 depicts the considered Lur'e-type system schematically.Although checking stability of forced LTI systems is well-understood, checking stability of (forced) solutions of nonlinear systems is non-trivial.Consequently, enforcing stability on the identified nonlinear model becomes challenging (Ljung, 2010).To address this challenge, we use the notion of convergent dynamics defined as follows.
Definition 1 (Pavlov et al., 2006).System (1) with input w(t) that is defined and bounded on t ∈ R, is said to be globally exponentially convergent if there exists a solution x(t) satisfying the following conditions: • x(t) is defined and bounded on t ∈ R, • x(t) is globally exponentially stable.The solution x(t) is called the steady-state solution and depends on the applied excitation signal w(t).For the considered system (1), the following theorem provides sufficient conditions for global exponential convergence.Here, vertical bars are used to denote the Euclidean norm.
Theorem 1 (Pavlov et al., 2013).Consider system (1) and suppose that C1 The matrix A is Hurwitz; C2 There exists a K > 0 such that the nonlinearity φ(y, w) Then, for any input w(t) that is defined and bounded on t ∈ R, system (1) is globally exponentially convergent according to Definition 1.Moreover, any other solution x(t) of system (1) subject to the same input w(t) satisfies Constants α > 0 and β > 0 are independent of the input w(t).
For notational convenience, we write exponentially convergent instead of exponentially convergent for any bounded input w(t) from here on.The conditions of Theorem 1 are used in Section 3 to enforce the identified Lur'e-type model to be globally exponentially convergent according to Definition 1. Globally exponentially convergent systems 'forget' their initial condition since, independent of the initial condition, all solutions x(t) converge exponentially to the steady-state solution x(t).The conditions of Theorem 1 give the even stronger property of inputto-state convergence, which implies that small variations (in the infinity norm) of the input signal w(t) lead to small variations (in the infinity norm) of the steady-state solution x(t), see Definition 3 in (Pavlov et al., 2006).Convergent systems also have the property that when excited by a periodic excitation signal w(t) with period-time T , the steady-state solution x(t) is also periodic with the same period-time T (Pavlov et al., 2006).In the literature, often the PISPO (Periodic Input, the same Period Output) system class is considered, which also possesses this property (Schoukens, Pintelon, Dobrowiecki, & Rolain, 2005;Schoukens & Tiels, 2017).This latter property facilitates the use of only steady-state data for identification, which is common practice in nonlinear system identification (Schoukens & Ljung, 2019).

Identification setting
This section formally introduces the identification problem.Hereto, we first introduce the considered model class.After that, we introduce the cost function central in the identification problem and show that the proposed estimator is statistically consistent.

Model parametrization
The considered model class is a copy of (1): parametrized by the vector θ ∈ R n θ .The state dimension n of the true system (1) is considered known.
We distinguish between two parametrization approaches, namely (i) the model parametrization following from firstprinciple modeling; and (ii) a black-box model parametrization.In approach (i), the model parameters of both the linear and nonlinear parts have a physical meaning.Such model parametrization is used in the experimental case study in Section 7. In approach (ii), the parameters do not have a physical meaning.Then, for the state-space matrices, a canonical form or a full parametrization can be taken, in which each element of the state-space matrices is one model parameter.In a black-box approach, the nonlinearity ϕ is typically parametrized by a set of basis function as follows: where n LTI is the number of parameters used to parametrize the LTI block.The user-defined basis functions f k (y(t, θ), w(t)) have associated parameters θ k .A black-box model parametrization is used in the simulation case study presented in Section 6.
One of the goals of this work is to guarantee that the identified model is globally exponentially convergent according to Definition 1.To this extent, we define the set Θ ⊂ R n θ with a nonempty interior as the set of parameter vectors θ for which model ( 2) satisfies conditions C1-C3 in Theorem 1 (with A, B, C, φ(y, w) replaced by A(θ ), B(θ ), C (θ ), ϕ(y, w, θ ), respectively).The set Θ serves as a constraint on the parameter vector θ to be respected in the identification process.We close this section by assuming that the to-be-identified true system (1) is also globally exponentially convergent, which is a stability property that the developed identification method aims to preserve for the identified model.
Assumption 1.The true system ( 1) is in the model class proposed in (2) with parameter vector θ 0 ∈ intΘ, and satisfies conditions C1-C3 of Theorem 1.Therefore, the true system (1) is globally exponentially convergent according to Definition 1.

Remark 2.
The nonlinearity ϕ of model ( 2) should be parametrized carefully in a black-box parametrization approach.In particular, the incremental sector bound imposed by condition C2 in Theorem 1 should be satisfied globally, i.e., for all y 1 , y 2 , w ∈ R. Inappropriate choices for the basis functions f k in (3) can lead to unstable models and to models that exhibit multiple co-existing (stable) steady-state solutions.This is encountered in, e.g., polynomial basis functions (Decuyper et al., 2018).Suitable nonlinearities are, e.g., piecewise linear maps, sigmoids, the arctangent function and the hyperbolic tangent function.

Cost function & consistency
Assumption 1 guarantees that the true system is globally exponentially convergent.Therefore, when excited by a periodic input w(t) with period-time T , the periodic steady-state output z0 (t) of system (1) has the same period-time T .In practice, only a discrete version of the measured steady-state output z0 (t) is available at sampling instants t k := t 0 + kt s with k = 1, . . ., N, where t s is the sampling interval and N = PT /t s is the total number of samples in P steady-state periods.We assume the following noise scenario.
The discrete-time noise source e(t k ) is zero-mean Gaussian white noise with finite variance σ 2 e and independent of the applied input signal w(t k ).
By guaranteeing that θ ∈ Θ, i.e., the candidate model is also globally exponentially convergent, model (2) subject to the same periodic input w(t) results in the periodic steady-state output z(t, θ) with the same period-time T .This leads to the definition of the steady-state error over P periods: The squared simulation error is then taken as cost function, which is defined as follows: (5) Given the data set {w(t k ), z(t k )} N k=1 , the identification objective is to find the parameter vector θN that globally minimizes J N (θ )   constrained to the set of convergent models characterized by the parameter set Θ: The following persistency of excitation assumption is required to study consistency of the estimator (6).

Assumption 3. The unique global minimum of
The function V N (θ ) can be interpreted as the steady-state mismatch between the model response z(t k , θ) and the noiseless system response z0 (t k ).Assumptions 1-3 allow us to formulate the following consistency result.

Theorem 2.
Under Assumptions 1-3, the estimate θN in (6) converges in probability to θ 0 , i.e., Proof.The proof can be found in Appendix A. ■ Remark 3. The model parametrization is non-unique (Schoukens, Gommé, Van Moer, & Rolain, 2008;Schoukens & Tiels, 2017) as a model subject to a similarity transformation of the LTI block, a gain exchange between the LTI block and the static nonlinear block or a loop-transformation (Khalil, 1996) produces the same steady-state output.Therefore, the true parameter vector θ 0 is a set and the consistency claim of Theorem 2 should be understood in the sense that θN converges in probability to the set θ 0 as N → ∞.

Cost function minimization strategy
The identification problem defined in ( 6) is a constrained optimization problem with the cost function in (5).As this cost function is generally non-convex, we propose a two-step optimization approach to solve (6): (i) initialization and (ii) gradient-based optimization.The result of step (i) should be a set of initial model parameters θ init in the vicinity of θ 0 , which globally minimizes J N (θ ).For applicability of the numerically efficient tools that we present in Section 5, the initial model with parameters θ init should satisfy the conditions of Theorem 1 and, therefore, be globally exponentially convergent.The gradient-based search of step (ii) is started from the set of initial parameters θ init and results in the parameter set θN that corresponds to nearest minimum of the cost function J N (θ ).Again, for applicability of the numerically efficient tools that we present in Section 5, in all optimization iterations of the gradient-based search, the conditions of Theorem 1 should be satisfied.

Initialization
Let us discuss three approaches for the initialization step that yield a set of initial model parameters θ init .Again, we want to stress that the model with parameters θ init should satisfy the conditions of Theorem 1.

Physical insight
This method is only applicable for models that are derived by first-principle modeling, implying that the parameters of these models represent physical quantities.The user appropriately chooses the initial parameters θ init based on its physical insights.However, it should be kept in mind that the model with parameters θ init should satisfy the conditions of Theorem 1, which can make the initial guess challenging, especially for models with a large number of parameters.
Best Linear Approximation This method relies on estimating the so-called Best Linear Approximation (BLA) of the nonlinear system (Schoukens et al., 2005).Hereto, the MATLAB implementations N4SID or TFEST, can be used, which can also enforce the A(θ )-matrix of the identified model to be Hurwitz.In addition, by imposing no feedback, i.e., choosing C (θ ) and D(θ ) both 0, B(θ ) = 0 or ϕ(y, w, θ ) = 0, the initial model satisfies the conditions of Theorem 1.The BLA framework has the advantages that it allows for fast estimation of initial model parameters, it provides a nonlinear distortion analysis, and it gives a rough estimate for the state dimension.
Global optimization This method employs global optimization routines such as genetic-type, swarm intelligent or Monte-Carlo routines (Locatelli & Schoen, 2013).Although global optimization is typically computationally expensive, we can achieve fast operation using the efficient numerical tools that we introduce in Section 5. Global optimization routines typically cluster around the global minimum and provide, therefore, an accurate set of initial model parameters θ init that correspond to a full nonlinear model rather than a linear model as in the BLA framework.To limit the search space, (rough) bounds on the parameters values should be known a priori.
Depending on the application at hand, either one of the three methods detailed above can be used to find a set of initial model parameters θ init .In the simulation case study in Section 6, we demonstrate the effectiveness of a global optimization routine.In the experimental case study in Section 7, we use available knowledge on the physical parameters for initialization.Numerous identification results with the BLA method are presented in the literature, see Schoukens and Tiels (2017) for an overview.
Remark 4. Global optimization can also be viewed as a standalone identification procedure without a local gradient-based search as a second step.To view it as a two-step approach, consider step one as a less accurate but global search and step two as a local refinement.

Gradient-based optimization
After an initial parameter vector θ init in the initialization step is obtained, a gradient-based search can be performed.Hereto, the gradient of the cost function in (5) with respect to the parameter vector θ is required: In ( 7), the steady-state error ε(t k , θ) defined in (4), as well as the gradient of this steady-state error ∂ ε(t k , θ)/∂θ = ∂ z(t k , θ)/∂θ , are required.We show in the next section that we can compute both ε(t k , θ) and ∂ ε(t k , θ)/∂θ in a computationally efficient and accurate way using the MTF algorithm.
To facilitate application of the MTF algorithm during a gradient-based search, the candidate model should be convergent in all iterations, i.e., θ ∈ Θ.In literature, there are many gradient-based methods 1 that deal with constrained optimization problems (Papalambros & Wilde, 2000).Alternatively, one can introduce an exterior penalty function ψ(θ), which turns the constrained optimization problem into an unconstrained one, where The modified cost function then reads as follows: 1 MATLAB provides built-in solvers such as fmincon that can handle constrained optimization problems.
with J(θ ) defined in (5).It can be minimized by any gradientbased optimization algorithm (such as the Levenberg-Marquardt algorithm) to yield θN,mod = arg min θ J mod (θ ). (9) The modified cost function J mod (θ ) is unbounded for models not inside the set of the convergent models, but leaves the cost function J(θ ) unaffected for models inside the set of convergent models.Consequently, as long as the set of convergent models is respected, the gradient of the cost function with respect to the model parameters remains unaffected.In contrast to interior penalty functions (also known as barrier functions), the exterior penalty function does not result in a bias in the parameters, i.e., θN in (6) equals θN,mod in (9).

Computation of steady-state responses & gradient information
As shown in the previous section, the steady-state output z(t, θ) of model ( 2) is required at sampling instances t k in the evaluation of the cost function ( 5) and the computation of its gradient (7).First, we introduce the MTF algorithm, which enables fast computation of the steady-state response of convergent Lur'e-type models.After that, we present a method to compute the gradient ∂ ε(t, θ)/∂θ , required in (7), again, in a computationally efficient manner using the MTF algorithm.

Mixed-time-frequency algorithm
To overcome the time-consuming drawbacks of computing model responses using numerical forward integration, Pavlov et al. ( 2013) developed the so-called MTF algorithm.If the underlying Lur'e-type model satisfies the conditions stated in Theorem 1, then this algorithm computes the steady-state model response z(t, θ) efficiently under a periodic excitation w(t).For notational convenience, the dependency on θ is dropped in this section.
The MTF algorithm is an iterative algorithm.In each iteration i, two mappings are involved: where the nonlinear steady-state operator F uy • y i := −ϕ(y i , w), and the linear steady-state operators F yu , F yw map periodic signals u i+1 and w to the periodic steady-state output y i+1 .The MTF algorithm relies on the fact that the composed operator F uy •F yu is a contraction operator acting from the space 2 L 2 (T ) to L 2 (T ), if the model satisfies the conditions of Theorem 1.For computational efficiency, the nonlinear mapping F uy • y i is evaluated in time domain, while the linear mappings F yu and F yw are evaluated in frequency domain according to the frequency response functions (FRFs): Ȳ and its time-domain representation is denoted by ȳM (t).
Having the steady-state ȳ(t), one can trivially compute ū(t) in the time domain using F uy and the output z(t) in the frequency domain using the following FRFs: The following theorem shows that the MTF algorithm converges to the unique 'true' steady-state model response and gives an accuracy bound for the obtained steady-state solution.
Theorem 5 (Pavlov et al., 2013).Under the conditions of Theorem 1, for any M > 0, there is a unique limit ȳM to the sequence y i , i = 1, 2, . .., resulting from the iterative process with truncation (10a), (10b).Moreover, Theorem 5 shows that can be made arbitrarily small.Hereto, observe that, for large M, γ M yu drops to zero (thanks to the transfer from u to y being strictly proper) and that drops to zero by the Riesz-Fischer theorem, see, e.g., Beals (2004).From here, also can be made arbitrarily small.Consider with F zu and F zw being linear steady-state operators that map the inputs u and w to the steady-state output z.By the Riesz-Fischer theorem, it is clear that both converge to zero for large M, which guarantees that also converges to zero for large M.
The MTF algorithm is presented in Algorithm 1.In lines 1 to 4, the contribution of the excitation signal w(t) in the steady-state output ȳ(t) is computed in the frequency domain and transformed to the time domain.After that, inside the while loop, first the nonlinearity is evaluated in the time domain in line 6 and its output is transformed to the frequency domain in line 7. Next, in lines 8 and 9, the LTI dynamics are evaluated in the frequency domain and the output y(t) is transformed to the time domain.
After convergence of the signal y(t), the model output z(t) is computed in the frequency domain in line 12 and transformed to the time domain in line 13.The termination criterion is the normalized mismatch between the Fourier coefficient Y i and those of the previous iteration Y i−1 , measured in the ℓ 2 signal norm, which is defined as The MTF algorithm with truncation can be considered as a multiharmonic variant of the describing function method, see Mees (1972) for a discussion for the autonomous case.The computed steady-state response z(t, θ) can be used together with the measured response z(t) to compute the output error (4), which can be subsequently used in ( 5) and ( 7) in an optimization routine to solve the underlying identification problem.

8:
Evaluate LTI dynamics in frequency domain 11: end 12: Evaluate LTI dynamics in frequency domain

Gradient computation
Any gradient-based optimization routine requires the gradient of the cost function with respect to the model parameters as in (7).It is well known that this gradient can be computed by simulation of a parameter sensitivity model (Keesman, 2011;Suykens, Demoor, & Vandewalle, 1995).
The sensitivity model is a continuous-time model with a similar structure as the proposed model structure.Therefore, in other works, its response is computed in the same way as the model response is computed, which is typically done by numerical forward integration.We show that the sensitivity model in our identification problem is again a convergent Lur'e-type model satisfying the conditions of Theorem 1, which facilitates the use of the efficient MTF algorithm to compute its steady-state response, i.e., the gradient information of the steady-state model response with respect to the parameters.The methodology extends that in Pavlov et al. (2013), where the to-be-optimized parameters only appeared in the static nonlinear block, whereas in the more generic case in (2), the parameters appear in both the LTI and nonlinear blocks.The theorem below presents the parameter sensitivity model, where x(t, θ), ū(t, θ) and ȳ(t, θ) of (2) enter as inputs.Partial derivatives with respect to θ are denoted by subscript θ, e.g., x θ := ∂x/∂θ .Theorem 6.Consider model (2).Under the conditions of Theorem 1, if the partial derivative ∂ϕ(y, w, θ )/∂y exists and is continuous for y ∈ R, θ ∈ Θ, and if the partial derivatives ϕ θ (y, w, θ exist and are continuous for y ∈ R, θ ∈ Θ, then the partial derivative ∂ ε(t, θ)/∂θ i is the unique T -periodic steady-state output zθ i (t, θ) of the following model: Furthermore, model (12) satisfies the conditions of Theorem 1 with A, B, C, φ(y, w) replaced respectively by A(θ ), B(θ ), C (θ ), ϕ y (ȳ(t, θ), w(t), θ) y θ i (t, θ) and model ( 12) is convergent according to Definition 1.

Proof. The proof can be found in Appendix B. ■
The above theorem shows that the gradient ∂ ε(t, θ)/∂θ i , for i = 1, . . ., n θ , is the steady-state response of a convergent Lur'etype model satisfying the conditions of Theorem 1.Hence, its steady-state response can be computed using the MTF algorithm for every parameter θ i in θ individually.We emphasize that the gradient information can be computed with arbitrary accuracy using the MTF algorithm, which can positively affect the convergence speed of the gradient-based search used to solve the identification problem in (6).
Remark 7. Once the steady-state output z(t, θ) of ( 2) is computed by the MTF algorithm, the steady-state state x(t, θ), required in (13), is computed in the frequency-domain from the linear part of (2).

True system & model class
To illustrate the effectiveness of the identification approach, we consider a simulation example on the double mass-springdamper system depicted schematically in Fig. 2. The dynamics are described by the Lur'e-type system in (1) with the matrices defined in (15) and the following nonlinearity: The static nonlinearity has the physical interpretation of a nonlinear spring between the fixed earth and the first cart.The system is characterized by masses m 1 = 0.15 kg, m 2 = 0.45 kg, linear damping constants d 1 = d 2 = 0.4 Nm/s, linear spring constants k 1 = 1000 N/m, k 2 = 2200 N/m, and parameters defining the nonlinear spring φ(y, w) in ( 14): The excitation w is the force exerted on the second cart and the output z is the position of the second cart.The Bode magnitude diagram of the transfer functions of the LTI part of the true system and the graph of the nonlinearity are included in Fig. 3.
The goal is to identify a black-box Lur'e-type model for this system.Hereto, consider the model class in (2) with state dimension n = 4 and consider a full parametrization for the LTI part, ϕ(y, w, θ ) , where each element of the vectors is a model parameter.This nonlinear function represents a specific type of a single hidden-layer feedforward neural network with hyperbolic tangent activation functions and a linear output layer (Suykens et al., 1995).This model is characterized by n θ = 44 parameters, which we collect in the vector θ ∈ R n θ .It can be shown that the conditions of Theorem 1 are satisfied for any where λ(A(θ)) denotes the eigenvalues of the matrix A(θ ), C <0 denotes the set of complex numbers with negative real part, R ≥0 denotes the set of nonnegative real numbers and K (θ ) := 1 2 W (θ ) [2] ⊤ W (θ ) [1] .For any θ ∈ Θ, the proposed model is globally exponentially convergent according to Definition 1.The true system parametrized by θ 0 is in the model class, i.e., θ 0 ∈ Θ.

Identification results
The excitation signal w is a random-phase multisine that excites frequencies between 0.5 Hz and 200 Hz with 0.5 Hz spacing and is sampled at 500 Hz, implying a period time of T = 2 s.
Following Assumption 2, discrete-time white noise is added to the steady-state output z0 , where the variance σ 2 e is selected such that a signal-to-noise (SNR) of 20, 40, 60 and ∞ dB is realized.
The SNR of ∞ corresponds to a deterministic setting.To find an initial parameter set θ init , we employ the so-called controlled random search (CRS) (Hendrix, Ortigosa, & García, 2001), which is a global, non-gradient-based optimization routine.In accordance to Remark 4, we terminate the CRS prematurely.In the second identification step, we perform a gradient-based search using Matlab's fmincon implementation of (Byrd, Gilbert, & Nocedal, 2000) to solve the constrained optimization problem (6).
The results are included in Table 1.Here, it can be seen that the cost function is successfully minimized to J( θN ) ≈ σ 2 e .For the specific case of a SNR of 60 dB, Fig. 3 presents the Bode magnitude diagram of the LTI part and the graph of the nonlinearity of the true system (θ 0 ), initial model (θ init ), and the identified model ( θN ).This figure evidences an accurate match between the identified model and the true system.For the 20 dB SNR case, the time-domain response is depicted in Fig. 4. In the top plot, it can be observed that the steady-state error ε(t, θN ) is significantly smaller than ε(t, θ init ).The bottom plot displays a validation test, where a new realization of the excitation signal (15) Fig. 3. Bode magnitude diagrams of the involved transfer functions and the graph of the nonlinearity of the true system (θ 0 ), initial model (θ init ) and identified model ( θN ).For visualization, the gain exchange between the LTI block and nonlinearity is fixed by normalization of W [1] , W [2] and a loop transformation is performed (Khalil, 1996).
w is used that is 10 times larger in amplitude than the w that is used in identification.This validation test shows that the model accurately describes the steady-state response of the true system.We would like to note that all identified models are guaranteed globally exponentially converging and, therefore, these identified models can safely be used for other excitation signals than the one used in identification, without showing instability issues.

Numerical efficiency
To illustrate the numerical efficiency 3 of our method, we consider the deterministic scenario with an SNR of ∞, i.e., σ 2 e = 0 in Assumption 2 The identification strategy turns out to yield an 'almost' perfect model for the true system evidenced by the cost J( θN ) being minimized up to Matlab's numerical 16-digit precision.In both identification steps, we compute the steadystate model responses z and the gradient zθ by two methods, namely by the MTF algorithm presented in Section 5.1 and by numerical forward integration (NFI).To guarantee a similar level of accuracy as in MTF, we simulate NFI for 20 periods and take the last period as steady-state.The result is presented in Table 2.It can be observed that in the first identification step, where no gradient information is required, the MTF algorithm reduces the computation time from over 1 and a half hour to 32 s.In the second identification step, where also gradients have to be computed, the MTF algorithm reduces the computation time from over 9 and an half hours to only 3 min and 16 s.In  total, the computation time is reduced with an factor of approximately 150.Such small computation time enables application of global optimization routines for initialization and also enables an  efficient and effective gradient-based search.Besides highlighting the numerical benefits of our approach, this example also illustrates one of reasons that literature around identification of nonlinear models is mainly focused on the discrete-time setting; namely, computing steady-state model responses and gradient information is computationally prohibitively expensive in the continuous-time case using NFI.

Experimental case study: Mechanical ventilation
The identification strategy developed in this paper is applied on a mechanical ventilation setup as depicted in Fig. 5. Mechanical ventilation is used in intensive-care units to assist or stimulate respiration of patients who are unable to breathe on their own.The treatment quality of mechanical ventilation depends on two aspects.The first aspect is the design of the breathing pressure profile that is fed to the patient (Amato et al., 2015;Lachmann, 1992).The second aspect is the ability to track the designed breathing pressure profile.Hereto, in the literature, pressure control strategies are proposed to ensure comfortable and stable air flow (Borrello, 2005;Hunnekens et al., 2018;Reinders et al., 2020;Tehrani et al., 2004).In both these aspects, knowledge on the health condition of the patient's lungs is crucial as incorrect treatment can lead to lung damage and, eventually, to patient mortality (Gajic et al., 2004;Slutsky & Ranieri, 2013).Therefore, by identification of a first-principle model, we aim to provide the patient's health information reflected by the patient's lung parameters in a fast way.

First-principle modeling & convergent dynamics
Fig. 6 gives a schematic view of the setup.Using conservation of flow and taking into account a (signed) quadratic hose resistance, the patient-hose dynamics can be written as follows (Heck, 2018;van de Kamp, 2019): where the scalar nonlinearity is given by ϕ(y, θ) The excitation signal w is the pressure p blower generated by the blower, the output z is the airway pressure p airway , the state x represents the lung pressure p lung , the output of the nonlinearity u is a scaled flow Q hose through the hose, and, finally, the input y of the nonlinearity is a scaled pressure drop over the hose p airway −p blower .The intended leakage flow Q leak flushes the system from CO 2 -rich exhaled air.The model is parametrized by the lung resistance θ 1 in mbar sec/ml, lung compliance θ 2 in ml/mbar, the linear hose resistance θ 3 in mbar sec/ml, the quadratic hose resistance θ 4 in mbar sec 2 /ml 2 and the linear leakage resistance θ 5 in mbar sec/ml.For any α > 0, the following nonlinear scaling of the parameters leaves the steady-state response z of model ( 17) unchanged for any bounded excitation w (van de Kamp, 2019): This scaling corresponds to a gain exchange between the LTI part of the model and the nonlinearity ϕ, see Remark 3. Consequently, parameters θ 1 , . . ., θ 5 cannot be uniquely identified on the basis of w and z only.As a remedy, we consider the parameter θ 5 to be known a priori and positive in value, since it is only leakage specific, i.e., independent of the patient and the hose, and it can be found relatively easily through calibration.Given knowledge on θ 5 , the to-be-identified parameters are θ 1 , . . ., θ 4 , which are collected in θ.It is shown in van de Kamp (2019) that model ( 17) satisfies all conditions of Theorem 1 for θ i > 0, for i = 1, . . ., 4, yielding the set Therefore, for any θ ∈ Θ, model ( 17) is globally exponentially convergent according to Definition 1.

Identification experiment
Identification of continuous-time models from discrete data points is a challenging task due to the unobserved intersample behavior of the true system.Analysis tools for errors due to intersample behavior and practical guidelines for the design of the excitation signal are presented (Goodwin, Aguero, Cea Garridos, Salgado, & Yuz, 2013;Rao & Unbehauen, 2006;Schoukens, Relan, & Schoukens, 2017;Soderstrom, Fan, Carlsson, & Bigi, 1997).As we are dealing with a medical application, there is no freedom in the design of the excitation signal; it is the signal that is also used during conventional operation.We can only set a standard target pressure profile p target for the blower, which is tracked by the blower (using an internal control loop) to realize the pressure  p blower at the blower side of the hose, see Fig. 7.The shape of the blower pressure represents a breathing cycle that starts with an inhalation phase and follows by an exhalation phase.The measured pressure p blower is used as the excitation signal w(t) for identification, implying that the internal control loop generating p blower from p target is not part of the to-be-identified dynamics.The period time of p target is T = 4 s, with a minimum pressure of 5 mbar and a maximum pressure of 30 mbar.Rather than using humans in these experiments, the ASL5000 breathing simulator is utilized, which is specifically designed to emulate the lung behavior of patients.Furthermore, we consider the case where the patient is completely sedated, implying no breathing activity from the patient's side.
We apply 15 periods of the excitation signal depicted in Fig. 7 and sample the realized blower pressure p blower and the patient's airway pressure p airway uniformly at a sampling frequency of 500 Hz.The average of the last 12 periods of the measured p blower is considered as the periodic input data w(t) and the average of the last 12 periods of p airway is considered as the steady-state output data z(t), both with a period time of T = 4 s.The total experiment time is only 1 min for each experiment, which is of crucial importance in practice as time is costly in this medical application.
The three hose-patient configurations considered are listed in Table 3, where the lung parameters θ 1 and θ 2 are set on the ASL5000 unit with an accuracy 4 of 10% and 5%, respectively.The lung resistance θ 1 remains unchanged over the experiments.To cover a set of different configurations, both the lung compliance θ 2 and the leakage resistance θ 5 are changed from experiment to experiment.The same hose with parameters θ 3 and θ 4 is used for all experiments.The considered configurations represent a healthy patient for configuration 1 & 3, and a patient with less stiff lungs for configuration 2.

Identification results
The identification problem in (6) is solved using the modified cost function in (8) with the constraint θ ∈ Θ with Θ defined in (18).The set of initial parameter reflects an 'average' configuration, see ISO Central Secretary (2020), van de Kamp (2019): This initialization method employs the prior knowledge that parameters are typically close to the 'average' values.The identified parameters θN obtained by solving the optimization problem in (6) are presented in Table 4.
The measured response z(t) for configurations 1, 2 and 3 is shown in Figs.8(a)-8(c), respectively.These figures also depict the steady-state output error ε(t, θ init ) and ε(t, θN ), produced by the initial and identified model, respectively, with ε defined in (4).For comparison, also an LTI transfer-function model is identified using the TFEST Matlab routine, initialized by the N4SID routine.The parameters θ LTI corresponding to the identified LTI models have no physical meaning.The steady-state output error ε(t, θ LTI ) is also included in the respective figures.
From Figs. 8(a)-8(c), we draw the following two conclusions.Firstly, the optimization problem ( 6) is successfully solved, as the error ε(t, θN ) of the identified nonlinear model is significantly smaller than the error ε(t, θ init ) of the initial model.Secondly, a linear model is insufficient to describe the nonlinear dynamics of this system, since the error ε(t, θN ) of the identified nonlin- ear model is significantly smaller than the error ε(t, θ LTI ) of the identified linear model.Both these conclusions are confirmed by Table 4, which shows that the cost J( θN ) of the identified nonlinear model is significantly smaller than both the costs J(θ init ) and J(θ LTI ).
The model ( 17) is a simplification of the system depicted in Fig. 5.For example, calibration experiments suggest that the leakage component exhibits a quadratic pressure-flow relation, whereas this relation is linear in model ( 17).Furthermore, the static hose model in model ( 17) does not account for hose dynamics.Consequently, a part of the unmodeled dynamics is compensated for by the identified model parameters, which is reflected as follows.Firstly, although the same hose is used in all configurations, the results in Table 4 show that its associated parameters θ 3 and θ 4 are different with respect to each other and also with respect to the ones in Table 3.However, the graph of the identified nonlinearities in each configuration in Fig. 9 shows that all identified models represent the same hose characteristics in the domain of interest.Besides that, the values for θ 3 and θ 4 in Table 3 are obtained by a static calibration measurement, which is not fully representative in case the hose is used in a dynamic identification experiment.Secondly, the identified parameters in Table 4 show that θ 1 and θ 2 are estimated slightly too large compared to the preset values as in Table 3.We are informed by medical personnel that an accuracy of 15% is sufficient for the lung parameters in the scope of the decision-making process of patient treatment.If the mean of the uncertain θ 1 and θ 2 in Table 3 are the true values, then we slightly exceed this accuracy requirement.We would like to note that a bias in model parameters is not a deficiency of the identification approach, but rather a matter of unmodeled dynamics.

Numerical efficiency
To illustrate the computational advantages of the proposed approach, the total computation time required to solve the optimization problem (6) is included in Table 4 in the column MTF.This table also includes the total computation time required to solve the same optimization problem with numerical forward integration (NFI).To guarantee a similar accuracy level for NFI as for MTF, ten periods of the excitation signal w(t) are applied to the model in the NFI simulations and the last period of z(t, θ) is taken as the steady-state model output z(t, θ).Table 4 shows that the total computation time when using the MTF algorithm is reduced from roughly 2 and a half minutes to only 14 s.The quick availability of the knowledge on the model parameters enables a high treatment quality that can start early.

Table 4
Identification results of the mechanical ventilation application.Column 1 represents the configuration number.Columns 2-5 contain the identified parameters θN after solving the identification problem (6).Columns 6-8 contain the value of the cost function (5) for the identified black-box LTI model (θ LTI ), initial model (θ init ) in ( 19) and identified model ( θN ).Columns 9-10 contain the elapsed computation time in seconds that took to solve the identification problem (6).

Conclusions
This paper presents an identification approach for convergent continuous-time Lur'e-type systems.The benefits of the proposed approach are that (i) it guarantees that the identified model preserves the convergence property, (ii) it is computationally attractive, and, (iii) it is applicable to a large class of block-oriented feedback systems.In a simulation example, we demonstrated the effectiveness and computational efficiency of our identification approach.Furthermore, in an experimental study on mechanical ventilation in hospitals, the identification approach has been shown to be effective to identify the parameters of a firstprinciple model in a fast way, which enables improved patient treatment.
In the final step, having (A.1) and (A.3) at hand, along the same lines as the proof of Theorem 2 in Söderström (1974), it is concluded that the parameter vector θN converges in probability to θ 0 as N → ∞, i.e., lim This completes the proof.

Appendix B. Proof of Theorem 6
The proof is given for the scalar case of θ and can be repeated analogously for each component of a vector-valued θ.First, consider the following property of Lur'e-type models (2) satisfying the conditions of Theorem 1.
Let us next prove that model ( 12) satisfies conditions C1-C3 of Theorem 1. Condition C1 is satisfied since the matrices A, B and C of the LTI block of ( 12) are the same as those of model (2).Condition C2 holds with the same K as for the model ( 2
The MTF algorithm iteratively evaluates the mappings (10b), (10a), while transforming the intermediate signals between the time-and frequency-domain using the (Inverse) Fast Fourier Transform ((I)FFT), truncated to only M frequency contributions.Hereto, denote by Ȳ [m] the sequence of Fourier coefficients of 2 The space L 2 (T ) denotes the space of piecewise-continuous real-valued T -periodic scalar functions y(t) satisfying ∥y∥ L 2 < ∞, where ∥y∥ 2 t)| 2 dt.ȳ(t).Then, the truncated version of Ȳ [m], denoted by Ȳ M[m] is defined as follows:

3
The computations of the examples in Sections 6 and 7 are carried out on an Intel Core i7-7700HQ, 2.8 GHz processor.

Fig. 4 .
Fig. 4. Measured z and the errors ε as defined in (4).Top: results for the identification excitation signal.Bottom: results for the validation excitation signal.

Table 1
Identification results for different SNRs.

Table 2
Comparison of computation time between MTF and NFI for both identification steps.The number of computed steady-state model responses is included in the second column.

Table 3
Considered experimental cases with uncertain parameters θ 1 and θ 2 .
Measured z and the errors ε as defined in (4).9.Identified nonlinearity ϕ(y) for each configuration in Table4over the domain of interest in y.