## Abstract

In this paper, a new technique is presented for parametrically studying the steady-state dynamics of piecewise-linear nonsmooth oscillators. This new method can be used as an efficient computational tool for analyzing the nonlinear behavior of dynamic systems with piecewise-linear nonlinearity. The new technique modifies and generalizes the bilinear amplitude approximation method, which was created for analyzing proportionally damped structural systems, to more general systems governed by state-space models; thus, the applicability of the method is expanded to many engineering disciplines. The new method utilizes the analytical solutions of the linear subsystems of the nonsmooth oscillators and uses a numerical optimization tool to construct the nonlinear periodic response of the oscillators. The method is validated both numerically and experimentally in this work. The proposed computational framework is demonstrated on a mechanical oscillator with contacting elements and an analog circuit with nonlinear resistance to show its broad applicability.

## 1 Introduction

Piecewise-linear (PWL) systems are often utilized to model the dynamic behavior of many physical systems with nonsmooth properties. The analysis of the dynamics of these systems is of importance in many fields including health monitoring of engineering structures [1–3], design of electrical circuits [4,5], and gene regulatory networks modeling in biological systems [6,7]. The analysis of the nonlinear behavior of these nonsmooth systems is usually computationally challenging, since the techniques that have been developed for analyzing nonlinear dynamical systems require much more computational cost than traditional linear methods [8,9]. Thus, it is crucial to develop efficient computational tools for predicting the dynamics of nonsmooth systems so that efficient parametric studies can be performed in a reasonable time frame.

The nonsmooth properties of engineering systems are usually triggered by discrete events in their dynamic responses. For instance, the intermittent contact observed in engineering structures due to cracks or delamination [10–13], or interface among components [14] often results in nonsmooth features. Also, the nonsmooth dynamics and mechanics are of great importance for modeling jointed structures [15]. A great amount of research has been carried out on the nonlinear modeling of mechanical joints over the last decade [16–18]. Furthermore, the nonsmooth property also occurs in many electrical systems when current-limit elements such as diodes are utilized to realize particular functions [4]. Since the nonlinearity induced by the discrete event usually dominates the dynamics of these nonsmooth systems, they are often modeled using PWL systems in which the dynamics are assumed to behave linearly in each of the systems' subregimes. Therefore, an effective computational tool that can capture the dynamic features of PWL systems is required for analyzing these nonsmooth systems.

Although PWL systems usually consist of multiple linear subsystems, traditional liner analysis techniques such as modal analysis [19] are not able to capture the event-triggered nonlinearity in these systems since they usually exhibit strong nonlinear features [5,20]. Therefore, numerical integration (NI) is usually required to compute the time response of these nonlinear systems. Common NI techniques include Runge–Kutta methods [21] and Newmark methods [22]. Although these methods are commonly employed in nonlinear system analyses, NI usually requires a fairly small time-step to predict the nonsmooth evolution of nonlinear variables in nonsmooth systems and hence reduces computational efficiency significantly. Recently, a new class of transient dynamic analysis tools has been developed to speed up time marching in the simulation process [23,24]. These techniques are referred to as the hybrid symbolic numeric computational (HSNC) method and have been shown to be more efficient than NI for simulating the dynamic response of PWL systems. The HSNC method employs the analytical solution of the system in each of its linear regimes and stitches these PWL responses together using a numerical incremental search process. Since the method only requires numerical computation at transition points where the system switches its state, the computational cost can be reduced significantly compared to traditional NI methods. The HSNC method is particularly useful for analyzing transient and chaotic responses of complex PWL systems. The method has been recently adapted to analyze the nonlinear rattling behavior of gear systems [25,26].

Although efficient transient dynamic analyses for PWL systems can be achieved using HSNC, a method that can directly capture the steady-state dynamics of these systems is preferred when parametric studies are needed. Thus, another branch of techniques based on the harmonic balance (HB) method has been developed to capture the steady-state vibrational response of PWL nonsmooth systems [27–31]. These methods have been applied to a broad array of nonlinear systems and are not limited to PWL applications. These HB-based methods approximate the periodic response of nonlinear systems using a truncated Fourier series with the unknown coefficients being solved using numerical methods. Although these methods are more efficient than transient dynamic analysis tools, numerous harmonics are usually required to capture the nonsmoothness in PWL systems in particular. Thus, these methods still require considerable computational cost and often incur convergence issues when dealing with PWL systems.

Recently, a new efficient method has been created to approximate the steady-state dynamic response of PWL systems. This method is referred to as the bilinear amplitude approximation (BAA) method [32] and is applicable for capturing the vibrational response of a subset of PWL nonsmooth systems whose dynamic behavior is dominated by two distinct linear states. The BAA method is based on a similar idea as that of HSNC; namely, the nonlinear dynamic response of PWL systems can be obtained by stitching together the analytical responses of the system in each of its linear regimes. The method derives a set of constraint equations to enforce the continuity condition of the PWL responses within a vibrational cycle. A numerical optimization solver is then used to solve for the unknown parameters in the analytical responses when combining each individual linear response to obtain the nonlinear vibrational response. It has been shown that the BAA method is more efficient than the transient dynamic analysis tools and HB-based techniques [32]. Furthermore, the BAA method does not have similar convergence issues as that of HB-based methods, since this method only requires a few linear modes to obtain an accurate approximated response. The BAA method has been successfully applied to investigate the nonlinear dynamics of complex structures with cracks [33]. Also, the method has been extended to develop an efficient technique that is able to capture the vibrational response of systems with Coulomb friction [34]. Note that the idea of utilizing a semi-analytic approach for computing the dynamic response of PWL systems was also used in some previous studies [35,36]. These studies are focused on investigating the nonlinear behaviors of single degree-of-freedom (DOF) oscillators in contrast to the BAA method, which was developed for enabling the analysis of high-dimensional structural systems that contain stationary dynamics.

Although the BAA method is an efficient and effective tool, it is developed based on the standard form of equations of motion of structural systems. It also requires the system to be proportionally damped since its formulation to obtain the analytic solution of the modal response is built off of second-order differential equations. Hence, the goal of this paper is to extend the original BAA method to more general PWL systems so that physical systems in a variety of disciplines (e.g., electrical engineering) modeled using state-space models can also be efficiently investigated. In the extended method, the complex analytical expression of PWL responses of bilinear state-space models is derived, and a numerical tool is used to find the parameters in the analytical expressions by enforcing appropriate compatibility conditions. Nonlinear vibrational responses of the system can then be obtained by coupling the analytical solutions of the system in each linear regime. The extended method is as efficient as the original BAA method since they both utilize efficient analytical approaches in determining the PWL responses. However, the extended method is also able to analyze state-space models whose modal properties are quantified using complex numbers. This capability significantly extends the applicability of the method. Also, the proposed method is validated by both numerical simulation and experimental measurement in this work.

The remainder of the paper is organized as follows. First, the proposed method is introduced. Next, the results of applying the proposed method for a spring-mass mechanical oscillator and a nonlinear electrical circuit are presented. Finally, conclusions are discussed.

## 2 Methodology

where $unl$ represents the nonlinear state variable that determines the system's state, $A1$ and $A2$ represent the state matrices in the first and second linear states, and $B1$ and $B2$ represent the input matrices in the first and second linear states. These matrices are all constant matrices since they represent the system matrices of the linear subsystems. Note that *a* in Eq. (2) represents a threshold that distinguishes the first and second linear states.

**A**'s can be expressed as

**A**'s. Note that if $\lambda j$ and $xj$ are a set of complex eigenvalue and eigenvector, then the complex conjugates $\lambda \xafj$ and $x\xafj$ are also a set of eigenvalue and eigenvector. Next, the adjoint eigenvalue problem of the transposed matrices $AT$ can be formulated as

**A**share the same eigenvalues, and the right eigenvectors $xj$'s and the left eigenvectors $yr$'s satisfy the following biorthogonality:

Note that $nj(t)=yjTBq(t)$.

where $\xi j,1(t)$ ($j=1,\u2026,n$) represents the modal response of the PWL system in the first state; $\xi r,2(t)$ ($r=1,\u2026,n$) represents the modal response of the system in the second state; $cj,1,\u2009cr,2\u2208\u2102$ are scalar complex amplitude coefficients of the transient responses; $\omega \u2208IR$ represents the frequency of excitation; and $\psi \u2208IR$ represents the phase difference between the excitation and the steady-state responses. Note that *ψ* is included as an additional phase shift so that the time origin can be reset when coupling these PWL responses. The coupling process will be explained shortly. Also note that the phrases “transient response” and “steady-state response” represent the responses of the system in each of its linear states and should not be confused with the total nonlinear response of the nonsmooth system. Furthermore, if *ξ _{j}* associated with the eigenvector $xj$ is a complex coordinate, its complex conjugate $\xi j\xaf$ is also a modal coordinate and the eigenvector associated with $\xi j\xaf$ is $xj\xaf$, the complex conjugate of $xj$. Thus, complex coefficients

*c*and $c\xafj$ exist in conjugate pairs. The physical coordinates of the system in its linear regimes can be obtained by applying the transformation $u1(t)=X1\xi 1(t)$ and $u2(t)=X2\xi 2(t)$, where the subscripts 1 and 2 represent the two states in Eq. (8). The key idea of the method is to couple $u1(t)$ and $u2(t)$ to construct the nonlinear vibrational response of the PWL system by applying appropriate compatibility conditions. Note that the modal expression Eq. (11) used in this work is derived for complete eigenproblems. For systems whose state matrices are defective, the generalized eigenvectors are required to complete the eigenbasis and Eq. (11) needs to be modified for those cases [37].

_{j}*T*

_{1}represents the period of time the system stays in the first linear state;

*T*

_{2}represents the period of time the system stays in the second linear state; and

*T*represents the full period of the nonlinear vibrational cycle. In order to construct this vibrational cycle, a set of compatibility conditions is derived and applied to obtain the unknowns $cj,1,\u2009cr,2$, and $\psi $ in Eq. (11)

*u*at the moments of transition. Note that the period of vibration

_{nl}*T*can be determined by the user. If the system responds at the excitation frequency, then the period of vibration can be calculated by $T=2\pi \omega $. If sub- or super-harmonic motions are of interest, then, the period is assumed to be a multiple or a factor of $T=2\pi \omega $. Furthermore,

*T*

_{1}in Eq. (12) is an additional unknown that needs to be obtained to construct the full nonlinear vibrational cycle and

*T*

_{2}can be calculated by $T\u2212T1$ once

*T*

_{1}is obtained. Next, Eq. (8) is employed in Eq. (12) to obtain the modal representation of the compatibility conditions

where $Xnl,1$ is the portion of the mode shapes corresponding to the nonlinear variable $unl,1$. In other words, $Xnl,1$ represents the *k*th row of $X1$, where the *k*th variable in $u(t)$ is the nonlinear state variable *u*_{nl}(*t*). The unknowns ($cj,1,\u2009cr,2,\u2009\psi $, *T*_{1}) in Eqs. (11) and (13) can be solved using numerical optimization tools. The optimization solver “lsqnonlin” in matlab [38] was used in this work to find these parameters by minimizing the residual in Eq. (13). The solver uses the trust-region-reflective algorithm to iteratively search for a local minimum of the objective function until a tolerance is achieved. Note that if the system is driven by $q0\u2009cos\u2009\omega t$, only the real part of Eq. (13) needs to be solved; if the system is excited by $q0\u2009sin\u2009\omega t$, only the imaginary part of Eq. (13) needs to be solved. In practice, one can arbitrarily choose either the real or imaginary part to solve since the phase delay in the excitation does not affect the steady-state vibration properties of the system. A steady-state vibrational cycle of the PWL system can then be constructed by coupling the solved linear responses of each linear state.

Since the proposed method is able to obtain the steady-state response of PWL systems, it is particularly useful for conducting parametric studies. A parametric sweep can be efficiently performed by starting from an initial parameter value. At this initial parameter value, multiple random initial guesses for ($cj,1,\u2009cr,2,\u2009\psi $, *T*_{1}) are provided to the nonlinear optimization solver and the solution with the minimum residual of Eq. (13) is chosen as the initial values for the second parameter point. In the remaining sweep process, the solution from the previous parameter value is used as initial values for computing the variables at the next parameter point. This sweep strategy ensures convergence while speeding up the analysis. Note that although the method is proposed for PWL systems with two distinct linear subsystems in this work, Eq. (13) can be expanded to include additional compatibility conditions for systems with more linear states (i.e., more general PWL systems that are not simply bilinear). Also note that the stability of the obtained solutions can be determined by the modified Floquet theory proposed by Leine [39].

## 3 Results

The proposed method is demonstrated on both mechanical and electrical systems to verify its broad capability of analyzing PWL nonsmooth systems in general fields. Results of applying the methodology to a mechanical oscillator with contacting elements and a nonlinear analog circuit with a diode are discussed.

### 3.1 Mechanical Oscillator With Contacting Masses.

*g*=

*0 m. Note that the damping coefficients are chosen such that the system is not proportionally damped and must be modeled using a state-space representation; thus, the system studied in this work represents a fairly general mechanical oscillator that the original BAA method [32] would not be able to handle. In this case study, a harmonic excitation $F(t)=f0\u2009cos(\omega t)$ with amplitude $f0=0.01$ N is applied on mass*

*m*

_{3}. In order to model the intermittent contact phenomenon between the masses

*m*

_{1}and

*m*

_{2}, two linear subsystems are used to model the dynamics of the nonsmooth oscillator in each of its linear states. The first subsystem shown in Fig. 2(b) represents the linear state of the oscillator when the gap is open; the second subsystem shown in Fig. 2(c) represents the linear state of the oscillator when the gap is closed. In the second subsystem, a set of contact stiffness $k*$ and damping $c*$ is added to the model for minimizing the penetration between the contacting masses. In this case study, the values $k*=1000$ N/m and $c*=50$ kg/s are used to model the contact stiffness and damping, which are at least two orders of magnitude larger than the structural stiffness and damping. The equations of motion of the nonsmooth oscillator can be expressed as

*m*

_{1},

*m*

_{2}, and

*m*

_{3}, respectively. The matrices in Eq. (14) can be expressed as

where **I** represents the identity matrix.

*m*

_{1}and

*m*

_{2}, the nonlinear variable for this mechanical oscillator is $unl=x1\u2212x2$ and the threshold value is set to

*a*=

*0. In order to apply the proposed algorithm, the original state vector*

**u**is converted to an appropriate coordinate system $u\xaf=[unl,x2,x3,x\u02d91,x\u02d92,x\u02d93]T$ by applying the transformation $u=\alpha u\xaf$, where

The system matrices associated with the converted state vector can also be obtained by applying the transformation $A\xaf1=\alpha TA1\alpha ,\u2009A\xaf2=\alpha TA2\alpha ,\u2009B\xaf1=\alpha TB1\alpha $, and $B\xaf2=\alpha TB2\alpha $. Note that *u _{nl}* is the only variable that determines the state of the system in the converted coordinate system. The proposed method is then applied to the transformed system, and the computed forced response is compared with the one obtained using NI in Fig. 3. For performing NI, the explicit Runge–Kutta method [21] and the event function in matlab [38] are employed in the computation. Figure 3 shows that the steady-state responses obtained using the proposed method and NI have an excellent agreement over the plotted frequency range in which multiple resonances are accurately captured by the extended BAA method. Furthermore, the proposed method only requires 0.039 s to obtain the steady-state vibrational response for a specific frequency. By contrast, the average CPU time required by NI is 2.557 s. The computations are conducted using a Dell XPS 15 laptop (Round Rock, TX) (2.60 GHz), and the extended BAA method only requires 1.53% of the CPU time of NI. The proposed method is expected to have increasing computational savings with respect to NI in the same way as the original BAA method as the system investigated becomes more and more complex [33,40].

### 3.2 Nonlinear Analog Circuit With a Nonsmooth Resistance.

*C*), an inductor (

*L*), several resistors (

*R*,

*R*

_{1},

*R*

_{2}), and a diode. Note that the resistor in the circuit is PWL since the diode restricts the current flow to one direction. If the circuit system is excited by an AC voltage source with frequency $\omega $ and amplitude $f0$, the system can be modeled as

*v*(

*t*) represents the voltage across the capacitor

*C*, and the state variable

*i*(

*t*) is the current across the inductor

*L*and the resistor

*R*. The system matrices in Eq. (18) can be expressed as

The steady-state response of the circuit is analyzed using the proposed method. In this study, the parameter values are chosen as $R=50\u2009\Omega $, $L$* *=* *10 mH, $R1=2500\u2009\Omega ,\u2009R2=500\u2009\Omega $, and $f0=9$ V. The phase portraits of the system for a few capacitor values ($C=0.01\u2009\mu $F, $1\u2009\mu $F, and $100\u2009\mu $F) and driving frequencies ($\omega =10\u2009$Hz, $100\u2009$Hz, $1\u2009$kHz, and $10\u2009$kHz) computed using the proposed technique and NI are compared in Fig. 5. In each plot, $(\u2212)$ and $(\u25ef)$ represent the solutions computed using the proposed technique and NI, respectively. Note that the plots in the same row represent phase portraits with the same driving frequency, and the plots in the same column represent phase portraits with the same capacitor values. Figure 5 shows an excellent agreement between the solutions computed using the proposed method and NI solutions. Furthermore, the system exhibits stronger nonlinearity as the capacitor value and/or the driving frequency is decreased since the trajectories in these phase portraits are more twisted.

Next, the extended BAA method is used to conduct a frequency sweep to validate its capability of capturing the dynamic response in a broad parameter range. The peak-to-peak voltage of the system for a couple of capacitor values ($C=1\u2009\mu $F and $C=100\u2009\mu $F) is shown in Figs. 6(a) and 6(b), respectively. The excitation frequency range in these case studies is chosen as $[10\u2009Hz,104.5\u2009Hz]$. The frequency sweep curves computed using the proposed method (–) are compared with the ones obtained by NI ($\xb0$) and experimental measurement (+). It is shown that the responses obtained using the proposed method, NI, and experimental measurement exhibits excellent agreement in the plotted frequency. The average CPU time required by the new method to compute the vibrational response for a specific frequency is 0.08 s. By contrast, NI requires 3.01 s to simulate to a steady-state response at each frequency. The new method only requires 2.7% of the CPU time of NI. Thus, the proposed method can be used as an efficient computational framework for designing circuit systems when current-limiting elements are required. Furthermore, this work provides the first experimental validation of the proposed hybrid analytic-numeric method and shows its excellent agreement for analyzing an electrical system. The original BAA method was only validated on a structural system [41].

Next, a multivariable parametric study is conducted using the proposed method to analyze the dynamic properties of the circuit. The driving frequency *ω* and the capacitor *C* are chosen as the control parameters in this study. Note that this type of parametric study can be efficiently performed using the proposed method and is generally computationally challenging using traditional NI. The peak-to-peak amplitude of voltage and current for the frequency range $[10\u2009Hz,104.5\u2009Hz]$ and capacitor range $[10\u22122\u2009\mu F,102\u2009\mu F]$ are computed and plotted in Figs. 7(a) and 7(b), respectively. Figure 7 shows that the circuit exhibits a low-pass behavior since the amplitude of voltage drops quickly after a cutoff frequency. Also, the bandwidth shrinks as the capacitor value increases. Furthermore, it is observed that the amplitude of current rises around the cutoff frequency. Although this nonlinear circuit exhibits similar low-pass behavior as a simple linear resistor-inductor-capacitor circuit, the system is different due to how the diode element affects the electrical load and the nonlinearity must be accounted for when analyzing it. The investigation of these properties is critical since the system exhibits a stronger nonlinear behavior in the lower frequency range as discussed above.

## 4 Conclusions

In this paper, an efficient methodology for analyzing the steady-state periodic response of piecewise-linear (PWL) nonsmooth systems in general fields is introduced. The technique extends the bilinear amplitude approximation method, which was developed previously for capturing the dynamic response of proportionally damped structural systems, to more general systems so that PWL oscillators modeled using state-space representations can also be efficiently investigated. The new method employs the eigen-decomposablility of a state-space model in distinct linear regimes to obtain the closed-form solution for each linear subsystem. The steady-state vibrational response of the PWL oscillator is then constructed by coupling the closed-form solutions with appropriate compatibility conditions. An efficient computation of the nonlinear dynamics can be achieved using the proposed hybrid analytical-numeric approach.

The method is successfully demonstrated on a mechanical oscillator with contacting elements and an electrical circuit with a diode to show its broad applicability. Also, the method is validated both numerically and experimentally. The proposed method shows the capability of capturing the nonlinear vibrational response of PWL nonsmooth oscillators for a wide parameter range and is an efficient computational framework for investigating the effects of multiple parameters on the system dynamics.

## Acknowledgment

This paper is based on work partially supported by the National Science Foundation (United States) under Grant No. 1902408, program manager Dr. Harry Dankowicz, and the Ministry of Science and Technology (Taiwan, R.O.C) under Grant No. MOST 109-2222-E-007-006-MY3. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation and the Ministry of Science and Technology.

## Funding Data

Ministry of Science and Technology, Taiwan (Grant No. MOST 109-2222-E-007-006-MY3 Funder ID: 10.13039/501100004663).

National Science Foundation (Grant No. 1902408; Funder ID: 10.13039/100000001).