A sliding-window approach for latent restoring force modeling

A sliding-windo w approac h f or latent res tor ing f orce modeling ⋆ Merijn Floren a,b , ∗ , Jan Swe vers a,b a KU Leuven, Department of Mec hanical Engineering, Leuven, 3001, Belgium b Flanders Make@KU Leuv en, MPRO Core Lab, Leuv en, 3001, Belgium A R T I C L E I N F O Keyw ords : Nonlinear system identification Restoring f orce surface Random-phase multisine ex citation State-space models Best linear appro ximation A B S T R A C T Restoring f orce surface (RFS) methods offer an attractiv e nonparametr ic framew ork for iden- tifying nonlinear restoring f orces directly from data, but their reliance on complete kinematic measurements at each degree of freedom limits scalability to multidimensional systems. The aim of t his paper is to o vercome t hese measurement limitations by pr oposing an identification framew ork with relaxed sensing requirements that exploits periodic multisine excitation. Starting from an initial linear model, a sliding-window feedbac k approach reconstructs latent states and nonlinear restoring f orces nonparametrically , enabling identification of the nonlinear component through linear-in-par ameters regression instead of highly non-conve x optimization. V alidation on synthetic and experimental datasets demonstrates high simulation accuracy and reliable recov ery of physical parameters under par tial sensing and noisy conditions. 1. Introduction Identifying state-space models of mechanical systems from input-output data is inherently challenging, as the system dynamics depend on latent states that are not directl y measurable. As a result, parameter estimation typically requires solving a highly nonlinear and non-con ve x optimization problem with an intrinsic recurrent str ucture imposed by the s tate e volution. How ev er , this difficulty can be significantly alleviated if an estimate of the latent state trajector y can be obtained separately from the parameter estimation step. In t hat case, the recur rent nature of the identification problem effectiv ely disappears, reducing parameter es timation to a considerably simpler regression problem. The abo ve procedur e illustrates the fundamental rationale underl ying restoring f orce surface (RFS) me thods. As a concrete example, consider the single degree of freedom (SDOF) sy stem studied in the original RFS formulation [ 12 ]: 𝑚  𝑦 ( 𝑡 ) + 𝑓  𝑦 ( 𝑡 ) ,  𝑦 ( 𝑡 )  = 𝑢 ( 𝑡 ) , (1) where 𝑚 denotes the mass, 𝑦 ( 𝑡 ) the displacement, 𝑢 ( 𝑡 ) t he external ex citation, and 𝑓 ( ⋅ ) the total restoring force, representing the combined linear and nonlinear stiffness and damping forces acting on the system. In the or iginal w ork [ 12 ], it is assumed that both the ex citation 𝑢 ( 𝑡 ) and t he acceleration  𝑦 ( 𝑡 ) are measured (or other wise a vailable), and that the mass 𝑚 is known or can be reliably estimated. Under these assumptions, the restoring f orce can be isolated through simple algebraic rearrangement: 𝑓  𝑦 ( 𝑡 ) ,  𝑦 ( 𝑡 )  = 𝑢 ( 𝑡 ) − 𝑚  𝑦 ( 𝑡 ) . (2) The r ight-hand side of ( 2 ) is theref ore directly comput able, enabling a fully nonparametric reconstruction of the restoring force behavior . The cor responding arguments of 𝑓 ( ⋅ ) are obt ained from the measured acceleration  𝑦 ( 𝑡 ) through numerical integration, effectivel y r econstr ucting the latent states needed to av oid recursive parameter estimation. With these quantities a vailable, the restoring f orce can be appro ximated using any suitable nonlinear function approximator 1 . Crucially , this modeling step is per formed in a static reg ression setting, where t he restoring force is treated as an explicit algebraic mapping from 𝑦 ( 𝑡 ) and  𝑦 ( 𝑡 ) to the computed f orce values, i.e., t he right-hand side of ( 2 ). As a result, identifying 𝑓 ( ⋅ ) reduces to a supervised function approximation problem wit h independent regression samples, rather than a recursive and non-con vex dynamical es timation problem. ⋆ This wor k is suppor ted by Flanders Make’ s IR V A projects ASSISStaNT and CoMoDO. ∗ Corresponding aut hor , e-mail: merijn.floren@kuleuven.be . OR CI D (s): 0000-0001-7265-7851 (M. Floren); 0000-0003-2034-5519 (J. Swev ers) 1 The original study [ 12 ] employ ed polynomial basis functions, while later work also considered neural-netw ork representations [ 13 ]. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier Page 1 of 20 Although first introduced in 1979 [ 12 ], RFS methods remain relev ant today; the survey in [ 15 ] recognizes them as a major time-domain identification techniq ue, valued for their simplicity and intuitive, visual inter pretations. RFS methods hav e, f or ex ample, been used to characterize t he complex nonlinear forces at the wing-to-pa yload mounting interface in an F-16 fighter jet [ 3 ]. Other applications include predicting the response of elastomer materials [ 19 ], analyzing the variable stiffness of an elastomagnetic suspension [ 1 ], and identifying nonlinear behavior in a robotic arm [ 6 ]. 1.1. Limitations of classical RFS methods Despite its success, the or iginal RFS method exhibits se veral limitations. First, it relies on kno w ledge of acceleration, velocity , and displacement. Ideally , all three quantities should be measured directly; how ev er, this is often impractical or prohibitiv ely e xpensiv e. In most applications, only one quantity is measured, while the remaining vari- ables are obt ained through numer ical processing, a procedure that can introduce significant errors [ 22 ]. In mechanical systems, acceleration is typically measured, with velocity and displacement obtained via numerical integ ration. This procedure tends to amplify low -frequency noise and sensor drif t o ver time, often requiring additional post-processing. Con versel y , if displacement is measured, v elocity and acceleration must be obtained through numerical differ entiation, which rapidl y becomes unreliable in t he presence of measurement noise. A detailed discussion of the implications of numerical integ ration and differentiation f or RFS-based approaches can be f ound in [ 22 ]. A second limitation of RFS methods arises from their demanding measurement requirements in higher -dimensional systems. Specifically , extending the SDOF restor ing-f orce isolation in ( 2 ) t o t he multi degree of fr eedom (MDOF) case requires measurements at ev er y degree of freedom [ 11 , 14 ]. In practice, for high-dimensional systems, instrumenting all deg rees of freedom quickl y becomes impractical, as sensors may alter system dynamics, compromise str uctural integr ity , be physicall y inaccessible, or result in ex cessive costs. Although numer ical post-processing and increased sensor cov erage can par tially alleviate t hese limitations, they do not fundament all y remo ve the str ict measurement dependence of RFS methods. Recent wor k has therefore sought to address t his issue more directly . In [ 18 ], a Ba yesian f or mulation is proposed in which a linear system of ordinar y differential equations (ODEs) is driven by a Gaussian process in time representing the unknown nonlinear component of the restoring force, enabling joint infer ence of both the latent state trajectory and the nonlinear restoring force. A related deterministic method was introduced in [ 5 ], where the nonlinear restoring f orce time ser ies is analyticall y reconstr ucted from an initially fitted linear model using a sliding-window procedure. Both approaches, howe ver , explicitly assume an SDOF system f or mulation and therefore do no t extend directl y to the MDOF setting. 1.2. Proposed approach This w ork extends the rationale introduced in [ 5 ] to the MDOF se tting, thereb y substantiall y relaxing the measurement requirements of RFS methods in higher-dimensional systems. In par ticular , no restrictions are imposed on the type of measured quantity (displacement, velocity , or acceleration), nor on the number or spatial locations of the measured degrees of freedom 2 . Follo wing [ 5 ], the only requirement is t hat the system is e xcited using random-phase multisine signals, which enable estimation of the best linear appro ximation (BLA) from frequency response measurements [ 16 ]. This initial linear model is subsequentl y employ ed within a sliding-window framew ork to analyticall y reconstruct both t he nonlinear restoring force and the associated latent state trajectory in t he time domain. While the spatial location of the nonlinear restoring f orce is assumed to be known, its functional f orm and magnitude remain unknown. The linear dynamics ar e assumed to f ollow a known ODE structure, whereas t he phy sical parameters of this linear subsystem are only appro ximatel y known through a vailable initial estimates. The objective is twof old: to recov er the tr ue phy sical parameters of the underlying linear sy stem and to obtain an accurate nonlinear model suitable f or simulation and control applications. 1.3. Notation The sets of real, integer , and natural numbers are denoted by ℝ , ℤ , and ℕ , respectivel y . The identity matr ix is denoted b y 𝐼 , and 0 denotes the zer o matrix, with dimensions clear from context. For a real-v alued v ector 𝑥 ∈ ℝ 𝑛 and a symmetric positive definite matr ix 𝑄 ∈ ℝ 𝑛 × 𝑛 , the squared weighted 2-norm is defined as  𝑥  2 𝑄 = 𝑥 T 𝑄𝑥 , where 𝑥 T denotes t he transpose of 𝑥 . When a unit w eighting is used, the notation simplifies to  𝑥  2 = 𝑥 T 𝑥 . For complex-v alued vect ors, the same definition applies ex cept that the Her mitian transpose, denoted by 𝑥 H , replaces the regular transpose. The imaginary unit is denoted by 𝑗 . The exponential of a square matrix 𝐴 is wr itten as e xp( 𝐴 ) . 2 Provided that the measured responses are dynamically coupled to the applied ex citation. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 2 of 20 Figure 1: The nonlinear linea r fractional representation (NL-LFR) model structure as a feedback interconnection of a linear time-inva riant (L TI) system and a static nonlinear mapping. 1.4. Paper outline The remainder of this paper is str uctured as follo ws. Section 2 defines the problem s tatement, after which Section 3 details the proposed sequential identification procedure. The methodology is validated using three case studies: a simulated SDOF system in Section 4.1 , a simulated MDOF system wit h unmeasured nonlinear locations in Section 4.2 , and experimental validation on an SDOF benc hmark setup in Section 5 . Conclusions are dra wn in Section 6 . 2. Problem statement W e consider the per iodic in, same period out (PISPO) system class, compr ising systems f or which a per iodic input signal yields a steady-s tate output with the same per iod 3 . The governing ODEs are assumed to f ollow a known structural f or m, wit h the ex ception of an unknown nonlinear restoring force, and are represented within the nonlinear linear fractional repr esentation (NL -LFR) framewor k, in which nonlinear systems are described as a f eedback interconnection between an linear time-in variant (L TI) dynamical system and a static nonlinear mapping, as illustrated in Fig. 1 . The cor responding continuous-time system dynamics are descr ibed by t he follo wing state-space equations:  𝑥 ( 𝑡 ) =  𝑥 ( 𝑡 ) +  𝑢 𝑢 ( 𝑡 ) +  𝑤 𝑤 ( 𝑡 ) , (3a) 𝑦 0 ( 𝑡 ) =  𝑦 𝑥 ( 𝑡 ) +  𝑦𝑢 𝑢 ( 𝑡 ) +  𝑦𝑤 𝑤 ( 𝑡 ) , (3b) 𝑧 ( 𝑡 ) =  𝑧 𝑥 ( 𝑡 ) , (3c) 𝑤 ( 𝑡 ) = 𝑓  𝑧 ( 𝑡 )  , (3d) where 𝑥 ( 𝑡 ) ∈ ℝ 𝑛 𝑥 is the latent state, 𝑢 ( 𝑡 ) ∈ ℝ 𝑛 𝑢 the external input, 𝑤 ( 𝑡 ) ∈ ℝ 𝑛 𝑤 the nonlinear feedbac k input 4 , 𝑦 0 ( 𝑡 ) ∈ ℝ 𝑛 𝑦 the exact output, and 𝑧 ( 𝑡 ) ∈ ℝ 𝑛 𝑧 the latent argument of t he nonlinear mapping 𝑓 ∶ ℝ 𝑛 𝑧 → ℝ 𝑛 𝑤 . The matrices  ∈ ℝ 𝑛 𝑥 × 𝑛 𝑥 ,  𝑢 ∈ ℝ 𝑛 𝑥 × 𝑛 𝑢 ,  𝑤 ∈ ℝ 𝑛 𝑥 × 𝑛 𝑤 ,  𝑦 ∈ ℝ 𝑛 𝑦 × 𝑛 𝑥 ,  𝑧 ∈ ℝ 𝑛 𝑧 × 𝑛 𝑥 ,  𝑦𝑢 ∈ ℝ 𝑛 𝑦 × 𝑛 𝑢 , and  𝑦𝑤 ∈ ℝ 𝑛 𝑦 × 𝑛 𝑤 describe the linear dynamics and input-output relationships. Since the linear system of ODEs is kno wn, the str ucture of t hese matr ices is also kno wn. Ho we v er, the v alues of the corresponding ph ysical parameters 𝜃 phys remain unknown, although an initial guess is a vailable, denoted as 𝜃 phys , 0 . The tr ue parameters 𝜃 phys cor respond to the underlying linear syst em , which ex clusivel y captures all linear input-output dynamics. Accordingly , the nonlinear mapping 𝑓 ( ⋅ ) in ( 3d ) is purely nonlinear 5 . Data are collected from ( 3 ) by e xciting the sy stem with random-phase multisine inputs [ 16 ]: 𝑢 ( 𝑛 ) ∶ = 2  𝑁 𝑁 ∕2−1  𝑘 =1 𝑈 𝑘 cos(2 𝜋 𝑘𝑓 0 𝑛 + 𝜑 𝑘 ) , (4) where 𝑈 𝑘 denotes the ex citation amplitude at frequency index 𝑘 , and 𝑓 0 = 𝑓 𝑠 ∕ 𝑁 is the frequency resolution, with 𝑓 𝑠 denoting the sampling frequency and 𝑁 the (e ven) number of samples per period. The cor responding frequency of the 𝑘 th harmonic is giv en by 𝑓 𝑘 = 𝑘𝑓 0 . The phases 𝜑 𝑘 are independent and uniforml y dis tr ibuted ov er the inter v al [0 , 2 𝜋 ) . 3 The PISPO class encompasses a broad rang e of systems, including those exhibiting amplitude-dependent resonance and nonlinearities such as saturation and dead zones, while ex cluding phenomena that produce subharmonics, such as bifurcations and chaotic behavior . 4 That is, the nonlinear restoring force. 5 For instance, 𝑓 ( 𝑧 ) cannot be expressed as 𝑓 ( 𝑧 ) = 𝐾 𝑧 + 𝑔 ( 𝑧 ) , where 𝐾 is a constant matrix and 𝑔 ( ⋅ ) is a nonlinear function. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 3 of 20 Assumption 1. The input signal 𝑢 ( 𝑛 ) is persistentl y ex citing of sufficiently high order, such that all relev ant system dynamics are ex cited over the freq uency band of interest. Assumption 2. The discrete-time signal ( 4 ) is applied to the continuous-time sys tem ( 3 ) via zero-order hold (ZOH), i.e., 𝑢 ZOH ( 𝑡 ) ∶ = 𝑢 ( 𝑛 ) , f or 𝑡 ∈ [ 𝑛𝑇 𝑠 , ( 𝑛 + 1) 𝑇 𝑠 ) , where 𝑇 𝑠 = 1∕ 𝑓 𝑠 denotes the sampling period. Assumption 3. The measured output is cor r upted by additive, zero-mean, stationar y noise 𝑣 ( 𝑛 ) wit h finite variance, i.e., 𝑦 ( 𝑛 ) ∶ = 𝑦 0 ( 𝑛 ) + 𝑣 ( 𝑛 ) . The noise ma y be colored and is uncor related with the input 𝑢 ( 𝑛 ) , which is ex actly kno wn. The resulting input-output dataset for parameter estimation consists of 𝑃 ≥ 1 steady-state periods of 𝑅 ≥ 𝑛 𝑢 realizations of a random-phase multisine ( 4 ):  =  𝑢 [ 𝑟,𝑝 ] ( 𝑛 ) , 𝑦 [ 𝑟,𝑝 ] ( 𝑛 )  𝑁 −1 ,𝑅 −1 ,𝑃 −1 𝑛 =0 ,𝑟 =0 ,𝑝 =0 , (5) where both the input and output signals are zero-mean. Based on  , the objective is to identify a discrete-time NL -LFR state-space model of the form: 𝑥 ( 𝑛 + 1) = 𝐴𝑥 ( 𝑛 ) + 𝐵 𝑢 𝑢 ( 𝑛 ) + 𝐵 𝑤 𝑤 ( 𝑛 ) , (6a) 𝑦 ( 𝑛 ) = 𝐶 𝑦 𝑥 ( 𝑛 ) + 𝐷 𝑦𝑢 𝑢 ( 𝑛 ) + 𝐷 𝑦𝑤 𝑤 ( 𝑛 ) , (6b) 𝑧 ( 𝑛 ) = 𝐶 𝑧 𝑥 ( 𝑛 ) , (6c) 𝑤 ( 𝑛 ) = 𝛽 T 𝜙  𝑧 ( 𝑛 )  , (6d) where t he discrete-time matrices in ( 6a ) are related to their continuous-time counter par ts in ( 3a ) through the ZOH discretizations: 𝐴 = exp(  𝑇 𝑠 ) , 𝐵 𝑢 = ∫ 𝑇 𝑠 0 exp(  𝜏 )  𝑢 𝑑 𝜏 , 𝐵 𝑤 = ∫ 𝑇 𝑠 0 exp(  𝜏 )  𝑤 𝑑 𝜏 , (7) while 𝐶 𝑦 =  𝑦 , 𝐶 𝑧 =  𝑧 , 𝐷 𝑦𝑢 =  𝑦𝑢 , and 𝐷 𝑦𝑤 =  𝑦𝑤 . Moreov er , 𝜙 ∶ ℝ 𝑛 𝑧 → ℝ 𝑛 𝜙 represents a polynomial feature mapping with 𝛽 ∈ ℝ 𝑛 𝜙 × 𝑛 𝑤 the cor responding coefficient matrix. Optimization is thus car ried out ov er the decision variables 𝜃 = ( 𝜃 phys , 𝛽 ) , which preser ves the physical inter pretability of t he underl ying linear sys tem, despite the final model being expressed in discrete time. An implication of t his c hoice is that the discrete-time matrices hav e to be recomputed at each iteration of the optimization routine that updates 𝜃 phys . For tunatel y , this can be done efficientl y using a single matrix exponential [ 20 ]:  𝐴 𝐵 0 𝐼  = exp    0 0  𝑇 𝑠  , (8) where  = [  𝑢  𝑤 ] and 𝐵 = [ 𝐵 𝑢 𝐵 𝑤 ] . In the remainder of t his paper , all discrete-time matr ices are obtained e xplicitl y via this formulation. Remar k 1. The ZOH discretization method in ( 7 ) assumes that both inputs 𝑢 ( 𝑡 ) and 𝑤 ( 𝑡 ) are piecewise constant ov er each sampling inter v al 𝑡 ∈ [ 𝑛𝑇 𝑠 , ( 𝑛 + 1) 𝑇 𝑠 ) f or all 𝑛 . Under ideal experimental conditions, this assumption holds f or the e xter nal e xcitation 𝑢 ( 𝑡 ) , but not for the inter nal feedbac k signal 𝑤 ( 𝑡 ) , which varies continuously . Consequently , the sampled continuous-time NL -LFR system in ( 3 ) cannot be represented exactl y by the discrete-time model in ( 6 ). Ne vertheless, b y c hoosing a sufficiently small sampling period 𝑇 𝑠 , the induced appro ximation error can be k ept within acceptable bounds. 3. A step-wise identification algor ithm A three-step identification algorithm is proposed. First, in Section 3.1 , the initial physical parameters 𝜃 phys are optimized to accurately capture the linear ized input-output beha vior . Second, in Section 3.2 , the nonlinear restor ing f orce is inf erred nonparametrically , after which the parameters 𝛽 are obtained by solving a linear sy stem of equations, while keeping 𝜃 phys fixed. Third, in Section 3.3 , all parameters 𝜃 are jointly optimized to compensate f or bias introduced in the preceding steps and to furt her impro ve the model’ s simulation accuracy . Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 4 of 20 Figure 2: Schematic overview of the BLA framewo rk. Fo r random excitation signals with a Gaussian distribution, the resp onse of a nonlinear system is replaced by the sum of an L TI appro ximation 𝐺 BLA ( 𝑞 ) , unmo deled nonlinear dynamics 𝑦 𝑆 ( 𝑛 ) , and a disturbing noise source 𝑣 ( 𝑛 ) . 3.1. Initial linear model The PISPO system class and the assumed experimental conditions in Section 2 allow us to use the theor y of t he BLA [ 16 , 4 ]. When applied with random ex citation signals that f ollow a Gaussian distr ibution, such as t he random- phase multisine ( 4 ), this framew ork represents t he response of a nonlinear dynamical sy stem as the sum of an L TI appro ximation, unmodeled nonlinear dynamics, and a disturbing noise source. A schematic representation of t his concept is pro vided in Fig. 2 . The BLA itself descr ibes a linearized dynamic relationship between the zero-mean input 𝑢 ( 𝑛 ) and the zero-mean output 𝑦 ( 𝑛 ) that is optimal in the mean-square sense: 𝐺 BLA ( 𝑞 ) ∶ = ar g min 𝐺 𝔼   𝑦 ( 𝑛 ) − 𝐺 ( 𝑞 ) 𝑢 ( 𝑛 )  2  , (9) where 𝔼 [ ⋅ ] denotes the expectation operator and 𝑞 the forw ard-shift operator . As shown in [ 16 ], for the single-input case and under multisine ex citation, the minimizer of ( 9 ) is equiv alent to 𝐺 BLA ( 𝑗 𝜔 𝑘 ) = 𝔼  𝑌 ( 𝑘 ) 𝑈 ( 𝑘 )  , (10) where 𝑌 ( 𝑘 ) and 𝑈 ( 𝑘 ) denote the leakage-free discrete Fourier transf orms (DFTs) of 𝑦 ( 𝑛 ) and 𝑢 ( 𝑛 ) , respectiv ely , at frequency line 𝑘 , theref ore 𝜔 𝑘 = 2 𝜋 𝑘𝑓 0 . The expectation is taken o ver different periods and realizations in  . In our setting, the nonparametric BLA provides an av eraged frequency response matr ix (FRM) that facilitates the subsequent parametr ization of the model components directly related to the NL -LFR input-output behavior , i.e., the state-space matr ices  ,  𝑢 ,  𝑦 , and  𝑦𝑢 . Furt hermore, the BLA estimation procedure yields an assessment of the nonlinear distortions and t he disturbing noise lev el, offer ing valuable insight into the sy stem’ s behavior . These quantities can also be incor porated into t he parameter estimation procedure as weighting functions. For more details on the BLA estimation procedure and its properties, we ref er to [ 16 , 4 ]. Let  𝐺 BLA ( 𝑗 𝜔 𝑘 ) denote the nonparametric estimate of the BLA at frequency 𝜔 𝑘 . The cor responding parametric BLA, parametrized by 𝜃 phys , is given b y 𝐺 ( 𝜁 𝑘 ∣ 𝜃 phys ) = 𝐶 𝑦  𝜁 𝑘 𝐼 − 𝐴  −1 𝐵 𝑢 + 𝐷 𝑦𝑢 , (11) where 𝜁 𝑘 denotes the z-transf or m variable ev aluated on the unit circle at frequency 𝑓 𝑘 = 𝑘𝑓 0 . Starting from the initial guess 𝜃 phys , 0 , the phy sical parameters 𝜃 phys are optimized by minimizing the weighted squared residuals between the nonparametric and parametric BLA estimates ov er the ex cited frequency lines: minimize 𝜃 phys 1     𝑘 ∈     𝑊 ( 𝑘 ) ⊙   𝐺 BLA ( 𝑗 𝜔 𝑘 ) − 𝐺 ( 𝜁 𝑘 ∣ 𝜃 phys )     2 𝐹 , (12) where  ⋅  2 𝐹 represents t he squared Frobenius nor m, ⊙ the element-wise Hadamard product, and  ⊆ {1 , … , 𝑁 ∕2 − 1} the set of e xcited frequency lines, with    denoting its cardinality . Moreover , 𝑊 ( 𝑘 ) denotes a frequency-dependent weighting matr ix, which, in pr inciple, can be any well-chosen real-valued matrix. Y et, we adopt the common choice of weighting the squared model residuals by t he reciprocal of the total variance of the nonparametric BLA estimates (see [ 16 , Ch. 4.3.1]). Such weighting effectively implements a whitening transformation that accounts for unequal variance for unequal variance across frequency lines. If the total variance estimate is not a vailable, i.e., in case of a single multisine experiment, a unity weighting is applied ins tead. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 5 of 20 Remar k 2. In the parametrization of the BLA in ( 12 ), each excited frequency line is treated independently wit hin the optimization problem. This decoupling eliminates the recursion-induced non-con ve xity typically encountered in time-domain parameter estimation, leading to a subs tantially more tractable optimization problem. Remar k 3. The nonparametric BLA is affected by a systematic er ror due to nonlinear distortions [ 16 ]. Conseq uently , the nonparametric BLA is biased, and optimizing 𝜃 phys by minimizing the model residual in ( 12 ) results in biased parameter estimates. This bias conflicts with the objective of identifying the tr ue physical parameters. N ev er theless, at this stage, the pr imar y objective is to achiev e accurate simulation per f or mance. The resulting bias will be addressed and compensated f or in the final step of the proposed algor ithm, described in Section 3.3 . 3.2. Modeling the nonlinear restoring force In this second step, t he nonlinear restoring force and the latent state trajector ies are identified f or fixed 𝜃 phys . To impro ve both the signal-to-noise ratio (SNR) and the computational efficiency we proceed wit h the sample means 𝑢 [ 𝑟 ] ( 𝑛 ) = 1 𝑃  𝑃 −1 𝑝 =0 𝑢 [ 𝑟,𝑝 ] ( 𝑛 ) and 𝑦 [ 𝑟 ] ( 𝑛 ) = 1 𝑃  𝑃 −1 𝑝 =0 𝑦 [ 𝑟,𝑝 ] ( 𝑛 ) . Note that under t he ZOH assumption, a veraging ov er the inputs is redundant as they are noise-free, i.e., 𝑢 [ 𝑟 ] ( 𝑛 ) = 𝑢 [ 𝑟,𝑝 ] ( 𝑛 ) . Nev er theless, we include this step here to accommodate the case where measured (and thus noisy) input data are used instead. 3.2.1. A nonparametr ic sliding window f eedback approach When simulating the output of t he parametric BLA model and compar ing it to the measured output, a mismatch arises due to t he unmodeled nonlinear dynamics. The idea of this step is to obtain a nonparametric estimate of the nonlinear f eedback force 𝑤 that improv es the ag reement between t he simulated and measured outputs. T o this end, we propose a sliding window strategy , in which the nonlinear feedbac k f orce is estimated by solving a local optimization problem ov er a finite time hor izon that is shifted forw ard at each sample instant. Specifically , the nonlinear feedbac k f orce is estimated over a window of length 𝐻 + 1 ≪ 𝑁 , with 𝐻 ∈ ℕ being the prediction horizon lengt h. Only the solution cor responding to the cur rent time instance is retained and used to shift the window forw ard in time by one sample. This pr ocedure is repeated until all 𝑁 time steps hav e been processed. Along the wa y , the full latent s tate 𝑥 is naturally inf er red as well. W e first define the stacked v ectors:  [ 𝑟 ] ( 𝑛 ) =      𝑢 [ 𝑟 ] ( 𝑛 ) 𝑢 [ 𝑟 ] ( 𝑛 + 1) ⋮ 𝑢 [ 𝑟 ] ( 𝑛 + 𝐻 )      ,  [ 𝑟 ] ( 𝑛 ) =      𝑦 [ 𝑟 ] ( 𝑛 ) 𝑦 [ 𝑟 ] ( 𝑛 + 1) ⋮ 𝑦 [ 𝑟 ] ( 𝑛 + 𝐻 )      ,  [ 𝑟 ] ( 𝑛 ) =      𝑤 [ 𝑟 ] ( 𝑛 ) 𝑤 [ 𝑟 ] ( 𝑛 + 1) ⋮ 𝑤 [ 𝑟 ] ( 𝑛 + 𝐻 )      , (13) with dimensions ℝ ( 𝐻 +1) 𝑛 𝑢 , ℝ ( 𝐻 +1) 𝑛 𝑦 , and ℝ ( 𝐻 +1) 𝑛 𝑤 , respectivel y . Then, the optimization problem t hat is solved for each time step 𝑛 ∈ {0 , … , 𝑁 − 1} and eac h realization 𝑟 ∈ {0 , … , 𝑅 − 1} is giv en by minimize  [ 𝑟 ] ( 𝑛 ) 1 2     [ 𝑟 ] ( 𝑛 ) −   [ 𝑟 ] ( 𝑛 )    2  + 𝜆 2     [ 𝑟 ] ( 𝑛 )    2 , (14a) subject t o   [ 𝑟 ] ( 𝑛 ) =  𝑥 𝑥 [ 𝑟 ] ∗ ( 𝑛 ) +  𝑢  [ 𝑟 ] ( 𝑛 ) +  𝑤  [ 𝑟 ] ( 𝑛 ) , (14b) where 𝜆 ∈ ℝ > 0 is a regularization parameter that controls the solution variability , and  is a block -diagonal matr ix f or med b y concatenating 𝐻 + 1 copies of the inv erse of the time-domain sample noise cov ariance matr ix defined in ( 29 ) 6 . Equation ( 14b ) expresses t he stacked simulated output as a function of the initial state 𝑥 [ 𝑟 ] ∗ ( 𝑛 ) , which will be defined shor tly , and the st ac ked inputs  [ 𝑟 ] ( 𝑛 ) and  [ 𝑟 ] ( 𝑛 ) ; the ev olution of t his sequence is described with the f ollowing bloc k -matr ices:  𝑥 =      𝐶 𝑦 𝐶 𝑦 𝐴 ⋮ 𝐶 𝑦 𝐴 𝐻      ,  𝑢 =      𝐷 𝑦𝑢 0 ⋯ 0 𝐶 𝑦 𝐵 𝑢 𝐷 𝑦𝑢 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 𝐶 𝑦 𝐴 𝐻 −1 𝐵 𝑢 𝐶 𝑦 𝐴 𝐻 −2 𝐵 𝑢 ⋯ 𝐷 𝑦𝑢      ,  𝑤 =      𝐷 𝑦𝑤 0 ⋯ 0 𝐶 𝑦 𝐵 𝑤 𝐷 𝑦𝑤 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 𝐶 𝑦 𝐴 𝐻 −1 𝐵 𝑤 𝐶 𝑦 𝐴 𝐻 −2 𝐵 𝑤 ⋯ 𝐷 𝑦𝑤      , (15) 6 In cases where t he sample noise covariance matr ix is not available,  is constructed as a diagonal matrix containing the inverses of the output variances. Doing so ensures that, in the multi-output case, the different outputs are properly scaled relative to their magnitudes. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 6 of 20 with dimensions ℝ 𝑛 𝑦 ( 𝐻 +1)× 𝑛 𝑥 , ℝ 𝑛 𝑦 ( 𝐻 +1)× 𝑛 𝑢 ( 𝐻 +1) , and ℝ 𝑛 𝑦 ( 𝐻 +1)× 𝑛 𝑤 ( 𝐻 +1) , respectiv ely . The optimization problem in ( 14 ) is conv e x. Theref ore, by substituting ( 14b ) into ( 14a ), setting the gradient with respect to  [ 𝑟 ] ( 𝑛 ) to zero, and sol ving for  [ 𝑟 ] ( 𝑛 ) , we obtain the closed-f or m solution:  [ 𝑟 ] ∗ ( 𝑛 ) = −  −1  T 𝑤    𝑥 𝑥 [ 𝑟 ] ∗ ( 𝑛 ) +  𝑢  [ 𝑟 ] ( 𝑛 ) −  [ 𝑟 ] ( 𝑛 )  , (16) where  =  T 𝑤  𝑤 + 𝜆𝐼 . Remar k 4. The regularization term in ( 14a ) is introduced to address two practical issues. First, the measured output data  [ 𝑟 ] ( 𝑛 ) are corr upted by noise, and wit hout regular ization (i.e., when 𝜆 = 0 ) the minimizer of ( 14 ) w ould be highly sensitive to this noise, resulting in larg e variability in the inferred  [ 𝑟 ] ∗ ( 𝑛 ) . Second, when 𝑛 𝑤 > 𝑛 𝑦 , the unregularized optimization problem is inherently ill-posed (i.e.,  is rank -deficient), leading to non-unique solutions; the regularization ter m in ( 14a ) ensures t hat a unique minimizer exists in such cases. The hyperparameter 𝜆 controls the trade-off between data fit and solution v ar iability . Remar k 5. The main motivation for using a prediction horizon lengt h greater than one, or a long er window in general, lies in its ability to mitig ate the effects of measurement noise. When the windo w is short,  [ 𝑟 ] ∗ ( 𝑛 ) ma y become ov erly sensitive to high-frequency noise in  [ 𝑟 ] ( 𝑛 ) , t hus resulting in a ner vous estimate. In contras t, long er windows naturally lead to smoother solutions, as t he individual noisy samples are av eraged out ov er time. The effect of the prediction horizon length 𝐻 and regular ization parameter 𝜆 , including t heir sensitivity to different noise le vels, is empirically studied in Section 4 . Equation ( 16 ) pr ovides the optimal solution o ver the defined window , and is naturally defined f or all 𝑅 realizations pro vided that 𝑛 < 𝑁 − 𝐻 . That is, when 𝑛 ≥ 𝑁 − 𝐻 , computing  [ 𝑟 ] ( 𝑛 ) and  [ 𝑟 ] ( 𝑛 ) requires “future” input-output data be yond the final sample instant. Howe ver , thanks to per iodicity , we are able to define the input-output data at ev ery possible time step 𝑚 ∈ ℤ in terms of t he first 𝑁 samples as follo ws: 𝑢 [ 𝑟 ] ( 𝑚 ) ∶ = 𝑢 [ 𝑟 ] ( 𝑚 mod 𝑁 ) , 𝑦 [ 𝑟 ] ( 𝑚 ) ∶ = 𝑦 [ 𝑟 ] ( 𝑚 mod 𝑁 ) , (17) where mod denotes t he modulus operator . Using ( 17 ), the optimal solution in ( 16 ) is well-defined at each time step and f or ev er y realization. How ev er, in ( 16 ), the initial state of eac h windo w , denoted 𝑥 [ 𝑟 ] ∗ ( 𝑛 ) , is required. The r ecursive ev olution of t his state gives rise to the “sliding” nature of the proposed approach, as detailed ne xt. T o advance to t he ne xt sample, we onl y use the first 𝑛 𝑤 elements of  [ 𝑟 ] ∗ ( 𝑛 ) , which w e denote with 𝑤 [ 𝑟 ] ∗ ( 𝑛 ) . These elements cor respond to the current time step, and are related to the state according to 𝑥 [ 𝑟 ] ∗ ( 𝑛 ) =  𝑥 [ 𝑟 ] 0 , if 𝑛 = 0 𝐴𝑥 [ 𝑟 ] ∗ ( 𝑛 − 1) + 𝐵 𝑢 𝑢 [ 𝑟 ] ( 𝑛 − 1) + 𝐵 𝑤 𝑤 [ 𝑟 ] ∗ ( 𝑛 − 1) , otherwise , (18) which sho ws t hat a value f or 𝑥 [ 𝑟 ] 0 is still needed at 𝑛 = 0 . This initial state is typically unkno wn, and arbitrarily setting it (e.g., to zero) introduces a transient response that is incompatible with the tr ue sys tem, and, therefore, unsuit able for parametric modeling of the nonlinear restoring force. While simpl y discarding the cor responding transient samples is possible, it also results in the loss of valuable inf ormation. Instead, the idea is to use the parametric BLA model to obtain a more inf or med initial state. In addition, the periodicity definition in ( 17 ) can be lev eraged to retain all 𝑁 samples. Specifically , the BLA state trajectories are computed in the frequency domain as 𝑋 [ 𝑟 ] BLA ( 𝑘 ) =  𝜁 𝑘 𝐼 − 𝐴  −1 𝐵 𝑢 𝑈 [ 𝑟 ] ( 𝑘 ) , (19) where 𝑈 [ 𝑟 ] ( 𝑘 ) denotes the DFT of 𝑢 [ 𝑟 ] ( 𝑛 ) at frequency line 𝑘 . T aking the in verse DFT (IDFT) of 𝑋 [ 𝑟 ] BLA ( 𝑘 ) yields 𝑥 [ 𝑟 ] BLA ( 𝑚 ) ∶ = 𝑥 [ 𝑟 ] BLA ( 𝑚 mod 𝑁 ) for 𝑚 ∈ ℤ . This per iodicity is essential, as it allow s t he simulation to be initialized at an arbitrar y phase of the same per iodic trajector y . Concretely , in ( 18 ) we define 𝑥 [ 𝑟 ] 0 ∶ = 𝑥 [ 𝑟 ] BLA (− 𝑁 0 ) and 𝑢 [ 𝑟 ] (0) ∶ = 𝑢 [ 𝑟 ] (− 𝑁 0 ) , where 𝑁 0 ∈ ℕ 0 denotes the number of offset samples required for the transient to decay . The simulation is then r un f or a total of 𝑁 + 𝑁 0 samples, after which the first 𝑁 0 samples of the resulting sequences are discarded. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 7 of 20 In summar y , the proposed procedure descr ibed in t his section yields an analytical expression for t he missing nonlinear restoring f orce, which can be e valuated efficiently ov er the entire time series. F or given values of the horizon length 𝐻 , regularization strength 𝜆 , and offset lengt h 𝑁 0 , advancing t he simulation one sample at a time yields a set of paired samples as  𝑤𝑧 =  𝑤 [ 𝑟 ] ∗ ( 𝑛 ) , 𝑧 [ 𝑟 ] ∗ ( 𝑛 )  𝑁 −1 ,𝑅 −1 𝑛 =0 ,𝑟 =0 , (20) where 𝑧 [ 𝑟 ] ∗ ( 𝑛 ) is easily der iv ed from 𝑥 [ 𝑟 ] ∗ ( 𝑛 ) through ( 6c ). This dat aset will be used in the next step f or parametric modeling of the nonlinear restoring force. 3.2.2. P arametric modeling of the nonlinear restoring force Given  𝑤𝑧 , estimating the nonlinear coefficient matrix 𝛽 reduces to solving a set of linear reg ression problems. Due to the black -box nature of the pol ynomial basis function model, the individual entr ies of 𝛽 do not necessar il y car r y a direct phy sical inter pretation. Ne v er theless, when 𝑛 𝑤 > 1 , it is important t hat 𝛽 respects the physical location of t he nonlinearities. T o enforce this str ucture, we introduce a binar y selection matr ix 𝑃 𝑧,𝑖 ∈ {0 , 1} 𝑛 𝑧,𝑖 × 𝑛 𝑧 , which, f or each location 𝑖 ∈ {1 , … , 𝑛 𝑤 } , extracts the cor rect subvector of dimension 𝑛 𝑧,𝑖 ≤ 𝑛 𝑧 from 𝑧 . The cor responding nonlinear restoring f orce term is then modeled as 𝑤 𝑖 ( 𝑛 ) = 𝛽 T 𝑖 𝜙 𝑖  𝑃 𝑧,𝑖 𝑧 ( 𝑛 )  , (21) where 𝜙 𝑖 ∶ ℝ 𝑛 𝑧,𝑖 → ℝ 𝑛 𝜙,𝑖 is a monomial f eature map, and 𝛽 𝑖 ∈ ℝ 𝑛 𝜙,𝑖 contains the associated coefficients. The full nonlinear mapping in ( 6d ) is subsequentl y defined as a decoupled combination of the individual nonlinearities t hrough 𝛽 ∶ =     𝛽 1 … 0 ⋮ ⋱ ⋮ 0 … 𝛽 𝑛 𝑤     ∈ ℝ 𝑛 𝜙 × 𝑛 𝑤 , 𝜙 ∶ =     𝜙 1 ⋮ 𝜙 𝑛 𝑤     ∶ ℝ 𝑛 𝑧 → ℝ 𝑛 𝜙 , (22) with 𝑛 𝜙 =  𝑛 𝑤 𝑖 =1 𝑛 𝜙,𝑖 . Since ( 21 ) is linear in the unkno wn parameters, each v ector 𝛽 𝑖 can be estimated analyticall y using ordinar y least sq uares on the input-t arg et pairs fr om  𝑤𝑧 , agg regated ov er all 𝑅 realizations. Doing so yields the estimated nonlinear coefficient matrix 𝛽 , which, tog ether with 𝜃 phys , fully initializes the NL-LFR model. Remar k 6. The bias in the 𝜃 phys estimate has an impor tant impact on t he str ucture of the monomials in 𝜙 𝑖 ( ⋅ ) . Specifically , the restoring for ce estimates 𝑤 ∗ inf er red through ( 14 ) absorb the discrepancy betw een the tr ue linear state tra jector ies and those reconstructed using t he es timated linear parameters. As a result, the inferred f orce contains an additional component that depends approximatel y linearl y on t he reconstructed state trajectories. Consequentl y , 𝜙 𝑖 ( ⋅ ) must include degree-one monomials to accuratel y represent the mapping from 𝑧 ∗ to 𝑤 ∗ , ev en though t he tr ue restoring f orce in ( 3d ) corresponds to a purely nonlinear mapping in 𝑧 . This discrepancy mo tivates t he bias-correction procedure and enables reco very of the tr ue linear parameters in the subsequent step. 3.3. Final optimization The previous steps pro vided fully initialized estimates of the parameters 𝜃 = ( 𝜃 phys , 𝛽 ) . In this final step, the goal is twof old: (i) to enhance the simulation per f or mance and (ii) to recov er the tr ue system parameters. These objectives are addressed by sol ving the f ollowing optimization problem: minimize 𝜃 1 𝑅𝑁 𝑅 −1  𝑟 =0 𝑁 ∕2  𝑘 =0  𝑌 [ 𝑟 ] ( 𝑘 ) −  𝑌 [ 𝑟 ] ( 𝑘 ∣ 𝜃 )  2 𝑊 ( 𝑘 ) + 𝛾  𝛽 [1] ( 𝜃 )  1 . (23) The first component in ( 23 ) quantifies t he simulation er ror in the frequency domain. Here, 𝑌 [ 𝑟 ] ( 𝑘 ) is the measured output spectrum, com puted b y appl ying the DFT to 𝑦 [ 𝑟 ] ( 𝑛 ) , while  𝑌 [ 𝑟 ] ( 𝑘 ∣ 𝜃 ) is the modeled output spectrum, obtained by simulating t he NL -LFR model in ( 6 ) wit h input 𝑢 [ 𝑟 ] ( 𝑛 ) and t hen applying t he DFT to t he resulting time-domain output. T o av oid spectral leakage caused by transients from an unknown initial state, we implement the same procedure described in Section 3.2.1 , i.e., simulating and subsequently discarding 𝑁 0 offset samples to ensure steady-state Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 8 of 20 conditions. The proposed frequency-domain e valuation allow s the weighting matrix 𝑊 ( 𝑘 ) to account f or the varying noise characteristics across frequencies. When av ailable, the in verse of the sample noise cov ariance matr ix in ( 30 ) is used at each frequency line. If this matr ix is not av ailable, 𝑊 ( 𝑘 ) is ins tead cons tr ucted as a diagonal matrix cont aining the inv erses of the output variances at each frequency , ensur ing t hat the samples are properl y balanced both across frequencies and, in the multi-output case, across different outputs. The second component in ( 23 ) serves as a regular ization ter m aimed at recov er ing the tr ue underlying linear parameters. As discussed in Remark 6 , t he initial estimates of both 𝜃 phys and 𝛽 include components associated with linear dynamics, which renders the former non-unique and theref ore not physicall y inter pretable. To address this ambiguity , the elements of 𝛽 t hat cor respond to first-degree monomials are stack ed into vector 𝛽 [1] ( 𝜃 ) and penalized using the sparsity-pr omoting 𝓁 1 -norm, denoted by  ⋅  1 . This strategy effectivel y redirects the contribution of 𝛽 [1] ( 𝜃 ) into 𝜃 phys , thereby encouraging the recov ery of the tr ue physical values. The strength of the proposed regular ization term is controlled through hyperparameter 𝛾 ∈ ℝ ≥ 0 . 3.4. Implementation details All algor ithms are implemented in Python. The optimization routines are formulated using JAX [ 2 ] toge ther with Equinox [ 7 ], enabling automatic differentiation and efficient ex ecution on both CPUs and GPUs. Nonlinear least-squares problems are solved using the Lev enberg-Marq uardt algorithm [ 9 , 10 ] as provided by the Optimistix library [ 17 ]. 4. Simulation studies T wo simulation examples are presented to demonstrate t he proposed step-wise algor ithm. The first considers an SDOF Duffing oscillator and analyzes the influence of the different training steps, hyperparameters, and noise lev els on model performance. The second examines a more com plex MDOF sys tem under limited sensing conditions. 4.1. SDOF Duffing oscillator In t his section we study in detail t he proposed step-wise algor ithm on a simulation example of a f orced Duffing oscillator , of which the dynamics are described by the follo wing ODE: 𝑚  𝑦 0 ( 𝑡 ) + 𝑐  𝑦 0 ( 𝑡 ) + 𝑘𝑦 0 ( 𝑡 ) + 𝑘 3 𝑦 3 0 ( 𝑡 ) = 𝑢 ( 𝑡 ) , (24) with 𝑚 , 𝑐 and 𝑘 t he linear mass, damping and stiffness parameters, respectiv ely , and 𝑘 3 the cubic stiffness parameter . Here, the exact noise-free output 𝑦 0 ( 𝑡 ) describes the displacement of the mass, while the input 𝑢 ( 𝑡 ) acts as an e xter nal f orcing term applied to the sys tem. If we define 𝑥 ( 𝑡 ) =  𝑦 0 ( 𝑡 ) ,  𝑦 0 ( 𝑡 )  T , t he Duffing dynamics in ( 24 ) adhere to the NL-LFR str ucture of ( 3 ), wit h  =  0 1 − 𝑘 ∕ 𝑚 − 𝑐 ∕ 𝑚  ,  𝑢 = −  𝑤 =  0 1∕ 𝑚  ,  𝑦 =  𝑧 =  1 0  ,  𝑦𝑢 =  𝑦𝑤 = 0 , (25) and t he typically unknown 𝑓  𝑧 ( 𝑡 )  = 𝑘 3 𝑧 3 ( 𝑡 ) . Within this formulation, the physical parameter v ector is defined as 𝜃 phys = ( 𝑚, 𝑐 , 𝑘 ) . W e generate synt hetic input-output dat a by solving the Duffing equation ( 24 ) using the f ourt h-order Runge-K utta (RK4) integration scheme. The parameters of the Duffing oscillator are set to 𝑚 = 1 k g , 𝑐 = 2 N s m −1 , 𝑘 = 100 N m −1 , and 𝑘 3 = 500 N m −3 . The training data consists of 𝑅 = 5 r ealizations of 𝑃 = 2 steady-state per iods of a random-phase multisine ( 4 ), with 𝑁 = 8192 samples each. The multisine input signal ex cites frequencies up to 10 Hz , is sampled at 𝑓 𝑠 = 128 Hz , and has individual amplitudes 𝑈 𝑘 chosen such that an ov erall root mean square (RMS) amplitude of 12 N is obtained. Star ting from noise-free samples, we generate t hree datasets with increasing output noise lev els, ranging from almost noise-free to strongl y cor r upted. Specifically , white Gaussian noise is added to t he noise-free dat a to achie ve SNRs of 60 dB , 40 dB and 20 dB . From these SNRs, low er bounds on t he achie vable relative simulation er rors can be derived, corresponding to 0 . 1 % , 1 % and 10 % , respectivel y 7 . 7 These bounds follo w directly fr om the definition SNR = 10 log 10 ( 𝑃 s ∕ 𝑃 n ) , wher e 𝑃 s and 𝑃 n denote t he signal and noise pow er, respectivel y . The low er bound is then defined as the noise-to-signal amplitude ratio, i.e.  𝑃 n ∕  𝑃 s = 10 − SNR ∕20 , expressed as a percentage. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 9 of 20 T able 1 P arameter estimates of the Duffing oscillator obtained using the BLA procedure for 100 different initializations in ( 26 ) . The estimates (mean ± standard deviation) a re highly consistent across initializations and noise levels, but bias is present, most noticeable in the stiffness pa rameter 𝑘 , which is systematically overestimated due to the hardening effect of the cubic sp ring. 𝑚 [ k g ] 𝑐 [ N s m −1 ] 𝑘 [ N m −1 ] true 1.00 2.00 100 SNR 60 dB 0 . 988 ± 5.77 ×10 −6 2.10 ± 6.34 ×10 −6 114 ± 6.44 ×10 −4 SNR 40 dB 0.988 ± 1.01 ×10 −5 2.10 ± 8.07 ×10 −6 114 ± 1.13 ×10 −3 SNR 20 dB 0.987 ± 1.53 ×10 −6 2.10 ± 5.77 ×10 −6 114 ± 1.76 ×10 −4 T able 2 Simulation output normalized ro ot mean square errors (NRMSEs) for the Duffing oscillato r at different stages of the step-wise algorithm; compa red to their resp ective theoretical low er b ounds. The initial NL-LFR mo dels already show a significant improvement over the BLA; for SNRs of 60 dB and 40 dB , the lo wer b ounds are not met, likely b ecause the ZOH discretization erro r exceeds the resp ective noise levels. NRMSEs [ % ] BLA initial NL-LFR optimized NL-LFR lo wer bound SNR 60 dB 20.5 4.38 1.08 0.10 SNR 40 dB 20.5 4.51 1.48 1.00 SNR 20 dB 22.8 10.8 10.0 10.0 The goal of this simulation example is t hreef old: first, to assess how eac h training s tep (BLA estimation, res tor ing f orce modeling, and final optimization) affects the parameter estimates and contr ibutes to the model’ s simulation performance; second, to study the influence of the regular ization strength 𝜆 and prediction hor izon lengt h 𝐻 on t he restoring f orce approac h; and third, to ev aluate the sensitivity of t he abo v e analy ses with respect to the different le vels of measurement noise. Step I: Best linear approximation The first step consists of estimating the parameters of the linear par t of the model, as outlined in Section 3.1 . The concrete outcome is a linearized state-space model with estimates of t he parameters 𝑚 , 𝑐 , and 𝑘 . To assess the sensitivity of the parameter estimates to initialization, w e repeat t he optimization procedure 100 times, each with a maximum of 100 Lev enberg-Mar quardt iterations, starting from different initial guesses generated as f ollow s:     𝑚 0 𝑐 0 𝑘 0     =     𝑚 (1 + 𝛿 𝑚 ) 𝑐 (1 + 𝛿 𝑐 ) 𝑘 (1 + 𝛿 𝑘 )     , where 𝛿 𝑚 , 𝛿 𝑐 , 𝛿 𝑘 ∼  (−0 . 9 , 0 . 9) , (26) with  (−0 . 9 , 0 . 9) denoting independent random variables sampled from a uniform distribution ov er the interval [−0 . 9 , 0 . 9] . In other w ords, each parameter is randoml y initialized within ±90 % of its tr ue v alue. The identification procedure was then carr ied out as described above, and the resulting parameter estimates are summarized in Table 1 . T wo main obser vations can be drawn from these results. First, the estimates are highly consistent across the 100 r uns and noise lev els, with small standard deviations indicating robust con ver gence to a unique solution. Second, the estimates are generall y close to t he tr ue values, but a bias is present, as expected. This bias is most noticeable in t he stiffness parameter 𝑘 , which is systematicall y ov erestimated due to the hardening effect of the cubic stiffness term 𝑘 3 . For furt her anal yses, we proceed with the mean v alues of the estimates in T able 1 . W e assess t he simulation perf or mance of the obt ained BLA model using the nor malized RMS er ror (NRMSE), computed here on the output signal 8 . The results, sho wn in Table 2 , indicate a consistent NRMSE of around 20 % across all noise le vels. This relativel y high error confir ms that the linear model f ails to adequatel y capture the Duffing oscillator’ s nonlinear beha vior . 8 In general, the NRMSE is computed as t he RMS value of the respective er ror signal divided by the RMS value of the ref erence signal, multiplied by 100 % , thereby providing a relative and easily inter pretable measure of performance. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 10 of 20 − 20 0 20 SNR 60 dB force w [N] H = 1 , λ = 1e − 6 NRMSE 37.2% H = 10 , λ = 1e − 4 NRMSE 24.2% H = 25 , λ = 1e − 2 NRMSE 62.0% − 20 0 20 SNR 40 dB force w [N] NRMSE 90.0% NRMSE 27.8% NRMSE 62.0% − 0 . 5 0 . 0 0 . 5 displacement z [m] − 20 0 20 SNR 20 dB force w [N] NRMSE 99.8% − 0 . 5 0 . 0 0 . 5 displacement z [m] NRMSE 81.5% − 0 . 5 0 . 0 0 . 5 displacement z [m] NRMSE 63.1% inferred v alues 3rd-degree p olynomial fit Figure 3: Visualization of  𝑤𝑧 and corresponding p olynomial fit for three hyperparameter combinations across all noise levels. The inferred values, obtained by solving ( 14 ) , generally suggest the correct cubic nature of the resto ring force, but the influence of noise and hyp erpa rameters is evident. Step II: Restoring force modeling In t his second step, we address two main questions: (i) can the sliding-window approach suggest a restoring f orce signal that dr iv es the simulated output close to the desired output, and (ii) can a static nonlinear mapping 𝑓 ∶ ℝ 𝑛 𝑧 → ℝ 𝑛 𝑤 be obtained from this signal? Par ticular attention is paid to the influence of the regularization parameter 𝜆 , the prediction hor izon lengt h 𝐻 , and t he impact of measurement noise on t he results. In the follo wing, the offset length is set to 𝑁 0 = 100 samples. Figure 3 visualizes the inf er red dataset  𝑤𝑧 (in blue) f or t hree hyperparameter settings across all noise lev els, where the nonparame tr ic restoring f orce 𝑤 is plotted as a function of the displacement 𝑧 . Although not visible from the plots, in each case the inf er red restoring f orce, obtained by solving ( 14 ), has successfull y driven the simulated output close to the desired output. When plotted against displacement, this nonparametric f orce generally rev eals the expected cubic nature of the underlying system nonlinear ity . Howe v er, the cor responding t hird-degree odd polynomial 9 fits (in orange) and their NRMSEs indicate that the relationship between 𝑧 and 𝑤 does not alw ays reflect a clear st atic mapping. This phenomenon is most pronounced in the bottom-left plot, cor responding to the highest noise lev el combined with the shortest window length and weakes t regularization. Here, the result lacks any meaningful structure, making it generally impossible to e xtract a reliable mapping from the signals in  𝑤𝑧 . A more detailed analy sis is car ried out by performing a g rid search ov er t he hyperparameters. The prediction horizon length 𝐻 is varied linearly from 1 to 50, while the regular ization parameter 𝜆 is varied logar ithmically from 10 −7 to 10 −1 , using 50 values in total. The results are summar ized in Fig. 4 , where the left column show s the 9 Specifically , 𝜙 ( 𝑧 ) = [ 𝑧, 𝑧 3 ] T , as this structure matches the tr ue underlying nonlinearity and includes the linear ter m required to compensate f or the bias in t he 𝜃 phys estimate (see Remark 6 ). Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 11 of 20 10 20 30 40 50 SNR 60 dB horizon H [-] output trac king p olynomial fit 10 20 30 40 50 SNR 40 dB horizon H [-] 10 − 6 10 − 4 10 − 2 regularization λ [N − 2 ] 10 20 30 40 50 SNR 20 dB horizon H [-] 10 − 6 10 − 4 10 − 2 regularization λ [N − 2 ] 5 10 15 20 NRMSE [%] 20 40 60 80 NRMSE [%] Figure 4: Grid search over the horizon length 𝐻 and regularization strength 𝜆 . The left column shows nonparametric NRMSEs betw een measured and simulated output signals, while the right column shows pa rametric NRMSEs of the co rresp onding polynomial fits. Red dots ma rk hyperparameter combinations where the nonparametric NRMSEs a re close to their lo wer b ounds; these combinations also yield the most accurate p olynomial fits. nonparametric NRMSEs betw een the measured and simulated outputs, while the right column presents the NRMSEs of third-degree odd polynomial fits to the inferred restoring f orce signals, similar to Fig. 3 . The red dots mark the hyperparameter combinations f or which the nonparametr ic output NRMSEs are considered close 10 to their lower bounds der ived from t he SNR. When ex amining 𝜆 first, Fig. 4 indicates that smaller values tend to reduce the nonparametric output NRMSE. Howe ver , this does not guarantee a good pol ynomial fit, since insufficient regularization causes the sliding-window algorit hm to ov erfit the noise. The best polynomial fits are achie ved at combinations that render t he nonparametr ic output NRMSEs close to their respective low er bounds. Here, regularization suppresses noise just enough without being too res tr ictiv e. As for the prediction hor izon length 𝐻 , its effect is less pronounced, but it is clear that in gener al the windo w lengt h should at leas t be a f ew samples long. Moreov er , it can be observ ed that longer window s are preferred for higher noise lev els. As discussed in R emark 5 , this is because longer window s allow the algorit hm to a verag e out the high-frequency noise, leading to smoother estimates of t he restoring force signal. Note that, in terms of computational cost, longer window s are somewhat more expensive to sol ve, but since the optimization problem admits an anal ytical solution, this increase is limited. Finall y , we assess t he simulation performance of the initial NL-LFR models obtained from t he hyperparameter combination 𝐻 = 10 and 𝜆 = 10 −4 , i.e., the second column in Fig. 3 . The corresponding simulation output NRMSEs, presented in T able 2 , sho w a per f or mance gain of more than a factor of f our f or t he SNRs of 60 dB and 40 dB , but furt her improv ement is s till required. At t he SNR of 20 dB , the simulation accuracy is already close to the low er bound 10 Defined using a relative tolerance of 5 % and an absolute tolerance of 0 . 02 % , with the larger of the two used as threshold. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 12 of 20 2 4 SNR 60 dB NRMSE [%] BLA NRMSE 20.5% 100 110 k [N/m] 2 4 SNR 40 dB NRMSE [%] BLA NRMSE 20.5% 100 110 k [N/m] 10 0 10 1 10 2 iteration [-] 10 . 0 10 . 5 SNR 20 dB NRMSE [%] BLA NRMSE 22.8% 100 110 k [N/m] Figure 5: Evolution of the output simulation NRMSE and the linea r stiffness parameter 𝑘 over the iterations of the final optimization step. It can b e seen that the optimization procedure first prio ritizes reducing the simulation erro r, follow ed b y correcting the bias in the stiffness parameter. T able 3 Physical parameter estimates of the Duffing oscillato r after the final optimization step. Compared with T able 1 , the bias in the stiffness parameter is almost entirely eliminated. 𝑚 [ k g ] 𝑐 [ N s m −1 ] 𝑘 [ N m −1 ] 𝑘 3 [ N m −3 ] true 1.00 2.00 100 500 SNR 60 dB 0.999 2.06 99.9 485 SNR 40 dB 0.999 2.06 100 485 SNR 20 dB 0.999 2.06 100 487 of 10 . 0 % , despite t he relativel y larg e er ror obser ved in t he polynomial fit of Fig. 3 . This pol ynomial model nonetheless captured the ke y cubic trend of the restor ing force, which pro ves sufficient f or accurate simulation. W e proceed wit h the abov e initial NL-LFR models to t he next identification s tep. Step III: F inal optimization The final identification step aims to impro ve the simulation performance of the initial NL -LFR models, while simultaneously recov er ing t he tr ue parameters of the Duffing oscillator . This goal is achie v ed by minimizing the dual-objective cost function in ( 23 ) using 100 Lev enberg-Mar quardt iterations with a regular ization parameter of 𝛾 = 5 × 10 −3 . Figure 5 show s the ev olution of the output NRMSEs and the stiffness parameter 𝑘 ov er the iterations. It can be obser v ed that the optimization procedure initially focuses on reducing t he NRMSE, and subsequentl y on cor recting the parameter bias. For all three SNRs, the NRMSE is significantly reduced, while the bias on 𝑘 is almost entirely remo ved, as also apparent from the final parameter estimates in T able 3 . The final output NRMSEs in T able 2 confir m a substantial impro vement in simulation per f or mance. Ho wev er, f or SNRs of 60 dB and 40 dB , the NRMSEs remain abo ve the ideal v alues e xpected from t he noise le vels. Since the model structure is correctly specified and noise is limited, this discrepancy is most plausibly due to the ZOH discretization, which does not ex actly replicate the underlying continuous-time ODE, as discussed in Remark 1 . This discretization mismatch ma y also explain wh y not all parameters in T able 3 are recov ered e xactly . Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 13 of 20 − 50 0 50 magnitude [dB] − 0 . 2 0 . 0 0 . 2 SNR 60 dB displacement [m] NRMSE 18.5 to 1.07% − 50 0 50 magnitude [dB] − 0 . 2 0 . 0 0 . 2 SNR 40 dB displacement [m] NRMSE 18.6 to 1.47% 0 5 10 frequency [Hz] − 50 0 50 magnitude [dB] 0 20 40 60 time [s] − 0 . 2 0 . 0 0 . 2 SNR 20 dB displacement [m] NRMSE 21.1 to 10.23% measured nonlinear distortion BLA residual noise distortion NL-LFR residual Figure 6: T est data: simulation p erfo rmance of the Duffing oscillator BLA and NL-LFR mo dels in time and frequency domains across different noise levels. Although the NL-LFR mo dels outp erfo rm the BLA, the noise levels are not fully met fo r SNR of 60 dB and 40 dB , most plausibly due to the ZOH discretization. P er formance on test data W e conclude by evaluating t he performance of the trained NL-LFR models on unseen test data. The test set is a different random-phase realization of a multisine wit h the same properties as the training dat a. Figure 6 presents the simulation results in both time and frequency domains f or all noise lev els. The nonlinear and noise distortion levels are computed from the training data via t he BLA estimation procedure described in [ 16 , Ch. 4.3.1]. The per f ormance is comparable to that on the training dat a, with time-domain output NRMSEs similar to those in Table 2 . In the frequency domain, optimal performance is, indeed, not achie ved f or SNRs of 60 dB and 40 dB , as t he noise lev el is not met. The frequency-domain plots provide an additional validation of the BLA approac h, as the residual between t he measured response and the fitted linear model aligns with the nonparametrically estimated variance estimates. This observation indicates that the BLA captures the full linear contr ibution of the system dynamics, with the remaining mismatch attributable to nonlinear distor tions and measurement noise. 4.2. MDOF mass-spring-damper system In the previous simulation ex ample, we considered an SDOF sys tem where the nonlinearity depended directly on the measured output. The follo wing simulation study inv estigates an MDOF system characterized by a nonlinear restoring f orce dependent on an unmeasured state. Under these conditions, classical RFS approac hes are not applicable. A schematic overview of the MDOF system is shown in Fig. 7 . The system comprises tw o masses, 𝑚 1 and 𝑚 2 , connected by linear springs and dampers, wit h the first mass linked to the ground via a nonlinear restoring f orce. The second mass ser ves as both the ex citation and measurement point. The cor responding physical parameter values are listed in T able 4 . Defining 𝑧 ( 𝑡 ) = [ 𝑥 1 ( 𝑡 ) ,  𝑥 1 ( 𝑡 )] T , the nonlinear restoring f orce is modeled as t he sum of a cubic spr ing Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 14 of 20 Figure 7: Schematic overview of the MDOF system. We excite and measure at the second mass, but the nonlinearit y is b et ween the first mass and the ground, which makes classical RFS approaches inapplicable. T able 4 T rue and estimated linea r pa rameters of the MDOF mass-spring-damper system shown in Fig. 7 . The value b et ween pa rentheses indicates the effective linear damping co efficient that takes into account the contribution of the nonlinear damping term a round zero velo city . 𝑚 1 [ k g ] 𝑚 2 [ k g ] 𝑐 1 [ N s m −1 ] 𝑐 2 [ N s m −1 ] 𝑘 1 [ N m −1 ] 𝑘 2 [ N m −1 ] true 2.00 1.00 5.00 (26.0) 2.00 800 600 estimated 1.97 0.998 24.2 1.97 796 600 and a damper with smooth saturation: 𝑓  𝑧 ( 𝑡 )  = 𝛼 1 t anh  𝛼 2 𝑧 2 ( 𝑡 )  + 𝛼 3 𝑧 3 1 ( 𝑡 ) , (27) with parameters 𝛼 1 = 7 . 0 , 𝛼 2 = 3 . 0 , and 𝛼 3 = 5 . 0 × 10 4 . The cor responding continuous-time state-space matr ices in ( 3 ) can be obtained directly from the schematic in F ig. 7 and are therefore omitted f or brevity . Synthetic, noiseless input-output data are generated by numerically solving the continuous-time system using the RK4 integration scheme. The training set comprises 𝑅 = 6 realizations, each containing 𝑃 = 1 steady-state period of a random-phase multisine ( 4 ), wit h 𝑁 = 8192 samples per realization. The multisine input signal ex cites frequencies up to 10 Hz , is sampled at 𝑓 𝑠 = 128 Hz , and its individual amplitudes 𝑈 𝑘 are chosen to yield an overall RMS am plitude of 10 N . The output is defined as the displacement of the second mass. The initial values of the linear parameters are defined as     𝑚 𝑖, 0 𝑐 𝑖, 0 𝑘 𝑖, 0     =     𝑚 𝑖 (1 + 𝛿 𝑚 𝑖 ) 𝑐 𝑖 (1 + 𝛿 𝑐 𝑖 ) 𝑘 𝑖 (1 + 𝛿 𝑘 𝑖 )     , wher e 𝛿 𝑚 𝑖 , 𝛿 𝑐 𝑖 , 𝛿 𝑘 𝑖 ∼  (−0 . 9 , 0 . 9) , (28) f or 𝑖 ∈ {1 , 2} . W e per f or m 10 optimization runs wit h a maximum of 100 Lev enberg-Marq uardt iterations each, starting from different initial guesses, and proceed to the subsequent steps with the best per f orming model. Next, f or the modeling of the res tor ing f orce, we use 𝐻 = 15 , 𝑁 0 = 100 , and 𝜆 = 10 −8 in the nonparametric inf erence step. An odd polynomial of deg ree 7, wit hout cross-terms, is then fitted to the inferred dat a. The results of both steps are shown in F ig. 8 . Remar k 7. Dur ing the nonparametric inf erence step, it is not necessar y to specify whether the nonlinear restoring f orce or iginates from the spring, the damper , or a combination of both. Only the location of the nonlinearity needs to be specified. The precise nature of the nonlinear ity can be decided later when fitting the polynomial model. This decision can be guided by visual inspection of a 3D plot such as Fig. 8 , or through cross-validation ov er different candidate input configurations. In the final optimization step, a maximum of 100 Lev enberg-Marq uardt iterations is per f ormed with bias regularization parameter 𝛾 = 10 −5 . The resulting parameter estimates are listed in T able 4 . Ov erall, the estimates are close to the tr ue v alues, e x cept f or t he dam ping coefficient 𝑐 1 , whic h appears to be ov erestimated by near ly a f actor of five. This beha vior is, ho we ver , e xpected given the explicit separation of t he restoring force into linear and nonlinear components. In par ticular , the nonlinear saturation ter m 𝐹 nl = 𝛼 1 t anh( 𝛼 2  𝑥 1 ) in ( 27 ) exhibits a non-negligible linear contribution around zero velocity . Linear izing 𝐹 nl with respect to  𝑥 1 at  𝑥 1 = 0 yields a slope of 21 . 0 N s m −1 , which, Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 15 of 20 − 0 . 05 0 . 00 0 . 05 displacement z 1 [m] − 1 0 1 velocity z 2 [m/s] − 40 − 20 0 20 force w [N] Figure 8: Nonlinea r resto ring fo rce modeling for the MDOF mass-sp ring-damp er system. The red dots rep resent the inferred nonpa rametric samples, while the smo oth manifold represents the fitted p olynomial mo del. − 1 0 1 velocity ˙ x 1 [m/s] − 10 0 10 damping force [N] F lin = c 1 ˙ x 1 ˆ F lin = ˆ c 1 ˙ x 1 F nl = α 1 tanh( α 2 ˙ x 1 ) ˆ F nl = ˆ β > φ ([0 , ˙ x 1 ] > ) F tot = F lin + F nl Figure 9: T otal damping fo rce decomp osed into linear and nonlinear contributions. The apparent overestimation of  𝑐 1 should b e interpreted as the effective linear damping around zero velo cit y , which includes contributions from b oth the linea r damp er and the nonlinear saturation term. when added to the tr ue linear damping 𝑐 1 = 5 . 0 N s m −1 , results in an effective linear damping of 26 . 0 N s m −1 . This value is significantly closer to t he estimated damping coefficient of 24 . 2 N s m −1 . This effect is fur ther illustrated in Fig. 9 , which depicts the total damping f orce 𝐹 t ot as t he sum of the linear contr ibution 𝐹 lin = 𝑐 1  𝑥 1 and t he nonlinear saturation term 𝐹 nl = 𝛼 1 t anh( 𝛼 2  𝑥 1 ) . The estimated effectiv e linear damping  𝑐 1 aligns with the t angent of 𝐹 t ot in the low -velocity regime. W e validate t he trained NL -LFR models on test data from a multisine realization wit h the same proper ties as the training set. Figure 10 show s the simulation results in both the time and frequency domains. The final NL -LFR model achiev es an improv ement of roughly a fact or of ten compared to the BLA model. The figure also includes the initial NL -LFR model from the restor ing force modeling step, which already performs close to t he final model, indicating that the final optimization mainly serves f or bias cor rection. Furt her impro vement is likel y hindered by the ZOH discretization errors discussed previousl y . Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 16 of 20 0 20 40 60 time [s] − 0 . 1 0 . 0 0 . 1 displacement [m] NRMSEs: 10.1, 1.08, and 0.90% 0 5 10 frequency [Hz] − 75 − 50 − 25 0 25 magnitude [dB] measured initial NL-LFR residual BLA residual optimized NL-LFR residual Figure 10: T est data: simulation p erfo rmance of the MDOF mass-spring damp er BLA and NL-LFR mo dels in time and frequency domains. The initial NL-LFR model from the resto ring fo rce modeling step already performs close to the final mo del, which cannot b e improved further due to mo deling errors originating from the ZOH discretization. 5. Experimental results on the Silv erbox benchmark system The proposed me thod is ev aluated on the e xperiment al data of the Sil verbo x benchmar k [ 21 ], which represents an electronic implement ation of the Duffing oscillator ( 24 ). In this experiment, the input voltage acts as an excitation analogous to mechanical force, whereas t he measured output voltag e represents t he displacement response. The recorded data f or m an ar row -shaped trajector y consisting of two distinct par ts and are sampled at approximatel y 610 Hz . The first part, ref er red to as the “ar ro whead” and used ex clusively for testing, contains 40 000 samples of a white Gaussian noise signal wit h gradually varying amplitude, filtered using a 9th-order Butterwort h filter wit h cutoff frequency 200 Hz . The remaining data cor respond to ten consecutive realizations of a random-phase multisine ex cit ation, together comprising 86 750 samples and e xciting only odd har monics up to 200 Hz . From t his multisine portion, the final 21 688 samples are reserved f or testing purposes. Successiv e realizations are separated b y shor t zer o- input intervals and slightl y ex ceed one full period length to accommodate transient samples, which are remo v ed pr ior to analy sis to ensure steady-state operation. The data can be considered nearly noise-free [ 21 ]; yet, it is not possible to compute the sample noise variance of the output, as each multisine realization consists of only a single per iod (see Appendix A ). The experiment al Silv erbox data deviate from ideal identification conditions in two respects. The input signal is measured af ter passing through a low -pass filter [ 21 ], such that it no longer satisfies the ZOH excitation of Assumption 2 . Furt hermore, the sampling frequency is only about t hree times higher than the maximum ex cited frequency . In sys tem identification practice, sampling rates betw een ten and twenty times t he highest frequency of interest are commonly recommended to accurately appro ximate continuous-time dynamics and to limit discretization er rors. The relativel y low sampling rate therefore introduces additional modeling inaccuracies. T o alleviate t hese effects, the data are upsampled by a f actor of twenty using cubic spline interpolation, after whic h simulation er rors are ev aluated on the downsampled signals. As the Silv erbo x implements the dynamics of a Duffing oscillator, we assume the model structure in ( 25 ). W e follo w the step-wise procedure using the same optimization settings as in the previous sections. For BLA estimation, the initial values of the linear parameters are taken from [ 8 ]. In the sliding-window approach, we use a prediction hor izon lengt h of 𝐻 = 10 , an offset length 𝑁 0 = 100 , and a regularization strength of 𝜆 = 10 −1 . The static nonlinear ity is modeled using a third-order polynomial, and a bias cor rection with 𝛾 = 10 −1 is applied dur ing the final optimization step. Both optimization steps used a maximum of 100 Le venber g-Marquardt iterations. Figure 11 presents the simulation results on the arrowhead test data. The lef t plot sho ws the time-domain performance of the BLA and NL-LFR models, while the r ight plot displays the cor responding nonlinear restoring f orce as simulated by the NL -LFR model, which is clearly cubic. The NL-LFR output simulation er ror is reduced by more than a factor of sixteen compared to the BLA model. The per f or mance increase is also apparent from T able 5 , which sho ws the simulation er rors of both models on both test datasets. It is wort h noting that the NL -LFR er rors are very consistent, e ven in the e xtrapolation region, thanks to t he correctly assumed nonlinear structure. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 17 of 20 0 20000 40000 sample index [-] − 0 . 2 0 . 0 0 . 2 displacement [V] NRMSE 28.7 to 1.73% − 0 . 2 0 . 0 0 . 2 displacement z [V] − 0 . 10 − 0 . 05 0 . 00 0 . 05 0 . 10 force w [V] measured nonlinear force BLA residual start extrapolation NL-LFR residual Figure 11: Simulation p erformance of the BLA and NL-LFR mo dels on the Silverb ox arro whead test data (left), and the co rresp onding nonlinear restoring force as simulated by the NL-LFR mo del (right). T able 5 BLA and NL-LFR simulation errors on the Silverb ox test datasets. The NL-LFR mo del remains ac curate even in the extrap olation region of the arro whead data thanks to its correctly assumed nonlinear structure. BLA NL-LFR NRMSE RMSE NRMSE RMSE multisine 16 . 1 % 8 . 72 mV 1 . 54 % 0 . 825 mV arro whead (full) 28 . 7 % 15 . 3 mV 1 . 74 % 0 . 927 mV arro whead (no extrap olation) 19 . 7 % 8 . 47 mV 1 . 61 % 0 . 684 mV T able 6 BLA and NL-LFR output simulation errors on the Silverb o x test datasets when the data are upsampled b y only a factor of five instead of t went y . Compared to T able 5 , the errors increase significantly , especially for the NL-LFR mo del, highlighting the imp o rtance of proper upsampling. BLA NL-LFR NRMSE RMSE NRMSE RMSE multisine 16 . 8 % 9 . 11 mV 5 . 25 % 2 . 85 mV arro whead (full) 29 . 1 % 15 . 4 mV 5 . 89 % 3 . 14 mV arro whead (no extrap olation) 20 . 3 % 8 . 70 mV 5 . 30 % 2 . 28 mV The results in T able 5 outper f or m repor ted results in studies that identify ph ysical model parameters from data, although there are not many physics-based methods t o compare wit h. In [ 18 ], a 2 . 1 mV output RMS error (RMSE) was reported on the arrowhead dat a, while [ 8 ] reports a 2 . 9 % output NRMSE on a fragment of t he multisine data. As all three methods adopt the same model str ucture, the discrepancies between t he proposed method and the results repor ted in [ 18 ] and [ 8 ] are likel y attributable to differences in the discretization sc hemes used during optimization. This claim is suppor ted by T able 6 , which reports simulation er rors on test data when t he training data were upsampled by only a factor of five instead of twenty . For t he NL-LFR model in par ticular, the er rors increased by more t han a factor of three, emphasizing that proper upsampling is essential f or simulation per f ormance. 6. Conclusion This wor k introduced a nov el approac h to identifying and modeling the nonlinear restor ing f orce in MDOF systems f or mulated as NL-LFR state-space models. In contras t to traditional RFS approaches, t he proposed method relax es t he measurement assumptions considerabl y . In particular, any output quantity can be measured (displacement, v elocity , or acceleration), only a subse t of the deg rees of freedom needs to be measured, and the measurements ma y be noisy . Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 18 of 20 Starting from an initial linear model, the nonlinear restoring force is reconstr ucted from measured data using a sliding-window approach. The resulting nonparametric estimate is subsequently used to identify the static nonlinear mapping within the NL -LFR structure by fitting a polynomial basis function model. Since both the nonparametric inf erence and nonlinear parametrization admit closed-f orm solutions, t he ov erall identification procedure remains computationally efficient and straightf or w ard to implement. A final optimization stage improv es simulation accuracy while compensating f or any bias in the estimated ph ysical parameters. V alidation on simulated e xamples and experimental Sil verbo x data confirms accurate nonlinear model identification, with remaining discrepancies pr imarily attributed to unav oidable discretization effects. Although the proposed algor ithm relies on per iodic multisine data, this requirement mainly suppor ts the frequency - domain estimation of the BLA and the final optimization stage. The central contr ibution of this w ork, the sliding- window estimation of the nonlinear restoring for ce, is not intrinsically tied to per iodic data. As a result, e xtending the method to arbitrar y input-output measurements is relativ ely straightf or ward b y car rying out the remaining estimation steps directl y in the time domain. A. Computing the sample noise cov ariance Given 𝑃 > 1 , we can compute sample noise covariance matr ices t hat quantify the disturbing noise source 𝑣 ( 𝑛 ) . These matr ices are useful as they (i) provide insight into the noise properties, (ii) form t he basis f or defining w eighting matrices in t he parameter estimation pr ocedure, and (iii) suppor t model validation. In the time domain, the sample noise cov ar iance matrix is computed as  Σ time 𝑦 = 1 𝑁 𝑅 ( 𝑃 − 1) 𝑁 −1  𝑛 =0 𝑅 −1  𝑟 =0 𝑃 −1  𝑝 =0  𝑦 [ 𝑟,𝑝 ] ( 𝑛 ) − 𝑦 [ 𝑟 ] ( 𝑛 )  𝑦 [ 𝑟,𝑝 ] ( 𝑛 ) − 𝑦 [ 𝑟 ] ( 𝑛 )  T , (29) where 𝑦 [ 𝑟 ] ( 𝑛 ) = 1 𝑃  𝑃 −1 𝑝 =0 𝑦 [ 𝑟,𝑝 ] ( 𝑛 ) denotes the sample mean ov er t he periods. In t he frequency domain, the sample noise cov ariance matr ix at frequency line 𝑘 is computed as  Σ freq 𝑦 ( 𝑘 ) = 1 𝑅 ( 𝑃 − 1) 𝑅 −1  𝑟 =0 𝑃 −1  𝑝 =0  𝑌 [ 𝑟,𝑝 ] ( 𝑘 ) − 𝑌 [ 𝑟 ] ( 𝑘 )  𝑌 [ 𝑟,𝑝 ] ( 𝑘 ) − 𝑌 [ 𝑟 ] ( 𝑘 )  H , (30) where 𝑌 [ 𝑟 ] ( 𝑘 ) = 1 𝑃  𝑃 −1 𝑝 =0 𝑌 [ 𝑟,𝑝 ] ( 𝑘 ) denotes the sample mean over the periods. Note t hat in t he time domain, we a verag e ov er all samples as t he noise is assumed stationar y , whereas in t he frequency domain, eac h frequency is treated separately since the noise ma y be frequency -dependent. CRedi T authorship contribution statement Merijn Floren: Conceptualization, For mal anal ysis, Me thodology , Softw are, V alidation, V isualization, W riting – original draf t preparation. Jan Sw ev ers: Funding acquisition, R esources, Super vision, W riting – review and editing. Declaration of Generativ e AI and AI-assisted technologies in the writing process During the preparation of this work, the authors used ChatGPT to assist with language refinement, grammatical cor rections, and improv ements to clarity and conciseness of the manuscr ipt te xt. After using this tool, the aut hors revie wed and edited all content as needed and take full responsibility f or t he content of the publication. Ref erences [1] Bonisoli, E., V igliani, A., 2007. Identification techniques applied to a passive elasto-magnetic suspension. Mechanical Sy stems and Signal Processing 21, 1479–1488. [2] Bradbury , J., Frostig, R., Hawkins, P ., Johnson, M.J., Leary , C., Maclaurin, D., Necula, G., Paszke, A., V anderPlas, J., W anderman-Milne, S., Zhang, Q., 2018. J AX: composable transf or mations of Python+NumPy programs. URL: http://github.com/google/jax . softw are a vailable from http://github.com/google/ jax . [3] Dossogne, T., Noel, J.P ., Kerschen, G., Peeters, B., Debille, J., V aes, M., Schoukens, J., 2015. Nonlinear g round vibration identification of an F-16 aircr aft. Part II: Understanding nonlinear beha viour in aerospace structures using sine-sw eep testing, in: International Forum of Aeroelasticity and Structural Design (IFASD), pp. 751–769. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 19 of 20 [4] Enqvist, M., Ljung, L., 2005. Linear approximations of nonlinear FIR systems for separable input processes. Automatica 41, 459–473. [5] Floren, M., Noël, J.P ., 2022. Nonlinear restoring f orce modelling using Gaussian processes and model predictive control, in: Conference Proceedings of ISMA2022-USD2022, pp. 2493–2498. [6] Gög e, D., Sinapius, J.M., 2006. Exper iences with dynamic load simulation by means of modal f orces in the presence of structural non- linearities. Aerospace Science and T echnology 10, 411–419. [7] Kidg er, P ., Garcia, C., 2021. Equinox: neural netw orks in JAX via callable PyTrees and filtered transformations. Differentiable Programming wor kshop at Neural Information Processing Systems 2021 . [8] K ocijan, J., 2018. Parameter estimation of a nonlinear benchmark system. Science, Engineering & Education 1, 3–10. [9] Le venber g, K., 1944. A method for the solution of cert ain non-linear problems in least squares. Quar terly of Applied Mathematics 2, 164–168. [10] Mar quardt, D.W ., 1963. An algor ithm for least-squares estimation of nonlinear parameters. Jour nal of the Society for Industrial and Applied Mathematics 11, 431–441. [11] Masri, S., Bekey , G., Sassi, H., Caughey , T ., 1982a. Non-parametric identification of a class of nonlinear multideg ree dynamic systems. Earthquake Engineering & Str uctural Dynamics 10, 1–30. [12] Masri, S., Caughe y , T., 1979. A nonparametric identification tec hnique for nonlinear dynamic problems. Journal of Applied Mec hanics 46, 433–447. [13] Masri, S., Chassiakos, A., Caughey , T., 1993. Identification of nonlinear dynamic systems using neural networ ks. Jour nal of Applied Mechanics 60, 123–133. [14] Masri, S., Sassi, H., Caughey , T., 1982b. Nonparametric identification of nearly arbitrary nonlinear systems. Journal of Applied Mec hanics 49, 619–628. [15] N oël, J.P ., Kerschen, G., 2017. Nonlinear system identification in structural dynamics: 10 more y ears of progress. Mechanical Systems and Signal Processing 83, 2–35. [16] Pintelon, R., Schoukens, J., 2012. System Identification: A Frequency Domain Approach. 2nd ed., John Wiley & Sons, Hoboken, NJ. [17] Rader , J., L yons, T., Kidger, P ., 2024. Optimistix: modular optimisation in JAX and Equino x. arXiv:2402.09983 . [18] R ogers, T.J., Fr iis, T ., 2022. A latent restoring force approach to nonlinear system identification. Mechanical Systems and Signal Processing 180, 109426. [19] Saad, P ., Al Majid, A., Thouverez, F., Duf our, R., 2006. Equivalent rheological and restoring force models for predicting the har monic response of elastomer specimens. Journal of Sound and Vibration 290, 619–639. [20] V an Loan, C., 2003. Computing integ rals invol ving the matrix exponential. IEEE Transactions on Automatic Control 23, 395–404. [21] W igren, T., Schoukens, J., 2013. Three free data sets for development and benchmarking in nonlinear system identification, in: 2013 European Control Conf erence (ECC), IEEE. pp. 2933–2938. [22] W orden, K., 1990. Data processing and experiment design f or the restoring force sur face method, part I: integration and differentiation of measured time data. Mechanical Systems and Signal Processing 4, 295–319. Merijn Flo ren and Jan Swevers: Preprint submitt ed to Elsevier P age 20 of 20