Data-driven forced response analysis with min-max representations of nonlinear restoring forces

Data-driv en forced resp onse analysis with min-max represen tations of nonlinear restoring forces ∗ Akira Saito † Hirom u F ujita ‡ Abstract This pap er discusses a no v el data-driven nonlinearit y identification method for mechanical systems with nonlinear restoring forces suc h as polynomial, piecewise-linear, and general displacemen t-dep enden t nonlinearities. The prop osed metho d is built upon the universal appro ximation theorem that states that a nonlinear function can b e approximated by a linear combination of activ ation functions in artificial neural net work framework. The prop osed approach utilizes piecewise linear springs with initial gaps to act as the activ ation functions of the neurons of artificial neural netw orks. A library of piecewise linear springs with initial gaps are constructed, and the contributions of the springs on the nonlinear restoring force are determined by solving the linear regression problems. The piecewise linear springs are realized b y combinations of min and max functions with biases. The prop osed metho d is applied to a Duffing oscillator with cubic stiffness, and a piecewise linear oscillator with a gap and their nonlinearities are successfully determined from their free resp onses. The obtained mo dels are then used for conducting forced resp onse analysis and the results match well with those of the original system. The metho d is then applied to exp erimen tally-obtained free resp onse data of a cantilev ered plate that is sub jected to magnetic restoring force, and successfully finds the piecewise linear representation of the magnetic force. It is also shown that the obtained mo del is capable of accurately capturing the steady-state resp onse of the system sub ject to harmonic base excitation . Keyw ords: Nonlinear system, Piecewise linear system, Universal Approximation Theorem, ReLU neural net work 1 In tro duction In recent years, data-driven analysis metho ds hav e b een widely applied in linear and nonlinear mec hanical systems for predicting, designing and controlling their dynamics [ 1 ]. Ko opman-operator based methods, such as the extended dynamic mo de decomp osition (DMD) has been applied to nonlinear oscillators [ 2 ]. V arious data-driv en DMD-based metho ds ha ve also b een developed to date to predict the dynamics of nonlinear systems. F or instance, Galerkin pro jection metho d [ 3 , 4 ] has b een applied to successfully construct reduced order mo dels of nonlinear equations of motion. DMD with con trol [ 5 ] has also b een applied to predict forced resp onse of the nonlinear oscillator [ 6 ]. A gro wing num b er of attempts ha v e b een made to represen t nonlinear dynamics using artificial neural netw orks (ANNs). F or instance, Kuptsov et al. [ 7 ] inv estigated the equiv alence of a netw ork mo del consisting of a single hidden lay er of p erceptron with sigmoid functions and nonlinear functions in dynamical systems. They examined a couple of nonlinear ODEs including Lorenz system, and show ed that the netw ork mo dels exhibit almost equiv alen t dynamical b eha vior including bifurcation diagrams and F ourier sp ectra. They also successfully applied ANN to the dynamics of a stiff system [ 8 ]. Moukhliss et al. prop osed an ANN-based surrogate mo deling framework for predicting the frequency of nonlinear oscillations of functionally graded materials [ 9 ]. Luo et al. prop osed a hybrid physics-data-driv en metho d ∗ This man uscript is a preprint and has b een submitted to a journal for possible publication. † Departmen t of Mec hanical Engineering, Meiji Univ ersity , Ka w asaki, Kanagaw a 214-8571, Japan. E-mail: asaito@meiji.ac.jp ‡ Departmen t of Mechanical Engineering, Meiji Univ ersity , Ka wasaki, Kanagaw a 214-8571, Japan. 1 to iden tify unknown structural nonlinear boundary condition of systems by using both finite elemen t mo dels and m ulti-lay er p erceptron (MLP) mo del [ 10 ]. It has b een successfully applied to systems with v arious nonlinearit y t yp es including hardening spring, anisotropic spring, and hysteretic nonlinearit y . These ANN-based nonlinear identification metho ds are promising b ecause they are versatile and capable of represen ting a broad class of nonlinear systems. Ho wev er, the identified mo del tends to b ecome a black-box and hence it is hard to interpret its physical meanings. Raissi et al. dev elop ed a h ybrid data-driv en and mo del-based metho d that ac hiev es system iden tification and prediction of its dynamics for nonlinear partial differen tial equations called physics-informed neural net works [ 11 ]. This class of metho ds enable the solution of nonlinear partial differential equations (PDEs) b y minimizing a cost function that is a combination of the PDE constraint and data errors with resp ect to the ANN’s parameters, which result in the ANNs capable of not only solving the PDEs but also achieving their consistency with the measurement data. This metho d, how ev er, still result in the black-box nonlinear dynamics mo del based on the ANNs. Brun ton et al. dev elop ed a sparse identification of nonlinear dynamical systems (SINDy)[ 12 ], which utilizes sequen tial thresholded linear least squares of nonlinear terms by a library of nonlinear functions. This giv es relatively clear symbolic representation of the nonlinear terms, which tends to provide interpretable nonlinear dynamics mo dels. T o obtain in terpretable nonlinear dynamics mo del with ANNs, Gonzalez and Lara prop osed a nonlinear system iden tification metho d for second-order nonlinear ordinary differen tial equations (ODEs) based on the concept of characteristic curves where the system dynamics are represen ted in terms of the nonlinear elemen ts in the ODEs [ 13 ]. They combined the concept with SINDy or ANN, and the characteristic curve metho d has b een successfully applied to v an der Pol equation and mechanical systems with dry-friction. Building up on these studies, this pap er attempts to prop ose a simple data-driven system iden tification and forecasting method for mechanical systems with nonlinear restoring forces with ANN-like system represen tation combined with a SINDy-like linear regression, which gives us interpretable netw ork mo del of nonlinear restoring forces as a netw ork of piecewise linear springs . The prop osed approac h do es not require explicit expressions of the nonlinear forces as in [ 13 ]. F undamen tally , the prop osed metho d is based on a piecewise linear approximation of nonlinear functions applied to nonlinear dynamics. Piecewise linear appro ximation of nonlinear forces in the dynamics has long b een studied b y man y researc hers for v arious purp oses. F or instance, Garg prop osed a piecewise linear approximation of a certain class of nonlinear functions to simplify the numerical integration inv olv ed in obtaining the aero dynamic forces and moments from the nonlinear functions [ 14 ]. Cao et al. [ 15 ] prop osed a piecewise linear appro ximation of a smo oth discontin uous nonlinear force generated by a link ed pair of inclined springs that are pinned to rigid supp orts. They deriv ed an analytical expression of a trilinear approximation for the nonlinear force using t wo piecewise linear functions. Xie and Shih [ 16 ] also prop osed a metho d to solv e for the resp onse of a single degree of freedom (DOF) system with nonlinear stiffness, b y utilizing the piecewise linear approximation of nonlinear force and finding the exact closed form solution at each time step. More recently , Li et al [ 17 ] developed a metho d for designing a piecewise linear oscillator to mimic nonlinear restoring forces by selecting the breakp oin ts and slop es of the piecewise linear oscillators. It w as sho wn that the designed piecewise linear oscillator successfully approximates the b eha vior of a nonlinear system. As it will b e shown, the principle b ehind their metho d is shared with the prop osed approach in this pap er. Inspired by these studies, the prop osed approac h in this pap er integrates such piecewise linear approximations of nonlinear forces in to a single-la yer p erceptron (SLP)-like ANN mo del, where the con tributions of the springs are determined from data through SINDy-like linear regression. This pap er is organized as follows. In Section 2 , mathematical foundation of the prop osed approac h is presen ted. Numerical examples that sho w the v alidity of the prop osed metho d are presented in Section 3 including a Duffing oscillator and a single DOF system with a piecewise linear stiffness. In Section 4 , the prop osed metho d is applied to exp erimen tally-obtained nonlinear vibration data where a plate with a p ermanen t magnet at its tip is sub ject to nonlinear magnetic force. Conclusions of the paper are 2 Figure 1: Sc hematic diagram of a single lay er p erceptron summarized in Section 5 . 2 Mathematical foundation 2.1 Min-max represen tation of nonlinear restoring forces Let us consider the dynamics of a single mass-spring-damp er system of the form, m ¨ x + c ˙ x + k x + f ( x ) = 0 (1) where m is the mass, c is the damping co efficien t, k is the spring constan t, x is the displacemen t. Also, let f ( x ) b e a nonlinear function that represents an arbitrary nonlinear force that is a function of x . Based on the univ ersal approximation theorem [ 18 , 19 , 20 ], whic h is the foundation of the ANN, assuming that f ( x ) b elongs to a class of temp ered distributions, its functional approximation can b e found f ( x ) ≈ Z X j =1 η ( a j x + b j ) (2) where η are activ ation functions, a j are the weigh ts, and b j are the biases, and Z is the num b er of neurons in the p erceptron. The schematic diagram of the SLP is sho wn in Fig. 1 . The activ ation function can b e simple functions, suc h as sigmoid function, hyperb olic tangen t, and max function (or so called rectifier linear unit (ReLU) in ANN). The definition of the max function is max (0 , x ) = x for x > 0 or 0 for x ⩽ 0. Also, the min function can b e an activ ation function where min (0 , x ) = x for x < 0 or 0 for x ⩾ 0. Note that since min(0 , x ) = − max(0 , − x ), they can b e related by a linear transformation. If the min and max functions are c hosen as the activ ation function, η ( a j x + b j ) = max (0 , a j x + b j ) = a j max (0 , x + b j /a j ) = a j η ( x + d j ) where d j ≜ b j /a j . Similarly , η ( a j x + b j ) = min (0 , a j x + b j ) = a j min (0 , x + b j /a j ) = a j η ( x + d j ) where d j ≜ b j /a j . Defining a ≜ [ a 1 , a 2 , . . . , a Z ] T , d ≜ [ d 1 , d 2 , . . . , d Z ] T , and η ( x, d ) ≜ [ η ( x + d 1 ) , η ( x + d 2 ) , . . . , η ( x + d Z )] T , then Eq ( 2 ) is rewritten as f ( x ) ≈ a T η ( x, d ). With these notations, Eq ( 1 ) b ecomes m ¨ x + c ˙ x + k x + a T η ( x, d ) = 0 . (3) If b oth min and max functions are chosen as the activ ation functions, Eq. ( 2 ) can b e interpreted as a linear com bination of forces exerted b y piecewise linear springs with initial gaps, i.e., f ( x ) ≈ M X j =1 ¯ k j min (0 , x − ¯ g j ) + N X j =1 k j max (0 , x − g j ) (4) where ¯ k j and k j are the equiv alen t spring constan ts of the piecewise linear springs that corresp ond to a j of Eq. ( 2 ) for min and max functions resp ectively , ¯ g j and g j are the equiv alen t initial gaps that corresp ond to 3 (a) Original nonlinear system (b) Equiv alent system represen tation using a net work of piece- wise linear springs Figure 2: Sc hematic diagram of multiple piecewise springs acting on the mass ¯ k j and k j and b j = − ¯ g j or b j = − g j . It is noted that b oth the spring constants and the initial gaps can b e negativ e. Theoretically , either max or min function can approximate any nonlinear function. How ev er, to ac hieve b etter con vergence p erformance for b oth increasing and decreasing functions in x , b oth min and max functions are used in the approximation. The schematic of this concept is shown in Fig. 2 . Namely , the nonlinear force exerted by the nonlinear stiffness in the original system is replaced by a set of piecewise linear springs with differen t initial gaps and spring constan ts. This means that any nonlinear force that is dep enden t only on x can b e mo deled as a sum of piecewise linear stiffnesses whose spring constants are ¯ k j ’s and k j ’s with the corresp onding initial gaps ¯ g j ’s and g j . Namely , the following is equiv alen t to Eq. ( 1 ) for sufficien tly large M and N , m ¨ x + c ˙ x + k x + M X j =1 ¯ k j min (0 , x − ¯ g j ) + N X j =1 k j max (0 , x − g j ) = 0 . (5) More compactly , m ¨ x + c ˙ x + k x + ¯ k T ¯ ξ ( x, ¯ g ) + k T ξ ( x, g ) = 0 , (6) where ¯ ξ ( x, ¯ g ) ≜ [ min (0 , x − ¯ g 1 ) , · · · , min (0 , x − ¯ g M )] T , ¯ k = [ ¯ k 1 , · · · , ¯ k M ] T , ξ ( x, g ) ≜ [ max (0 , x − g 1 ) , · · · , max (0 , x − g N )] T , and k = [ k 1 , · · · , k N ] T . The authors hav e shown that Eq. ( 5 ) holds for a piecewise linear system with an initial gap [ 21 ], and the spring constants k for the max functions can b e successfully obtained from measurement data b y using SINDy . Namely , a piecewise linear system with an initial gap can b e appro ximated by a sum of piecewise linear systems with initial gaps. This pap er further explores this concept by examining the equiv alence of the system b eha vior of nonlinear systems with that of the appro ximated nonlinear system with Eq. ( 6 ). 2.2 Determination of equiv alen t spring constan ts and gaps In general ANN formation, the weigh ts and biases of the activ ation functions are determined based on the minimization of a loss function that is typically defined as a squared error b et ween the v alues of the functions and those of the approximations by optimization solvers suc h as sto c hastic gradient metho ds. Since the prop osed metho d uses a linear combination of nonlinear functions of a single lay er, simple linear regression of the nonlinear force with a list of nonlinear functions, or a library can b e applied as in SINDy . Tw o scenarios are considered. One is the case where the data of the nonlinear function of interest is a v ailable and direct fit of the nonlinear function is p ossible. Suc h cases are discussed in 2.2.1 and referred to as direct metho d. The other scenario is where the data of the nonlinear function of interest is not a v ailable but measurements data of x is av ailable. This case is discussed in 2.2.2 and referred to as indirect metho d. 4 2.2.1 Direct metho d Assuming that the oscillator mov es within a certain range, e.g., x L ⩽ x ⩽ x R , and the nonlinear restoring force f ( x ) is measured at finite and discrete measurement p oints ˜ x = [ x 1 , · · · , x n ] T , i.e., f ( ˜ x ) = [ f ( x 1 ) , · · · , f ( x n )] T where x L = x 1 and x R = x n and x 1 < · · · < x n . Then, the equiv alen t spring constan ts of min and max functions in Eq. ( 4 ) can b e obtained simply by solving the following linear least squares problem. Namely , we seek ˜ k and k such that minimize κ : ∥ L ( ˜ x , ¯ g , g ) κ − f ( ˜ x ) ∥ 2 , (7) where κ = [ ¯ k T , k T ] T con tains a list of unknown equiv alen t spring constants to b e determined. ¯ g and g con tain lists of M and N p ossible gap v alues that need to b e set a priori. L ( ˜ x , ¯ g , g ) contains a library of min and max functions ev aluated at x 1 , · · · , x n , i.e., L ( ˜ x , ¯ g , g ) = [ ¯ L ( ˜ x , ¯ g ) , L ( ˜ x , g )], where ¯ L ( ˜ x , ¯ g ) =    min(0 , x 1 − ¯ g 1 ) , · · · , min(0 , x 1 − ¯ g M ) . . . min(0 , x n − ¯ g 1 ) , · · · , min(0 , x n − ¯ g M )    (8) and L ( ˜ x , g ) =    max(0 , x 1 − g 1 ) , · · · , max(0 , x 1 − g N ) . . . max(0 , x n − g 1 ) , · · · , max(0 , x n − g M )    . (9) The minimizer of Eq. ( 7 ) can b e obtained as κ = [ L ( ˜ x , ¯ g , g )] † f ( ˜ x ) , (10) where † denotes the Mo ore-P enrose pseudo-in verse of a rectangular matrix. Obviously , this pro cess of obtaining the contributions from the min and max functions do es not require sophisticated gradien t-based optimization metho ds such as sto c hastic gradient descent or Adam [ 22 ] that are used in the construction of mo dern ANNs. This comes at the exp ense of forming a list of p ossible gaps ( ¯ g i ’s and g i ’s) that corresp ond to the biases of the neurons in ANNs, which are typically obtained by the optimization solvers. Instead, the contributions from the piecewise linear springs that are in the list ( ¯ L and L ) are obtained as the spring constan ts ( κ ), which corresp ond to the weigh ts of the activ ation functions of the neurons, b y the regression Eq. ( 10 ) . As long as the moving range of the oscillator can b e estimated a priori, the list of gaps can b e created relatively easily from the moving range. Unnecessary gap v alues are then naturally eliminated in the pro cess of linear regression. 2.2.2 Indirect metho d This case also assumes that the oscillator mov es within a certain range x L ⩽ x ⩽ x R . This time, it is assumed that the direct measuremen t of nonlinear force is not feasible and hence its data is unav ailable. Instead, displacement, velocity , and acceleration of the nonlinear system are measured at discrete time instan ts t = [ t 1 , t 2 , · · · t n ], and they are used to indirectly iden tify the nonlinear restoring force. F or this purp ose, regression of time series of the a v ailable data is conducted as follo ws. Rearranging the terms in Eq. ( 6 ), we obtain, ¨ x + 2 ζ ω n ˙ x + ω 2 n x = − ¯ k ′ T ¯ ξ ( x, ¯ g ) − k ′ T ξ ( x, g ) , (11) where ζ = c/ 2 √ mk , ω n = p k /m , ¯ k ′ = ¯ k /m , and k ′ = k /m . Note that ω n is the natural frequency and ζ is the damping ratio with the absence of the nonlinear terms on the right hand side of Eq. ( 11 ). Let ˜ a = [ ¨ x ( t 1 ) , · · · , ¨ x ( t n )] T = [ a 1 , · · · , a n ] T , ˜ v = [ ˙ x ( t 1 ) , · · · , ˙ x ( t n )] T = [ v 1 , · · · , v n ] T , ˜ x = [ x ( t 1 ) , · · · , x ( t n )] T = [ x 1 , · · · , x n ] T . 5 Figure 3: Con vergence of the approximated function with resp ect to the n umber of functions. p 2 = 10. The regression was conducted for data in the filled region ( − 5 ⩽ x ⩽ 5). W e then seek ¯ k and k such that minimize κ ′ : ∥L ( ˜ x , ¯ g , g ) κ ′ − ˜ a ℓ ∥ 2 , (12) where κ ′ = [ ¯ k ′ T , k ′ T ] T , and ˜ a ℓ corresp onds to the linear p ortion of the equation of motion ev aluated with the resp onses of the nonlinear system, or ˜ a ℓ ≜ ˜ a + 2 ζ ω n ˜ v + ω 2 n ˜ x . Therefore, the minimizer of Eq. ( 12 ) is obtained as, κ ′ = [ L ( ˜ x , ¯ g , g )] † ¯ a ℓ . (13) With the obtained κ ′ , the right hand side of Eq. ( 11 ) can b e ev aluated for an y x . 3 Numerical examples 3.1 Duffing oscillator The prop osed metho d is applied to a single DOF system with a cubic restoring force, or a Duffing oscillator of the form: m ¨ x + c ˙ x + p 1 x + p 2 x 3 = 0 . (14) where p 1 and p 2 are arbitrarily chosen stiffness parameters. First, the equiv alen t form of the nonlinear term p 2 x 3 represen ted by Eq. ( 4 ) needs to b e found. Namely , the spring constants k j need to b e computed, or they need to b e trained using known data of f ( x ). In this case, the spring constants can b e obtained by simple regression using the direct metho d as describ ed in 2.2.1 . 3.1.1 Piecewise linear appro ximation of the cubic term Figure 3 sho ws the results of the approximations using the min-max functions for ( M , N ) = (8 , 8) , (32 , 32), and (128 , 128) where p 2 = 10. The spring constants ¯ k and k w ere obtained for the v alues in − 5 ⩽ x ⩽ 5. As can b e seen, the approximation is p oor for ( M , N ) = (8 , 8). Ho wev er, as M and N increase, the appro ximation b ecomes more accurate. Indeed, when ( M , N ) = (128 , 128), as can b e seen, the original function is quite well represented b y the piecewise linear functions esp ecially within − 5 ⩽ x ⩽ 5. The appro ximation quality is degraded outside the range of − 5 ⩽ x ⩽ 5, but it is exp ected b ecause the spring constan ts ¯ k and k were trained by the dataset in − 5 ⩽ x ⩽ 5 and the capability of the appro ximation to extrap olate the function outside the range is not guaranteed. 6 (a) ( M , N ) = (8 , 8) (b) ( M , N ) = (32 , 32) (c) ( M , N ) = (128 , 128) Figure 4: Spring constants versus the initial gap v alues Figure 5: Time history of the oscillator with the original cubic nonlinearit y and with the min-max represen tations In addition, Fig. 4 sho ws the obtained equiv alent spring constan ts versus the gap v alues for different v alues of M and N . As seen in Fig. 4 (a), the springs constan ts are large for large initial gaps near the lo wer and upp er b ound of the range. This makes sense b ecause the cubic nonlinearity pro duces the largest v alues at the b ounds and hence the springs with large spring constan ts with large initial gaps need to b e added to represent suc h a nonlinear force. 3.1.2 F ree resp onse The original equation of motion Eq. ( 14 ) , and the ones with the nonlinear term replaced with the min-max represen tations with the trained ¯ k and k hav e b een solved by time integration for the initial condition of x (0) = − 4 . 0, and ˙ x (0) = 0. The time histories are sho wn in Fig. 5 . As can b e seen, the predicted time history of x ( t ) with all approximation cases matc h well with the original data until approximately t = 5s. Ho wev er, prediction accuracy gradually decreases, and the delay from the original data b ecomes visible, esp ecially for ( M , N ) = (8 , 8). On the other hand, the one with ( M , N ) = (32 , 32) do es not sho w significant dela y until around t = 25s. The approximation with ( M , N ) = (128 , 128) shows almost negligible dela y 7 (a) With original cubic stiffness (b) With ( M , N ) = (8 , 8) (c) With ( M , N ) = (128 , 128) Figure 6: Comparison of the time-frequency relationships using scalogram of the resp onse T able 1: RMSE v alues for the approximations ( M , N ) (8 , 8) (32 , 32) (128 , 128) RMSE (m) 2 . 38 0 . 497 0 . 215 b efore t = 30s. T able. 1 sho ws the ro ot mean square error (RMSE) v alues b et ween the approximations and the original resp onse. The RMSE v alue decreases as ( M , N ) increases. This also supp orts the con vergence of the prop osed approach with resp ect to M and N . T o see the frequency con tents in the resp onse more clearly , wa v elet transform has b een applied to the original resp onse, the ones with ( M , N ) = (8 , 8), and ( M , N ) = (128 , 128). The obtained scalograms are shown in Fig. 6 . Ov erall, the trend of the resp onse of the original system can b e well captured b y the approximations. Ho w ever, for the case with ( M , N ) = (8 , 8), the dominant frequency comp onen ts near t = 30s b ecome different from those of the original system. On the other hand, we can see that the frequency con tents can b e well captured by the approximation with ( M , N ) = (128 , 128). 3.1.3 F orced resp onse T o examine the applicability of the prop osed approach to forced resp onse problems, a forcing term has b een added to Eq. ( 14 ), i.e., m ¨ x + c ˙ x + p 1 x + p 2 x 3 = F sin(2 π f t ) , (15) where F is the amplitude of forcing, f is the excitation frequency . Also, the same forcing term is added to the equiv alen t equation of motion with the cubic term replaced by its equiv alen t min-max representation, i.e., m ¨ x + c ˙ x + p 1 x + ¯ k T ¯ ξ ( x, ¯ g ) + k T ξ ( x, g ) = F sin(2 π f t ) . (16) The steady-state solutions of the equations Eqs. ( 15 ) and ( 16 ) ha ve b een sought for v arious frequencies o ver 0 ⩽ f ⩽ 5 H z b y using time integration for 200 cycles of a p eriod of excitation. The results are sho wn 8 Backbone curve Cubic stif fness Approximate solution (M,N)=(128,128) (M,N)=(32,32) (M,N)=(8,8) Figure 7: Comparison of the forced resp onse results with single cubic stiffness and multiple PWL springs, along with classical approximate solutions and backbone curve in Fig. 7 for the steady-state resp onses of Eqs. ( 15 ) and ( 16 ) for ( M , N ) = (8 , 8), (32 , 32), and (128 , 128). When solving the equations for the steady-state, the excitation frequency f w as swept up ward or do wnw ard to obtain solutions of m ultiple stable branc hes. Backbone curve for the undamp ed system ( c = 0) and the appro ximate solutions b y a classical p erturbation metho d [ 23 ] ha ve also b een plotted in Fig. 7 . Note that f is normalized with resp ect to the linear natural frequency f n = p k /m/ (2 π ), and the amplitude of oscillation is normalized with resp ect to the static equilibrium X st that satisfies p 1 X st + p 2 X 3 st − F = 0. Note that it w as numerically obtained as X st = 0 . 3930. As can b e seen, the resp onse curve of the original system with the cubic stiffness shows a typical stiffening effect where the resp onse curv e is b en t tow ard high frequency and shows the primary resonance. Also, there are secondary resonances b elo w f /f n = 1 including 3:2 sup erharmonic resonance. W e can now see that these characteristics are well captured by the appro ximate system with ( M , N ) = (128 , 128), including the primary and the secondary resonances. On the other hand, for the case with ( M , N ) = (32 , 32), the v alues tend to b e sligh tly off from those of the original system or the appro ximate system with ( M , N ) = (128 , 128). F urthermore, if the n umbers of min-max functions are further reduced down to ( M , N ) = (8 , 8), the system fails to predict the accurate steady-state resp onse for almost the entire frequency range considered. Interestingly , it sho ws a strongly nonlinear yet typical resp onse curv e of piecewise linear systems (see, e.g., [ 24 ]) at around X/X st = 3 . 2, i.e., X ≈ 1 . 25. This approximately corresp onds to the smallest gap v alue of the activ e min and max functions, i.e., the system acts almost as a piecewise linear system with the smallest gap v alue, which do es not w ell represen t the Duffing oscillator. 3.2 Piecewise linear system Next, a single DOF system with a piecewise linear stiffness sub jected to harmonic forcing is considered. Suc h a piecewise linear system has b een studied extensively due to its simple yet strong nonlinearit y [ 25 ]. The sc hematic diagram of the system is sho wn in Fig. 8 . The piecewise linear spring is fixed at the right end and the initial gap b et w een the spring and the mass is L . The equation of motion of the system is written as, m ¨ x + c ˙ x + p 1 x + p 2 max(0 , x − L ) = F sin(2 π f e t ) , (17) where p 2 is the spring constant of the piecewise linear stiffness. The piecewise linear spring with gap L is appro ximated using the min and max functions with different gap v alues (without L ) using the pro cedure 9 AAACoHichVG7SgNBFD1Z3/EVtRHShARFmzAromIVEMROoyaKRsLuOomD+2J3EojBwtYfsLBSEBEbrfwAG3/AIp8glgo2Ft5sFkRFvcPMnDlzz505XN01hS8Za0SUtvaOzq7unmhvX//AYGxoOO87Fc/gOcMxHW9T13xuCpvnpJAm33Q9rlm6yTf0/YXm/UaVe75w7HVZc/mOpZVtURKGJokqxuKLBV/YE1MFVyRKCTlZ2PNdzeB1lVuHxViKpVkQiZ9ADUEKYaw4sTsUsAsHBiqwwGFDEjahwaexDRUMLnE7qBPnERLBPcchoqStUBanDI3YfVrLdNoOWZvOzZp+oDboFZOmR8oExtgju2Iv7IFdsyf2/mutelCj+Zca7XpLy93i4PHo2tu/Kot2ib1P1R8KnbL/9iRRwlzgRZA3N2CaLo1W/erBycva/OpYfZyds2fyd8Ya7J4c2tVX4yLLV08RpQap39vxE+Sn0upMWs1OpzLJsFXdiCOJCerHLDJYwgpy9O4RLnGDWyWpLCnLSraVqkRCzQi+hLL1AdWSmhs= F sin(2 ω ft ) Figure 8: Schematic diagram of the piecewise linear system -10 -5 0 5 10 x -20 0 20 40 60 80 100 f ( x ) p 2 max(0, x - L ) (M,N)=(256,256) (M,N)=(32,32) (M,N)=(8,8) Figure 9: Con vergence of the approximated function with resp ect to the n umber of functions. p 2 = 10. The regression was conducted for data in the filled region ( − 5 ⩽ x ⩽ 5). describ ed in 2.2 . Note that the equiv alence of a single max function with a gap and multiple max functions with different gap lengths w as considered in Ref. [ 21 ] b oth n umerically and exp erimen tally , and free resp onse was well captured by the max functions. With the assumption that the data of nonlinear force is av ailable, direct metho d describ ed in 2.2.1 can b e applied. F or this numerical exp erimen t, L = 0 . 5 is assumed. The linear regression was conducted using data taken within − 5 ⩽ x ⩽ 5 for ( M , N ) = (8 , 8), (32 , 32), and (256 , 256). The approximated functions as well as the original piecewise linear stiffness are sho wn in Fig. 9 . As can b e seen, the appro ximation is p o or for ( M , N ) = (8 , 8) ev en within − 5 ⩽ x ⩽ 5. On the other hand, for ( M , N ) = (32 , 32), the min and max functions well represent the original function. With ( M , N ) = (256 , 256), the difference is negligible for − 5 ⩽ x ⩽ 5. The obtained spring constan ts for the corresp onding gap v alues are plotted in Fig. 10 . As can b e seen, ev en though the original function consists only of the max function, b oth min and max functions contribute to the representation of a single max function. As the num b er of min and max functions increases, we can see that the contributions of the min and max functions whose gap v alues are close to the gap v alue of the original max function ( L ) b ecome dominant. This app ears as a spik e in the plot as shown in Fig. 10 (c). This c haracteristics can b e exploited to identify unknown L , as it has already b een discussed in Ref. [ 21 ]. 3.3 F ree resp onse F ree resp onse analyses hav e b een conducted for the original piecewise linear system and the system with the multiple min and max functions for x (0) = − 2 and ˙ x (0) = 0. Namely , F = 0 in Eq. ( 17 ) . Other parameter v alues are set to m = 1, c = 0 . 1, and k = 1. The results are sho wn in Fig. 11 as phase p ortrait format to clearly visualize the changes in the dynamics as the system stiffness changes. The vertical dashed 10 (a) ( M , N ) = (8 , 8) (b) ( M , N ) = (32 , 32) (c) ( M , N ) = (256 , 256) Figure 10: Spring constants versus the initial gap v alues for the piecewise linear stiffness ( L = 0 . 5) T able 2: RMSE v alues for the approximations ( M , N ) (8 , 8) (32 , 32) (256 , 256) RMSE (m) 0 . 554 0 . 229 0 . 028 line designates L = 0 . 5 where the mass collides with the piecewise linear stiffness. The phase p ortrait with ( M , N ) = (8 , 8) p o orly predicts the correct tra jectory . The one with ( M , N ) = (32 , 32) predicts that of the original piecewise linear system well, but there are noticeable differences esp ecially near the gap. On the other hand, ( M , N ) = (256 , 256) predicts quite well the one with the original system even near the gap region. T able 2 shows the RMSE v alues for the appro ximations, where the reduction in RMSE can b e observ ed for increasing M and N . This indicates the high accuracy of the prop osed system representation of the nonlinear function even when the nonlinearity is piecewise linear. 3.4 F orced resp onse F orced resp onse analyses hav e b een conducted for 0 ⩽ f e ⩽ 0 . 5, or 0 ⩽ f e /f n ⩽ π where f n = 1 / 2 π . At eac h frequency , the equation of motions w ere solv ed n umerically for x (0) = 0 and ˙ x (0) = 0 until they reach the steady-state solutions. The amplitude of the resp onse w as computed and normalized with resp ect to the static displacement. It is plotted for the normalized excitation frequency and shown as Fig. 12 for the original piecewise linear system, the approximate systems with ( M , N ) = (8 , 8), (32 , 32), and (256 , 256). As can b e seen, there is a primary resonance at f e /f n = 1 . 47, whereas the largest secondary resonance can b e seen at f e /f n = 0 . 622. The prop osed mo del with ( M , N ) = (8 , 8) fails to predict b oth primary and secondary resonances. The prop osed mo del with ( M , N ) = (32 , 32) succeeded in predicting the primary resonance, ho wev er, it failed to capture the secondary resonance and resp onses at higher frequency range than the primary resonance. The prop osed mo del with ( M , N ) = (256 , 256) captures b oth the primary and the secondary resonances. In addition to the forced resp onse results, backbone curve is shown in Fig. 12 . The backbone curve was extracted from the instantaneous frequency and amplitude of the free resp onse of the system with the identified multiple PWL stiffnesses using the initial condition of x (0) = 2 . 0 and ˙ x (0) = 0. The instan taneous amplitude of the resp onse was defined as the mean v alue of the upp er and lo wer env elop e of the resp onse. The instantaneous frequency of the resp onse was computed from the 11 Figure 11: Comparison of the free resp onse results with a single piecewise linear stiffness and multiple PWL stiffnesses. time instants at which the displacemen t crossed the instan taneous center of oscillation. The instantaneous cen ter of oscillation was defined as the midp oint b et w een the instantaneous upp er and low er env elop es. As seen, the bac kb one curv e of the obtained system captures the amplitude dep endence of the resonan t frequency quite well. T o further explore the v alidit y of the prop osed approac h, tra jectory of the solutions as w ell as their F ourier sp ectra w ere examined. Figure 13 shows the phase p ortraits and the F ourier sp ectra computed b y F ast F ourier T ransform (FFT). The horizon tal axis of the phase p ortrait is normalized with resp ect to the static displacement x st and the vertical axis is normalized with resp ect to v n = ω n x st where ω n = p k /m = 1. The horizontal axes of the FFT sp ectra are normalized with resp ect to the frequency of the external forcing f e , whereas the vertical axes are normalized with resp ect to x st . At the primary resonance shown in Fig. 13 (a), the resp onse is dominated by the harmonic comp onen t of the excitation frequency . Other frequency comp onen ts that are integer m ultiples of the excitation frequency are muc h smaller than the one with f /f e = 1. On the other hand, at the secondary resonance, the frequency comp onen t with twice the excitation frequency is as large as that with the excitation frequency , as can b e seen in Fig. 13 (b). This indicates the existence of the sup erharmonic motion. Indeed, the lo op in the phase p ortrait has a knot, indicating the sup erharmonic motion at t wice the excitation frequency . This arises b ecause the excitation frequency approximately equals to twice the natural frequency of the system, resulting in 2:1 sup erharmonic resonance. 4 Exp erimen tal v alidation In this section, the prop osed approac h is applied to a cantilev ered plate with magnetic-force induced nonlinearit y from exp erimen tally obtained dataset and its applicability to real measurement data is discussed. First, the nonlinear restoring force acting on the system is identified by the prop osed approac h from free resp onse data. Second, forced resp onse of the system under base excitation is predicted by adding the excitation term in the obtained equation of motion and its v alidit y is discussed by comparing it with the measurement. 12 Original Figure 12: Comparison of the forced resp onse results with a single piecewise linear stiffness and m ultiple PWL springs (a) Primary resonance: f e /f n = 1 . 470 (b) Secondary resonance: f e /f n = 0 . 622 Figure 13: Phase p ortraits and F ourier sp ectra at primary and secondary resonances 4.1 Exp erimen tal setup The system studied in this section is sho wn in Fig. 14 . A can tilevered plate with a p ermanent magnet at its tip is clamp ed and fixed to a surface plate. The plate is made from Titanium Alloy and its dimension is 175 mm × 50 mm × 0 . 5 mm . There is another p ermanen t magnet that is fixed to an aluminum blo c k rigidly attac hed to the surface plate, and is placed near the tip of the cantilev ered plate. The distance b et ween the magnet attached to the cantilev ered plate and the other magnet was set to 12 mm. The directions of the p olarit y of the magnets are b oth in the direction p erp endicular to the surface of the plate, but in the opp osite directions with each other. The magnet attac hed to the cantilev ered plate exp eriences repulsiv e magnetic force if its p osition is not aligned with the fixed magnet, as shown in Fig. 14 (b). Since the plate also exp eriences the restoring force due to the elasticity of the plate, there is a stable equilibrium p osition where the restoring elastic force balances the repulsive force. The same applies to the situation when the plate is b en t tow ard the opp osite side. Therefore, this oscillator has t wo stable equilibria and hence b ecomes a bistable oscillator. This class of oscillator has b een used as a type of energy harvesting devices [ 26 ] and nonlinear energy sinks [ 27 ]. The exp ected single DOF approximation of the equation of motion of the system is given as ¨ x + 2 ζ ω n ˙ x + ω 2 n x = f t ( x ) (18) 13 Cantilevered beam Magnet (fixed) Magnet Clamped Magnet (fixed) Clamped (a) Side view (b) Overhead view Figure 14: T est equipment for free vibration test. The plate is in equilibrium state. where x is the displacemen t of the plate at the tip, f t ( x ) is the nonlinear restoring force that stem from the magnetic repulsiv e force generated b etw een the fixed magnet and the magnet at the tip of the plate. Note that the equiv alen t stiffness co efficien t k mainly represen ts the linear restoring force that comes from the elasticit y of the plate. W e consider the identification of such nonlinear restoring force exerted by the magnets, and aim to dev elop mathematical mo del that can well represent the dynamics of the system from measured data. In this case, the magnetic nonlinear force is not directly measured. Instead, displacement is measured by laser displacemen t sensor (LDS), and the indirect metho d introduced in Section 2.2.2 is applied. 4.2 P oten tial constraints on the basis functions F or the identification of the magnetic force, the follo wing constraints on the min and max functions are considered. T he magnetic force acting b et ween the magnets is dep enden t on the angle, alignment, and the relativ e distance of the magnets. In the present setup, the amplitude of vibration is sufficiently small so that the changes in the angle and the alignmen t are negligible. Therefore, it is reasonable to assume that the magnetic force is dep endent solely on the relative distance of the magnets. Since one of the magnets is fixed to the base, the relative distance b et w een the magnets is dep enden t only on the displacement of the magnet that is attached to the tip of the cantilev er plate. As a consequence, the magnetic force can b e regarded as a conserv ativ e force, and asso ciated p oten tial function can b e defined. Namely , we h yp othesize that the p otential can b e approximated as a linear combination of a linear function, quadratic function, and piecewise quadratic functions of the form: V t ( x, g ) = V c ( x ) + V ℓ ( x ) + V p ( x, g ) (19) where V c ( x ) = q 1 x , V ℓ ( x ) = 1 2 q 2 x 2 , V p ( x, g ) = κ T ϕ ( x, g ), κ = [ κ 1 , . . . , κ N ] T , ϕ ( x, g ) = [ ϕ ( x, g 1 ) , . . . , ϕ ( x, g N )], and ϕ ( x, g i ) ≜ 1 2 x 2 + 1 2 min(0 , x + g i ) 2 + 1 2 max(0 , x − g i ) 2 . (20) q 1 , q 2 , and κ ′ i s are the contributions from these terms determined b y the regression. The nonlinear force is then obtained as: f t ( x, g ) = − d V t d x ( x, g ) = f c ( x ) + f ℓ ( x ) + f p ( x, g ) (21) 14 (a) F ree response of the plate tip (b) Bac kb one curv e obtained from the free resp onse Figure 15: Time resp onse of the displacement of the plate tip and bac kb one curve obtained from the free resp onse where f c ( x ) = − q 1 , f ℓ ( x ) = − q 2 x , f p = − κ T f ( x, g ), f ( x, g ) = [ f ( x, g 1 ) , . . . , f ( x, g N )] T , and f i ( x, g i ) = − d d x ϕ ( x, g i ) = − { x + min(0 , x + g i ) + max(0 , x − g i ) } . (22) The expression on the right hand side is defined as, ψ ( x, g i ) ≜ d d x ϕ ( x, g i ) = x + min(0 , x + g i ) + max(0 , x − g i ) (23) and ψ is herein referred to as basis functions and used for the formation of the library . With these assumptions, the physics b ehind the data can b e enforced as constrain ts on the regression, i.e., p oten tial function asso ciated with the iden tified force can b e defined, otherwise the min-max functions do not resp ect the ph ysics and hence ma y pro duce mo dels that do not follo w physical laws. 4.3 F ree resp onse With the system shown in Fig. 14 , free resp onse exp erimen t has b een conducted. The pro cedure of the exp erimen t is stated as follows. The tip of the plate was displaced along transverse direction. The tip was then released with zero initial v elo cit y . The displacement of the plate tip where the magnet is attac hed was measured b y an LDS (IL-100, KEYENCE, Japan). The sampling rate was 1kHz. T o obtain the velocity , deriv ative of the obtained discrete time history of the displacement was tak en n umerically . The deriv ativ e of the velocity was then tak en to obtain the acceleration. Note that a low-pass filter with cut-off frequency of 100Hz has b een applied to the measured displacemen t to eliminate excessive noise in the deriv atives of the displacemen t. Result of the free resp onse measurement is sho wn in Fig. 15 (a). The plate tip vibrates with appro xi- mately 13Hz. The amplitude of the displacement gradually decrease due to damping. As a result, the frequency of oscillation slightly deviates as the amplitude of oscillation decreases. T o further examine this, the backbone curve was extracted from the data of Fig. 15 (a) by the instantaneous frequency and amplitude, and is shown in Fig. 15 (b). As can b e seen, clear softening effect can b e observed. Namely , the bac kb one curv e is b ent tow ard left, which means that the resonan t frequency decreases as the amplitude of oscillation increases, or the resonant frequency increases as the amplitude of oscillation decreases. This mak es sense b ecause if the tip displacement amplitude increases, the magnetic repulsive force w eakens b ecause the distance b et ween the p ole of the magnet at the tip and that of the other magnet fixed at the base increases. 15 (a) Poten tial (b) Magnetic force Figure 16: Identified p oten tial and its asso ciated magnetic force With the free resp onse shown in Fig. 15 , the prop osed metho d has b een applied. First, the damping ratio ζ and the natural frequency ω n corresp onding to m , c and k in Eq. ( 18 ) w ere iden tified from a free resp onse of the same plate without the magnetic force, i.e., fixed magnet was remov ed from the fixture. The identified v alues w ere ω n = 65 . 2 rad / s , which corresp ond to 10.4Hz, and ζ = 0 . 0054. Second, the indirect metho d has b een applied on the free resp onse of the system with the magnetic force. The library used for the regression of the nonlinear force consists of a constan t term and the piecewise linear forces defined in Eq. ( 21 ), i.e., L ( ˜ x , g ) =    1 , ψ ( x 1 , g 1 ) , . . . , ψ ( x 1 , g M ) . . . 1 , ψ ( x n , g 1 ) , . . . , ψ ( x n , g M )    . (24) where g 1 , . . . , g M w ere equally sampled M p oin ts within − 0 . 02 ⩽ x ⩽ 0 . 02. The iden tified p otential function and the asso ciated magnetic force are sho wn in Fig. 16 . In Fig. 16 (a), the obtained total p oten tial energy V t ( x ) is shown along with the constituent p oten tial energies with resp ect to displacement x , where V c ( x ) is the p oten tial corresp onding to the constant force, V ℓ ( x ) is the p oten tial corresp onding to linear force, ϕ i ( x )’s are the piecewise quadratic p oten tial defined as Eq. ( 19 ) . As seen in Fig. 16 , V t ( x ) is a con vex function with tw o inflection p oin ts that are almost symmetric with resp ect to the v ertical axis. In Fig. 16 (b), the obtained total magnetic force f t ( x ) is sho wn with the constan t term f c ( x ), linear term f ℓ ( x ), and piecewise linear terms f i ( x ). As can b e seen in the figure, the magnetic force is almost symmetric with resp ect to the origin. The slop e of f t is slightly steep er than the linear p ortion of the force f ℓ due to the con tributions from quadratic and piecewise quadratic terms. A t the inflection p oints of the p oten tial V t ( x ), f t ( x ) tak es its extrema. Figure 17 shows the identified spring constants, or the contributions of the basis functions ψ ( x, g ) with resp ect to g i for M = 20, 32, and 40. As M increases, the p ositiv e contribution near g i = 0 . 005 tends to b ecome the largest whereas the negative contribution near g i = − 0 . 005 tends to b ecome the smallest. This generally implies that the effects of the piecewise linear springs tend to b ecome larger as the distance from the origin increases. The springs with p ositiv e gaps tend to ha ve p ositiv e spring constants whereas those with negativ e gaps tend to hav e negative spring constan ts. The combined effect of these springs represen t the observed softening nature of the nonlinearit y . The time histories of the measured displacemen t, which is the same data as Fig. 15 (a), velocity , and 16 g i -0.02 -0.01 0 0.01 0.02 k i 0 2000 4000 (a) M = 20 g i -0.02 -0.01 0 0.01 0.02 k i 0 2000 4000 (b) M = 32 g i -0.02 -0.01 0 0.01 0.02 k i 0 2000 4000 (c) M = 40 Figure 17: Spring constants versus the initial gap v alues for the min-max functions acceleration are shown in Fig. 18 for 0s ⩽ t ⩽ 5s. The data for 0s ⩽ t ⩽ 3s was used for obtaining magnetic restoring force shown in Fig. 16 (b) by the regression, and the resulting nonlinear gov erning equations w ere solved by n umerical time integration to forecast the time histories of displacement, velocity , and acceleration for 3s ⩽ t ⩽ 5s. As seen, the forecasted displacement, v elo cit y , and acceleration agree very w ell with the measurements. T o further examine the characteristics of the obtained mo del, con vergence study with resp ect to the n umber of ϕ ( x, g i ), or M , has b een conducted. The RMSE v alue b et ween the measured acceleration and the fit result was computed for 0 ⩽ t ⩽ 3. Also, the RMSE b et ween the measured and the predicted v alues for 3 ⩽ t ⩽ 5 w as computed. This means that the data 0 ⩽ t ⩽ 3 w as used for training the mo del, and that for 3 ⩽ t ⩽ 5 was used for v alidation. The results are shown in Fig. 19 . The RMSE v alue for the fit decreases monotonically with resp ect to M . On the other hand, the RMSE for the prediction decreases as M increases up to around 30. F or M ⩾ 30, the RMSE for the prediction increases. Considering that the RMSE for the fit decreases for M ⩾ 30, this clearly shows the consequence of o verfitting. This shows that con vergence study needs to b e conducted with resp ect to M , or cross-v alidation needs to b e conducted to find the optimal M . F urthermore, con vergence study with resp ect to the duration of the fit to the total time has also b een conducted. Fit duration was v aried from 0.5s to 4.5s while the total time w as fixed to 5.0s. Namely , the fit duration w as v aried from 9.98% to 90.0%. M w as fixed to 32. The prediction RMSE decreases as fit duration increases up to 70%. It then slightly increases for fit duration larger than 70%. This trend implies that ev en though the prediction accuracy increases as fit duration increases for mo derately small fit duration, if the fit duration exceeds 70%, it suffers from ov erfitting. This also sho ws that the conv ergence study with resp ect to the prediction accuracy should also b e conducted to obtain the optimal fit duration. 17 Figure 18: Results of free resp onse ( M = 32). The regression was conducted for data in the filled region (0 ⩽ t ⩽ 3). 4.4 F orced resp onse Next, forced resp onse of the oscillator is examined. Namely , the prediction capability of the mo del obtained in the previous section is examined under harmonic excitation. The obtained results are compared with measuremen ts. The forced resp onse test w as conducted by using an electro dynamic shaker. The test setup is shown in Fig. 21 . The oscillator w as fixed at the moving head of the electro dynamic shaker (K2004E01, The Mo dal shop, USA). A sinusoidal signal was supplied from a function generator (WF1947, NF Corp oration, Japan) to the electro dynamic shak er to generate the base excitation to the oscillator. The displacement was measured at the mo ving head of the electro dynamic shaker where the plate is attac hed (herein referred to as base) and at the tip of the plate (herein referred to as tip) using tw o of the LDSs. With this setup, forced resp onse tests ha ve b een conducted. The sinusoidal excitation with a fixed frequency has b een applied to the plate, and its vibration resp onse w as measured at the tip and the base until it reaches the steady-state. The measurement w as rep eated for different frequencies from 5Hz to 30Hz with 1Hz increment. T o examine the amplitude-dep endence of the system, the shaker was op erated with tw o vibration amplitude lev els, i.e., 0.27mm and 0.83mm, which approximately corresp ond to the input DC voltages to the shaker of 0.1V and 0.3V, resp ectiv ely . As a representativ e result, Fig. 22 (a) shows the results of forced resp onse at the steady-state for the sinusoidal excitation with its frequency of 15Hz. The displacements measured at the tip and the base are shown. The resp onse was extracted from the measurements for an arbitrary time window of 0.6s in the steady-state. As can b e seen, the sin usoidal motion of the base induces almost sin usoidal motion of the tip with phase shift. The phase shift b et w een the base and the tip motions were computed to b e 178.20deg. F orced resp onse calculation has b een conducted with the mo del identified from the free resp onse. The 18 M 0 10 20 30 40 RMSE(m/s 2 ) 2 3 4 5 6 7 8 9 Fit Prediction Figure 19: Conv ergence study with resp ect to the num ber of ϕ ( x, g i ) (Fit duration=50%) Fit d uration (% o f t otal) 0 20 40 60 80 100 PredictionRMSE(m/s 2 ) 2 4 6 8 10 12 14 16 Figure 20: Conv ergence study with resp ect to the fit duration to total duration ( M = 32) equation of motion Eq. ( 18 ), which is unforced, was conv erted to that with the base excitation, i.e., ¨ x + 2 ζ ω n ˙ y + ω 2 n y = f t ( y , g ) (25) with y ≜ x − X b sin(2 π f e t ) (26) where y is the relativ e displacement of the oscillator measured from the base, X b is the amplitude of the base excitation, f e is the excitation frequency of the base. F or the ev aluation of the nonlinear force f t ( y , g ), the parameters q 1 , q 2 , and κ in Eq ( 21 ) obtained from the free resp onse analysis were used. Figure 22 (b) shows the steady-state resp onse predicted b y the mo del generated based on the free- resp onse. The resp onse was extracted from the results of numerical integration for an arbitrary time windo w of 0.6s in the steady-state. As can b e seen, the tip displacemen t prediction shown in Fig. 22 (b) agree well with the measurement shown in Fig. 22 (a), in terms of the amplitude of oscillation and the phase shift from the base displacemen t. The phase shift was computed to b e 175.86deg, which is 1.3% error from the measurement. Figure 23 shows the results of forced resp onse near the resonance in terms of transmissibility and frequency for b oth exp erimen ts and simulations. The transmissibility w as defined as the ratio of the steady-state amplitude of oscillation of the tip to that of the base ( X b ). The plots clearly show that the system has softening nonlinearit y , i.e., the resonan t p eak b ends tow ard low er frequency . The graph shows that the resonant frequency is low er for X b = 0 . 83mm than that for X b = 0 . 27mm. On the other hand, resonan t p eak amplitude is higher for X b = 0 . 83mm than that for X b = 0 . 27mm. These trends are w ell captured b y the mo del obtained by the prop osed approac h. 19 Magnet (fixed) Cantilever beam Magnet LDS to measure base displacement LDS to measure tip displacement Electrodynamic shaker Figure 21: T est equipment for forced vibration test 5 Conclusions In this pap er, a nov el data-driven nonlinear identification metho d for mechanical systems under nonlinear restoring forces has b een prop osed. The metho d is based up on the netw ork of piecewise linear springs represen ted by min and max functions. The contributions from the springs on the nonlinear restoring forces are obtained by solving a regression problem of the nonlinear forces with a library of min-max functions. The metho d has b een applied to a Duffing oscillator and a piecewise linear oscillator. It w as sho wn that the resulting mo dels can capture the free and forced resp onse under harmonic excitations. The metho d w as also applied to identify the magnetic restoring forces exerted on a cantilev ered plate. T o exploit the conserv ative nature of the magnetic forces, p oten tial constrain ts are applied to the min-max represen tation of the nonlinear force suc h that they represent conserv ativ e forces. The prop osed metho d successfully identified the nonlinear restoring forces exerted b y the p ermanen t magnet on the cantilev ered plate and the asso ciated p otential function and pro duced a mo del that can accurately predict free and steady-state forced resp onse of the system under harmonic excitation. References [1] S. L. Brunton and J. N. Kutz. Data-driven scienc e and engine ering . Cambridge Univ ersity Press, 2019. [2] Qianxiao Li, F elix Dietrich, Erik M. Bollt, and Ioannis G. Kevrekidis. Extended dynamic mo de decomp osition with dictionary learning: A data-driv en adaptive sp ectral decomp osition of the Ko opman op erator. Chaos: A n Inter disciplinary Journal of Nonline ar Scienc e , 27(10):103111, 10 2017. [3] Alessandro Alla and J. Nathan Kutz. Nonlinear mo del order reduction via dynamic mo de decomp osi- tion. SIAM Journal on Scientific Computing , 39(5):B778–B796, 2017. [4] Akira Saito and Masato T anak a. Data-driven mo del order reduction for structures with piecewise linear nonlinearit y using dynamic mo de decomp osition. Nonline ar Dynamics , 111(22):20597–20616, 2023. [5] Josh ua L. Pro ctor, Steven L. Brunton, and J. Nathan Kutz. Dynamic mo de decomp osition with con trol. SIAM Journal on Applie d Dynamic al Systems , 15(1):142–161, 2016. 20 T ime(s) 0 0.1 0.2 0.3 0.4 0.5 0.6 -2 0 2 Displacement(m) # 10 -3 T ip T ime(s) 0 0.1 0.2 0.3 0.4 0.5 0.6 -4 -2 0 2 4 Displacement(m) # 10 -4 Base (a) Exp eriment (phase shift: 178.20deg T ime(s) 0 0.1 0.2 0.3 0.4 0.5 0.6 -2 0 2 Displacement(m) # 10 -3 T ip T ime(s) 0 0.1 0.2 0.3 0.4 0.5 0.6 -4 -2 0 2 4 Displacement(m) # 10 -4 Base (b) Simulation (phase shift: 175.86deg) Figure 22: Comparison of the measured forced resp onse time histories and those computed with the prop osed approach at 15Hz. [6] Hideaki Namiki and Akira Saito. Short-time dynamic mo de decomp osition: a data-driven forced resp onse analysis of nonlinear systems under swept-sine excitation. Nonline ar Dynamics , 113(24):33109– 33126, 2025. [7] P . V. Kuptsov, A. V. Kuptso v a, and N. V. Stankevic h. Artificial neural netw ork as a universal mo del of nonlinear dynamical systems. R ussian Journal of Nonline ar Dynamics , 17(1):5–21, 2021. [8] P av el V. Kuptsov, Nataliy a V. Stankevic h, and Elmira R. Bagautdinov a. Disco vering dynamical features of Ho dgkin–Huxley-type mo del of physiological neuron using artificial neural netw ork. Chaos, Solitons & F r actals , 167:113027, 2023. [9] Anass Moukhliss, Elmahdi Ezzoubaidi, Nassima Ayoub, ab dellah Amouch, ihsane tikonab, Mohcine Cha jdi, Ab dellatif Rahmouni, and Rhali Benamar. Deep neural netw ork-based surrogate mo deling for nonlinear vibrations of functionally graded stepp ed b eams informed b y a discrete mo del. Disc over Me chanic al Engine ering , 5(1):8, 2026. [10] Lanxin Luo, Limin Sun, Yixian Li, and Y ong Xia. Structural nonlinear b oundary condition iden- tification using a hybrid physics data-driven approach. Nonline ar Dynamics , 113(9):9605–9623, 2025. [11] M. Raissi, P . Perdik aris, and G. E. Karniadakis. Physics-informed neural netw orks: A deep learning framew ork for solving forward and in verse problems in volving nonlinear partial differential equations. Journal of Computational Physics , 378:686–707, 2019. [12] Stev en L. Brunton, Joshua L. Pro ctor, and J. Nathan Kutz. Disco vering gov erning equations from data b y sparse identification of nonlinear dynamical systems. Pr o c e e dings of the National A c ademy of Scienc es , 113(15):3932–3937, 2016. 21 Figure 23: Comparison of the measured forced resp onse results and those computed with the prop osed approac h [13] F ederico J. Gonzalez and Luis P . Lara. Interpretable neural net work system iden tification metho d for t wo families of second-order systems based on characteristic curves. Nonline ar Dynamics , 113(24):33063–33086, 2025. [14] D. P . Garg. Piecewise linear approximation of nonlinear forces and moments resulting from crosswind gusts applied to magnetically levitated v ehicles. A dvanc es in A er osp ac e Structur es and Materials , pages 181–188, 1981. [15] Qing jie Cao, Marian Wiercigro c h, Ek aterina E P avlo vsk aia, J. Michael T Thompson, and Celso Greb ogi. Piecewise linear approach to an arc hetypal oscillator for smo oth and discontin uous dynamics. Philosophic al T r ansactions of the R oyal So ciety A: Mathematic al, Physic al and Engine ering Scienc es , 366(1865):635–652, August 2008. [16] P .-S. Xie and P .-J. Shih. A unique formulation of piecewise exact metho d to analyze a nonlinear spring system under harmonic excitation. Journal of Me chanics , 31(3):337–344, 2014. [17] Haining Li, Kefu Liu, and Jian Deng. Using a piecewise linear spring to approximate an essentially nonlinear spring: design and v alidation. T r ansactions of the Canadian So ciety for Me chanic al Engine ering , 49(2):192–206, 2024. [18] H.N Mhask ar and Charles A Micc helli. Approximation by sup erposition of sigmoidal and radial basis functions. A dvanc es in Applie d Mathematics , 13(3):350–373, 1992. [19] Kurt Hornik, Maxw ell Stinchcom be, and Halb ert White. Multilay er feedforw ard net works are univ ersal appro ximators. Neur al Networks , 2(5):359–366, 1989. [20] G. Cyb enk o. Approximation by sup erpositions of a sigmoidal function. Mathematics of Contr ol, Signals and Systems , 2(4):303–314, 1989. [21] Ry osuke Kanki and Akira Saito. Data-driven initial gap iden tification of piecewise-linear systems using sparse regression and universal approximation theorem. Journal of Computational and Nonline ar Dynamics , 19(6):061003, 05 2024. 22 [22] Diederik Kingma and Jimmy Ba. Adam: A metho d for sto c hastic optimization. International Confer enc e on L e arning R epr esentations , 2014. [23] Leonard Meiro vitch. F undamentals of Vibr ations . McGraw-Hill, 2001. [24] Meng-Hsuan Tien, Keng-Y en Lee, and Shih-Chun Huang. Analyzing the backbone curve of piecewise- linear non-smo oth systems using a generalized bilinear frequency approximation metho d. Me chanic al Systems and Signal Pr o c essing , 204:110765, 2023. [25] S. W. Shaw and P . J. Holmes. A p eriodically forced piecewise linear oscillator. Journal of Sound and Vibr ation , 90(1):129–155, 1983. [26] Huaxia Deng, Y u Du, Zhemin W ang, Jingchang Y e, Jin Zhang, Mengchao Ma, and Xiang Zhong. P oly-stable energy harvesting based on synergetic multistable vibration. Communic ations Physics , 2(21), 2019. [27] Haining Li, Kefu Liu, and Jian Deng. A magnetically enhanced piecewise-linear nonlinear energy sink: T ransient resp onses. Journal of Sound and Vibr ation , 626:119622, 2026. Statemen ts and Declarations Comp eting interests The authors hav e no relev an t financial or non-financial interest to disclose. Author contributions A. Saito: Conceptualization, Metho dology , F ormal analysis, Inv estigation (n u- merical), W riting – original draft, W riting – review & editing. H. F ujita: Resources, Inv estigation (exp erimen tal), Data curation, V alidation, W riting – review & editing. All authors appro ved the final man uscript. Data a v ailabilit y The dataset generated during and/or analyzed during the curren t study are a v ailable from the corresp onding author up on reasonable request. 23