Model reduction of a periodically forced slow–fast continuous piecewise linear system

In this paper, singular perturbation theory is exploited to obtain a reduced-order model of a slow–fast piecewise linear 2-DOF oscillator subjected to harmonic excitation. The nonsmooth nonlinearity of piecewise linear nature is studied in the case of bilinear damping as well as with bilinear stiffness characteristics. We propose a continuous matching of the locally invariant slow manifolds obtained in each subregion of the state space, which yields a reduced-order model of the same nature as the full dynamics. The frequency-response curves obtained from the full system and the reduced-order models suggest that the proposed reduction method can capture nonlinear behaviors such as super- and subharmonic resonances.


Introduction
Piecewise linear (PWL) systems constitute an important class of nonsmooth dynamical systems and are commonly used to model complex physical phenomena.Typical applications are found in a wide range of fields, such as control systems [1], neuroscience [2] and mechanical engineering.For instance, the presence of effects such as dry friction [3], intermittent contact or bilinear stiffness characteristics [4] yields PWL mechanical models consisting of multiple linear subsystems.A typical example is the quarter car model, which is widely used to study the vertical dynamics of the suspension of a single wheel under road excitation [5].A variation of this simplified low-dimensional example including bilinear damping is investigated in the current work.
The dynamics of PWL systems can exhibit very interesting dynamical effects, which are impossible to observe in smooth systems.A prominent example of this rich dynamic behavior was reported in [6], where the continuous matching of two stable linear subsystems can result in an unstable dynamics.The regularization of PWL nonsmoothness using smooth functions is often used to circumvent the typical problems emanating from the nonsmooth nature of the models.However, these smooth representations may not be well suited for the study of general PWL systems and are not always able to capture their typical nonlinear behaviors, as pointed out in [7] and [8].Therefore, PWL systems cannot directly be studied with classical methods from the theory of smooth nonlinear systems.Instead, these methods need to be extended to take the nonsmoothness into account.
In the past three decades, a number of researchers have sought to determine the qualitative behavior of PWL mechanical systems.Shaw et al. [9] provide a detailed analysis of a periodically forced singledegree-of-freedom PWL oscillator, where the limiting case of an infinite stiffness slope represents a rigid impact oscillator.Wang et al. [10] proposed an analytical method to determine the amplitude frequency and phase frequency characteristics for single-degreeof-freedom forced oscillators with nonlinear restoring forces, which can be smooth or of PWL nature.More recently, Stefani et al. [11] carried out a numerical bifurcation analysis of a symmetrically constrained single-degree-of-freedom vibro-impact oscillator, taking effects into account such as impact, clearance and unilaterality of the constraints.Path-following techniques were used by Chavez et al. [12] to study the dynamics of a PWL capsule system.A semi-analytical method to study a single-degree-of-freedom impact oscillator with PWL stiffness and damping under harmonic excitation was proposed by Hernández Rocha et al. [13].This approach is based on a mapping technique using the closed-form solutions available in each region of the PWL system and allows to describe any periodic motion and study its stability and bifurcations.For autonomous PWL multi-degree-of-freedom systems, which model cracked rotating shafts and cracked beams, Zuo et al. [14] were the first to investigate the nonlinear modal behavior in the gyroscopic and nongyroscopic case.Chen et al. [15] introduced a method to construct the nonlinear normal modes (NNMs) for a more general class of undamped autonomous PWL systems using the concepts of Poincaré maps and invariant manifolds.An extension of this study was proposed by Jiang et al. [16] using Galerkin projections to account for large vibration amplitudes.Chati et al. [17] used perturbation methods to obtain the nonlinear normal modes of an autonomous conservative two-degree-offreedom PWL oscillator.The NNMs are based on the bilinear frequency formula and are accurate only when the difference between the linear subsystems is small.Similarly, Butcher [4] investigated the effects of a clearance on the normal mode frequencies of undamped ndimensional systems with bilinear stiffness.For systems in which the linear part of the model dominates in the dynamics, Butcher et al. [18] proposed an approach for model order reduction in the presence of PWL stiffness.An analytical approach was proposed by Bellizi et al. [19] to predict the periodic and quasiperiodic response regimes of a two-degree-of-freedom system with piecewise nonlinear damping and time delay.
From the above review, it becomes clear that most of the existing techniques to study PWL systems are either restricted to low-dimensional systems or do not include effects such as damping and external forcing.However, the analysis of real-world applications requires highdimensional multibody mechanical models, which in turn present an additional challenge to the study of the PWL dynamics.Hence, adequate model order reduction tools are needed to investigate the dynamics efficiently, aiming at reducing the computational efforts while also providing an insight into the main features of the global dynamics.
For smooth mechanical systems, perturbation methods have been used to derive reduced-order models by assuming a global slow-fast decomposition of the coordinates, such as in [20,21] and [22].A more recent order reduction methodology is provided by the local theory of spectral submanifolds [23], which presents a unified mathematical approach to NNMs of general non-autonomous dissipative systems.This technique relies on a slow-fast decomposition of the state variables that results in a lower-dimensional manifold, on which the slow dynamics of the full system can be investigated.This allows to obtain a global reduced model capturing the key characteristics of the steadystate system behavior.However, this approach relies on smoothness properties of the system and its starting point consists in the study of the underlying linear dynamics.In the case of PWL systems, a proper linearization does not exist, in the sense that an equilibrium located on a switching manifold does not admit a neighborhood around which the dynamics can be linearized.Nonetheless, a slow-fast decomposition analogous to smooth systems may provide a promising first step toward a reduction to lower-dimensional invariant manifolds (or their nonsmooth generalizations).Therefore, the investigation of slow-fast PWL systems using perturbative approximations could pave the way for the development of novel reduction methods accounting for PWL nonlinearities.In the smooth case, slow-fast systems have been studied extensively in the context of bistable oscillators and mixed-mode bursting and relaxation oscillations, which occur in mechanical systems among other applications, such as in the work by Kovacic [24] and Rakaric [25].In this study, the focus lies rather on exploiting the slow-fast decomposition of the state variables for the purpose of model order reduction of PWL mechanical systems.
In this paper, the original approach of the authors proposed first in [26] for the use of singular perturbation theory on slow-fast continuous PWL systems is extended from the autonomous configuration in R 3 to non-autonomous slow-fast oscillators with two degrees of freedom.In the three-dimensional state space as described in [26], the reduction method consisted in coupling only one fast variable to two slow variables.However, the addition of a second fast variable introduces an additional complexity since the interplay between the two fast variables at the switching manifold affects the proposed continuous matching procedure.Furthermore, in the autonomous case studied in [26], the trajectories of the full and reduced systems converge toward the equilibrium along timeindependent half-manifolds.Extending the analysis to the non-autonomous setting results, however, in explicitly time-dependent linear half-manifolds, which are continuously matched to approximate a periodic solution.The novelty of our proposed continuous matching technique is that the resulting reduced model inherits the continuity property of the original system.This feature is important for two reasons.First, a general PWL reduced model can exhibit the complex dynamic behaviors mentioned above.Therefore, determining its stability properties is more complex than that of the underlying original continuous system.Furthermore, the continuity property ensures that trajectories cross the switching manifold from one subregion to the other, thus excluding sliding along the switching hyperplane as may occur in Filippov systems [27].A discontinuous reduced model may, however, exhibit sliding solutions, which are not present in the original continuous system.Moreover, the presence of sliding behavior in the reduced model can cause numerical problems (it necessitates the numerical solution of a differential inclusion) and increases the computational effort, which is not desirable from a reduction point of view.
Motivated by the quarter car model with bilinear damping characteristics, two examples of slow-fast oscillators subjected to a harmonic excitation are used to illustrate the reduction.
The remainder of this paper is organized as follows.In Sect.2, we carry out the direct singular perturbation reduction on the quarter car model with bilinear damping characteristic to illustrate the emerging challenges of matching the slow dynamics of each subsystem at the switching manifold, motivating the need for a continuous matching procedure of the slow half-manifolds.In Sect.3, such a continuous matching approach is proposed in a general setting aiming to obtain a continuous reduced-order model.In Sects. 4 and 5, the reduction using the continuous matching approach is applied to two continuous PWL systems with bilinear damping and stiffness characteristics, respectively.Reduced-order models are obtained and the resulting frequency-response curves are compared to the dynamics of the full systems.Finally, conclusions are drawn in Sect.6.

Motivating example for singular perturbation theory on continuous PWL systems
In order to enhance the driving comfort and vehicle safety, various suspension models have been developed with the aim of preventing up-swings of the vehicle and providing a good road grip.A widely used model to investigate the dynamics of suspension systems is the quarter car model with two degrees of freedom.Variations of this system with nonlinear elements, such as bilinear dampers and unilateral springs, have been studied extensively.For instance, Silveira et al. [5] analyzed a 2-DOF model with bilinear damping under harmonic excitation using the harmonic balance method.
A numerical stability analysis of a similar 2-DOF system has been conducted by Szabó [31].Moreover, the incremental harmonic balance method has been utilized by Wang et al. [32] to study a suspension model with piecewise linear stiffness and damping.A 2-DOF quarter car model with bilinear damping is illustrated in Fig. 1a and consists of two masses, describing the car body and the wheel, respectively, while the tire of the wheel is simplified by a linear spring.In a typical car, the natural frequency of the tire is much larger than the natural frequency of the car body, which leads to the introduction of a small parameter responsible for a time-scale separation in the behavior of the system variables [33].Motivated by this effect and in order to simplify the embedding of the model in a singularly perturbed framework, we choose to model the system with the set of scaled parameters as shown in Fig. 1a.The car body has the mass m and its displacement is described by the absolute coordinate x, whereas the vertical position of the wheel of small mass εm is denoted by z.Both masses are connected by a suspension of stiffness k and a bilinear damping element with switching positive damping coefficients c + and c − .The characteristic of the damping force is illustrated in Fig. 1b, to account for a varying damping depending on the relative velocity.The road profile is modeled as a small harmonic excitation εr (t) of the form r (t) = R sin( t).The equations of motion of the system read as where the damping coefficient is defined as By introducing a set of new coordinates the system equations can be rewritten in the singularly perturbed first-order form with the slow and fast variables defined by x = x 1 x 2 T and y = y 1 y 2 T , respectively, the switching function defined as h(x, y) = y 2 and the continuous linear functions f ± = f ± (x, y, t; ε) and g ± = g ± (x, y, t; ε) given by The switching hyperplane is denoted by with state vector x = x 1 x 2 y 1 y 2 T .The motivating goal for the current work is to obtain a singular perturbation reduction for continuous PWL (hereafter, CPWL) systems similar to smooth systems, also taking external harmonic forcing into account.Before we embark on the application of the reduction method to the quarter car model with bilinear damping, we recall the main idea of singular perturbation theory.
For an n-dimensional smooth system having s slow variables and a small perturbation parameter ε, which is responsible for a time-scale separation, classical geometric singular perturbation theory can be used to obtain a reduced-order model [28].The limiting case ε = 0 gives an s-dimensional critical manifold M c .According to Fenichel's theorem [28], if M c is normally hyperbolic, then there exists an s-dimensional slow invariant manifold M s , on which the dynamics is a regular perturbation of the dynamics on M c .This theorem can be applied to slow-fast CPWL systems only on the subsets of the state space that do not include the switching manifold.This yields two linear locally invariant slow half-manifolds M ± s for the linear subsystems.Furthermore, a forward invariant neighborhood enveloping the linear critical manifold, which is continuous at the switching hyperplane, has been shown to exist under suitable conditions [29].For this, the linear critical manifolds M ± c , meeting continuously at , need to have an attracting fast dynamics (for ε > 0) within their respective linear subsystems.The next section presents a proof and discussion of the continuity property of the critical manifolds within the scope of CPWL systems.
Setting the small parameter ε = 0 in (5) yields the critical dynamics ẋ = f ± (x, y, t; 0) (8a) where the ± switch is governed by h(x, y) ≷ 0. The algebraic constraints (8b) admit one isolated solution y = h ± c (x, t) each, which describe the behavior of the fast variables as linear functions of the slow variables x along the critical manifolds of each linear subsystem.In the following, we define the critical "half-manifolds" as For the quarter car model governed by the slow and fast dynamics given in (6), solving the algebraic equations for y on each subregion yields the expressions where the external forcing is only present in the approximation of the first fast variable y 1 .The critical manifold is therefore defined as It is important to note at this stage that the original switching condition y 2 ≷ 0 simplifies to x 2 ≷ 0 in view of the second line in (9).The continuity of M c (t) at the switching manifold easily follows from h + c = h − c for y 2 = x 2 = 0 and is treated in the next section for general CPWL systems.The reduced dynamics on the critical manifold reads as The critical manifold gives a zero-order approximation of the slow manifold, where the fast variables are coupled to the slow variables by the expressions h ± c (x, t).The stability of the critical manifold can be obtained by the Jacobians of the fast dynamics which are given by Therefore, evaluating the trace and the determinant of B ± yields the stability of the fast dynamics for positive parameter values k, m, c ± > 0. From this point, the existence of two locally invariant slow manifolds for the + and − subsystems is established.As shown in [26], the zeroth-order approximation may result in conservative dynamics for specific parameter values, which does not reflect the dissipative nature of the original system and therefore cannot be used as a qualitative approximation of the full dynamics.In the present case of the CPWL quarter car model, setting ε = 0 and substituting the critical dynamics (11) yields an autonomous oscillator since all timedependent forcing terms only appear in h ± c,1 , which is canceled out by ε = 0. Thus, the reduced dynamics along the critical manifold cannot be used to approximate the full system and first-order terms have to be included in order to approximate the slow locally invariant manifolds from each subsystem and obtain a suitable reduced-order model.
The two-dimensional, locally invariant slow halfmanifolds are defined in the corresponding regions of the state space as Since the state space is decomposed into two linear parts, the invariance property of M ± s (t) must be understood in a local way.Inserting y = h ± s and ẏ = Within the linear regions h(x, y) ≥ 0 and h(x, y) < 0, the asymptotic expansions given by are substituted into their corresponding invariance Eq. (13).By equating the coefficients of powers of ε, it follows that h ± c (x, t) = h ± 0 (x, t), which means that in each linear region the critical manifold is the zeroorder approximation of the slow manifold.Moreover, the dynamics on the slow manifold is a regular perturbation of the critical dynamics.Furthermore, the firstorder terms h ± 1 (x, t) are obtained as where [•] 0 ± denotes evaluation at y = h ± 0 (x, t) and ε = 0.
By applying this procedure on both subsystems of the bilinear quarter car model, one can obtain an approximation up to order O(ε 2 ) of the two linear locally invariant slow manifolds.From Eq. ( 15), the first-order terms of the slow manifolds read as with Using Eq. ( 14), the approximation of the linear locally invariant slow half-manifolds up to O(ε 2 ) is obtained as where The switching condition in this case is not trivial anymore, since the choice of the damping coefficients c ± depends on h s (x, t) itself.As shown in [26], this dependency leads to regions in the state space, where either none or both of the switching conditions are valid.This may lead to numerical problems while evaluating the reduced dynamics on the slow manifolds at the crossing region.To circumvent this problem of switching between h ± s , one could take x 2 = 0 as a switching hyperplane, which corresponds to the switching condition on the critical manifold (see Eq. ( 10)).At this modified switching hyperplane, the linear slow manifolds, obtained as two half-planes , are askew and meet only at the equilibrium point.This leads to a reduced system containing a jump at x 2 = 0 given as follows: This discontinuous reduced-order model (hereafter, ROM) is a Filippov system characterized by a vector field which is bounded but discontinuous on the switching manifold , defined here as The state space is therefore split into two regions in each of which the system is determined by a smooth vector field.
At the switching manifold, the Filippov system can exhibit sliding motions, as opposed to a CPWL system, in which only direct crossing behavior can occur.For a general PWL Filippov system, the switching manifold is subdivided into different regions, depending on the interaction between the neighboring vector fields.For a point x * ∈ on the switching manifold at a time instant t * , the projections of the vector fields on the normal vector to pointing into the region x 2 > 0 are given by e T 2 f + disc (x * , t * ; ε) and e T 2 f − disc (x * , t * ; ε), respectively.The direct crossing set is the set of all points of , where both vector fields have nontrivial projections on in the same direction and is hence defined as where f ± disc = f ± disc (x * , t * ; ε).The sliding set, where both vector fields are pointing toward or away from , reads as On the sliding set, the solution concept needs to be extended to Filippov's solution concept [34], i.e., the discontinuous differential equation is replaced by a differential inclusion.The sliding set may be decomposed into the attractive sliding mode, for which the vector fields point toward each other, and the repulsive sliding mode, for which the vector fields point away.On the attractive sliding mode, the solution of the differential inclusion is uniquely given by a sliding solution along .On the repulsive sliding mode, the solution is non-unique as it may slide along or leave on either side.Various types of numerical methods exist to simulate Filippov systems (or differential inclusions).Time-stepping methods [27,35] directly discretize the differential inclusion using an implicit scheme.Eventdriven methods integrate an ODE up to the switching manifold, detect a possible sliding mode and then solve a differential algebraic equation on the sliding mode.
In the case of a CPWL system, it can easily be shown that trajectories can only cross and that therefore sliding behavior is excluded.However, for a discontinuous PWL system such as the reduced dynamics obtained in (20), sliding motion can occur, which presents a dynamical behavior that cannot be encountered in the original CPWL full system (5).In order to simulate the discontinuous system (20) numerically, a specialized event-driven simulation method capable of accurately localizing the switching points and handling sliding solutions is needed.In this work, the MATLAB package DISODE45, developed by Calvo et.al [36] based on an adaptive Runge-Kutta scheme to solve Filippov systems is used.In order to show that sliding motion can indeed occur in the discontinuous reduced system, as opposed to the original full system, the time histories of the slow variables are compared in Fig. 2 for the parameter set m = k = R = 1, c + = 1.8, c − = 0.01 and ε = 0.01 under an external excitation with the frequency = 0.94.The dotted black line representing the solution of the full system is not matched by the solution of the discontinuous reduced system shown in yellow.In Fig. 2b, the sliding behavior in the discontinuous system is apparent at the time instants where x 2 = 0 and the solution slides along the switching manifold for a finite time interval, before leaving it and crossing to the other linear region.The blue and orange lines obtained by the continuous reduced models, which are derived in Sect.4, capture the behavior of the full system accurately, avoiding the sliding behavior.Obviously, for other parameter values, the discontinuous reduced system can also exhibit direct crossing behavior such as the CPWL full system.However, sliding cannot a priori be excluded without ensuring the continuity of the reduced model.Hence, this motivates the construction of a reduced system which retains the continuity property of the full system.

Model reduction by a continuous matching approach
This section presents a general treatment of the critical and slow manifolds of CPWL systems with regard to their continuity properties.First, the critical manifold obtained as a concatenation of two critical "halfmanifolds" is shown to be continuous at the switching manifold for arbitrary CPWL systems.Then, a continuous matching approach is proposed in order to circumvent the discontinuity problem arising in the approximation of the slow locally invariant half-manifolds.Without loss of generality, we consider nonautonomous slow-fast CPWL systems with a single switching hyperplane described by Eq. ( 5), where 0 < ε 1 is the small perturbation parameter, dt denotes the derivative with respect to the slow time scale t and h is a linear scalar switching function given by h(x, y) = ẽT j x, with ∇h nonzero and h(0, 0) = 0 where x ∈ R n is the state vector and ẽ j is the jth basis vector of R n .Since we assume that the system is CPWL, the functions f ± and g ± are linear with respect to x and y and continuous at the switching manifold.The state vector x ∈ R n contains all state variables, which can be divided into a set of slow vari- and n = s + f .We note that such a decomposition of the state variables into fast and slow states requires a special choice of coordinates and is therefore not necessarily at hand for a general slow-fast system.For the interested reader, the book by Wechselberger [30] provides a comprehensive review of geometric singular perturbation theory in the more general coordinate-free setup.

Continuity property of the critical manifold
In the following, we assume that the switching manifold is not tangent to all fast directions, that is, ∇h is nonzero in at least one of its last f components.This yields a more general switching condition compared to the case where the switching does not involve the fast variables.By reordering the components of y, one can assume that ∂h ∂ y 1 = 0. Following [29], to rewrite system (5) in a simpler form, an invertible transformation (x, y) → (x, y) is introduced, with the same slow coordinates x = x and a new set of fast coordinates y = y1 • • • y f T , defined as To check that the new state variable y1 is indeed fast, its dynamics can be easily obtained as with ∂h ∂ y 1 = 0. Therefore, system (5) can be rewritten using the new set of coordinates and by replacing the switching condition h(x, y) ≷ 0 simply by y 1 ≷ 0 as where the ( •) are dropped for simplicity and g represents the fast dynamics in the new coordinates.For a CPWL system of the form (26), the linear functions f ± and g ± read as with constant matrices depending on ε, with D ± assumed to have full rank.The vectors r 1 (t) and r 2 (t) contain the time-dependent forcing terms.
Proposition 1 Consider a slow-fast CPWL system in the form (26) governed by the slow and fast dynamics given in (27), with the matrices D ± being full rank.Let be the critical half-manifolds obtained for the limit ε = 0 from the ± subsystems, with e 1 being the first vector of the standard basis of R f .Then, the critical manifold M c defined as the union of the critical halfmanifolds The formal proof of this proposition is given in Appendix A. From an intuitive point of view, this result is not surprising.For a CPWL system written in matrix form as ẋ = Ã± x+r(t), the system matrices describing both subsystems are identical with the exception of the column that contains the coefficients of the variable y 1 , which in turn dictates the switching condition.At the switching manifold, y 1 = 0 holds and the vector fields of both subsystems evaluated at the switching manifold are equal.Henceforth, the dynamics of system (26) for the limit ε = 0 takes place on a critical "manifold" consisting of two linear parts meeting continuously along a kink on the switching manifold.In Fig. 3a, the continu-Fig.3 Sketches of the critical manifold M c and the slow half-planes M ± s ous critical manifold with a kink is depicted graphically in a three-dimensional space.Assuming the attractivity of the critical half-manifolds in their corresponding subregions, it follows that the dynamics of each subsystem can be reduced to a slow locally invariant manifold.This yields two slow half-manifolds for each region of the system, similar to the critical half-manifolds.However, the continuity property characterizing the critical manifold is automatically lost for approximations of order O(ε) and higher, as shown in the example of the CPWL quarter car model.The resulting slow halfmanifolds intersect the switching manifold along different hypersurfaces and only meet at the origin.This is shown in Fig. 3b.Therefore, reducing the dynamics of each subsystem to an approximation of the slow manifold without a proper matching procedure at the discontinuity would yield a reduced model containing a jump in the vector field at the switching manifold (i.e., a Filippov system).This, in return, results in a general PWL reduced dynamics, capable of exhibiting more complex dynamical behavior than the original CPWL system, such as sliding behavior as was shown for the example in the previous section.Hence, the aim is to circumvent the discontinuity at the switching by continuously matching the slow half-manifolds up to their intersection, in order to obtain a slow "manifold" with a kink analog to the critical manifold.The resulting reduced dynamics must be continuous to avoid the possibility of sliding motion along the switching manifold, thereby preserving the nature of the original system.In the next subsection, the continuous matching approach is developed for arbitrary problems.

Continuous reduced slow dynamics
For a non-autonomous CPWL system of the form ( 5) with a single switching manifold , the state space can be divided into two subregions, which we refer to by the + and − region, where h(x, y) ≥ 0 and h(x, y) < 0 hold, respectively.In each of these regions, we expect the trajectories to converge toward the corresponding slow locally invariant manifold.In order to obtain a continuous reduced system, the two linear locally invariant slow manifolds M ± s (t) can be continued up to their intersection M + s (t) ∩ M − s (t).By construction, these manifolds are the graphs of the vectorvalued linear functions which can be obtained in closed form from Eqs. ( 14) and ( 15) following the same procedure described in Sect. 2. For simplicity, we assume s = 2, i.e., there are two slow variables x 1 and x 2 , as was the case for the quarter car model with bilinear damping, whereby a generalization for arbitrary s is straightforward.For the purpose of illustrating the continuation approach, it is helpful to consider a single fast variable y i and its approximation on the slow half-manifolds given by the scalar functions h ± s,i (x, t), which allows to handle planes in R 3 , which we denote by M ± s,i (t).The manifolds are planes since ∇h ± s,i are constants, i.e., h ± s,i (x, t) = (n ± s,i ) T x + ν ± (t), where ν ± (t) are timedependent terms shifting the planes in time and n ± s,i are the corresponding normal vectors.The planes M ± s,i (t) cross the (x 1 , x 2 )-plane along two zero contour lines given by Fig. 4 Manifold continuation up to the intersection line C s,i (t).Projection to the (x 1 , x 2 )-plane showing the different regions obtained by the inequality signs of h ± s,i The aim is to obtain a rule defining the approximation of y i using the expressions h + s,i and h − s,i , which is continuous in the slow coordinates x.For an arbitrary pair of linear functions h ± s,i that map the slow variables x to the fast variable y i , the space of slow coordinates can be divided into the following six regions, as depicted in Fig. 4, depending on the inequality signs Each edge between these regions corresponds to an equality: The edge between regions (i) and (ii) as well as between (iv) and (v) is given by h + s,i = h − s,i and is denoted by C s,i , depicted in red.The zero-level h + s,i = 0 corresponds to the edge between (i) and (vi) as well as between (iii) and (iv), denoted by C + s,i and depicted in green.Lastly, the zero-level h − s,i = 0 gives the edge C − s,i shown as an orange line between regions (ii) and (iii) as well as (v) and (vi).A continuous rule defining the approximation of y i using the expressions h ± s,i must be given as a quantity that remains unchanged across an edge or one that can only change continuously between these expressions.Therefore, a switch from the expression h + s,i to h − s,i or vice versa is only allowed across the edge C s,i where the equality h + s,i = h − s,i holds, which presents the projection of the intersection line M + s,i (t) ∩ M − s,i (t) on the (x 1 , x 2 )-plane and is used as a new switching condition in the reduced slow dynamics.As a consequence, the choice of h +/− s,i in region (i) forces the same choice h +/− s,i in regions (vi) and (v).The same holds on the other side of C s,i .Thus, the only possible continuous rules defining y i using h ± s,i are given as follows: 1. Choose y i = h + s,i everywhere.

Choose y i = h −
s,i everywhere.

Choose y i = h +
s,i in regions (i), (vi) and (v) and choose y i = h − s,i in regions (ii), (iii) and (iv).This rule can be expressed by 4. Choose y i = h − s,i in regions (i), (vi) and (v) and choose y i = h + s,i in regions (ii), (iii) and (iv).This rule can be expressed by Obviously, the first two options correspond to a reduction to the slow invariant manifold obtained from one subsystem and ignore the piecewise nature of the CPWL system.Therefore, and in order to account for the switching, one is left with options 3 and 4 using the min and max functions, respectively.Graphically, options 3 and 4 can be regarded as extending the slow half-plane M + s,i (t) across the switching plane to continuously meet M − s,i (t) and vice versa.The resulting continuous approximations of the fast variables are then substituted in the slow dynamics f ± , thus yielding a continuous ROM without jump at the switching manifold.In the following two sections, we apply the continuous matching approach to two examples of harmonically forced CPWL systems.The first example is the quarter car model with bilinear damping from Sect. 2, where the switching manifold is given in terms of a fast variable.In the second example, the case of bilinear stiffness is investigated with a simpler switching manifold given in terms of a slow variable.The derivation of continuous ROMs as well as the obtained numerical results including frequency-response curves are discussed.

Application: quarter car model with bilinear damping
As illustrated in Sect.2, the reduced-order model obtained from the critical dynamics of the quarter car model is continuous but does not include the explicit time dependence, and hence fails to capture the nonautonomous character of the original system.Considering an ε order approximation of the slow half-manifolds results in a discontinuous ROM which may exhibit sliding behavior.Aiming to avoid this sliding scenario and obtain a CPWL ROM, the continuous matching approach proposed in Sect. 3 is applied to the quarter car model with bilinear damping.

Continuous ROM with bilinear damping
Going back to the approximations of the slow locally invariant half-manifolds given in ( 18) and ( 19), a continuous rule for the fast variables following the approach from Sect. 3 can be obtained using the min and max functions.In order to estimate which of those two options would yield the better approximation of the fast variables, one could investigate the far field behavior by considering the expressions of h ± s obtained from the critical manifold ( 9) and the first-order terms from (17).In the following and without loss of generality, we assume that the damping coefficients fulfill the condition c + < c − .Since the switching manifold is given by y 2 = 0, we first consider the continuous matching of the approximations h ± s,2 .For ε = 0, the min and max functions yield where the full expressions of are obtained directly using Eqs.( 9) and (17).In view of these full expressions and for finite x 1 , small but nonzero ε and sufficiently large-in-magnitude and positive y 2 , both h + s,2 and h − s,2 imply that x 2 is large in magnitude and positive, such that x 2 ≥ −k c + +c − x 1 .In this case, Eqs.(34) yield Note that these expressions are based on the assumption that y 2 is large and positive.In this region of the state space, the trajectories are attracted to the slow manifold of the + subsystem and therefore h + s,2 should be the better approximation, which is given by the min function.Analogously, a large-in-magnitude and negative y 2 implies also a large-in-magnitude and negative x 2 , and thus, one obtains of which only the expression obtained from the min function agrees with the expected far field behavior.
Regarding the remaining fast variable y 1 , similar arguments can be used to a priori determine which of the approximations is more suitable for the reconstruction of y 1 .For simplicity, taking the limit ε = 0 in h ± s,1 leads to Therefore, one obtains Since for y 2 > 0 and thus x 2 > 0 the trajectories are attracted to the manifold M + s (t) from the + subsystem, one can expect the approximation using min to be a more suitable choice than max, also for the reconstruction of the fast variable y 1 .
In summary, the resulting continuous reduced dynamics can be obtained in an explicit form, which reads as with for i = 1, 2 and the varying damping coefficient defined by Table 1 Choice combinations for h s,i = min(h In view of Eqs. ( 39), ( 40) and ( 41), the continuous reduced model is a CPWL system with eight linear subsystems.The continuous matching of M ± s,1 (t) and M ± s,2 (t) yields two different possibilities each, depending on the use of max or min.For instance, in the case of a reduction using the max operator, one would substitute h + s,1 in the reduced dynamics (39) as long as h + s,1 > h − s,1 holds, and correspondingly, h − s,1 is used wherever h − s,1 > h + s,1 holds.This works analogously for h s,2 .However, a sign change in the term h s,2 implies automatically an additional switch in the damping coefficient c(h s,2 ).Therefore, one ends up with eight different subsystems, which are linear due to the linearity of the approximations h ± s,1 and h ± s,2 .The various sign combinations and the corresponding choices for h s1,2 and for the damping coefficient c(h s,2 ) to be substituted in (39) are summarized in Tables 1 and 2 for the cases where the approximation of both fast variables is obtained using the min and max operators, respectively.

Numerical results
To illustrate the accuracy of our continuous matching approach in approximating the steady-state response of the quarter car model with bilinear damping, we Fig. 5 Instantaneous projection of the four-dimensional dynamics in the three-dimensional space (x 1 , x 2 , y 1 ) of the computed PWL time-dependent manifolds (here depicted for t = t * ) for ε = 0.01 and the excitation frequency = 0.94 consider the numerical integration of the full CPWL model ( 5) and the CPWL reduced models given by Eqs. ( 39) and (40) using the following values of the system parameters First, the small parameter is chosen as ε = 0.01.As illustrated in Fig. 5, a trajectory of the full system initialized on a random point of the state space converges toward the computed PWL manifolds with a kink.Shown is the complete history of the trajectory as dashed line in the three-dimensional space (x 1 , x 2 , y 1 ), as well as the solution at a particular time instant t * as a black dot.The instantaneous projections of the PWL manifolds obtained at t * are shown in pink and sky blue for the reductions using the max and min operators, respectively.These projections are obviously time-dependent and the steady-state solution of the full system lies in their immediate vicinity at every time instant.The dynamics of the CPWL ROM obtained by (39) converges to a periodic solution (red line) which approximates the steady state of the full system very accurately.The difference between the periodic orbits obtained using max and min was found to be negligible and only one solution is illustrated for clarity (red line).This can be explained as follows.For ε small enough, the orientations of the different slow half-manifolds M ± s,2 (t) become very similar and their relative angle vanishes for ε = 0, which can be easily seen by evaluating their normal vectors.Since the approximation h s,1 is multiplied by ε in the slow dynamics, the kink observed in M ± s,2 (t) is the determining factor in the vanishing difference between the reductions using max and min.Concerning the fast variables y 1 and y 2 , the exponentially fast convergence of the solution of the full system toward the reduced dynamics is shown in Fig. 6.
Next, a higher value of the small parameter ε = 0.1 is investigated, for which the forced-response curves (FRC) around the harmonic resonance frequency using the full model and the continuous reduced models have been obtained.The results are shown in Fig. 7a and 7b, where the maximal displacement and velocity are plotted against the excitation frequency, respectively.
Figures 7a and 7b show that the reduced models are able to capture the forced response of the full system with a great accuracy, even with a significant nonlinearity caused by the difference between the values of the damping coefficients c + = 0.1 and c − = 0.5 and a relatively high perturbation parameter ε = 0.1.The continuous slow reduced dynamics, both using the functions min (blue) and max (orange), show a good agreement with the full system.Looking at the peak amplitudes in the FRCs and the time histories at the resonance in Fig. 8a and 8b, it is apparent that the con-tinuous reduced model using the max function gives a slightly better agreement with the full system in approximating the slow variables around the main resonance frequency.The inverse holds for the approximation of the fast variables, where min gives the better overall reconstruction of y 1 and y 2 , as expected from the discussion in the previous subsection.As of yet, there is no a priori criterion for the difference in the accuracy between max and min depending on whether the slow or fast variables are approximated.However, it should be noted that this trend appears only if the perturbation parameter ε is rather high, which automatically contradicts the assumed slow-fast nature of the system.By taking a closer look at the error dynamics between the two ROMs and the full dynamics, this inverse trend may be linked to a resonance behavior in the error dynamics.Nevertheless, the error dynamics has to be investigated closely in an effort to determine a priori the more suitable choice for the ROM, which remains to be done in a future work.
Another feature that stands out from Fig. 7a and 7b is that the system with a bilinear positive damping dis-123 Fig. 8 Time histories of the slow variables for = 0.94 and ε = 0.1 (the same legend as in Fig. 7) plays a similar frequency response as a linear system.This fact can be explained by the convergence property of the system.In the work of Pavlov et al. [37], it was shown that for a CPWL system with a time-dependent input, the existence of a common quadratic Lyapunov function is equivalent to a quadratic convergence of the system.The convergence property means that all solutions of the system converge to a unique steady-state which is period-1, i.e., the system exhibits the same period time as the input r (t).Therefore, this property implies the absence of coexistence of periodic solutions or other limit sets, which in turn excludes any nonlinear sub-or superharmonic resonance.
By applying the results from [37], a linear matrix inequality can be solved to obtain a common quadratic Lyapunov function, which was found for every combination of positive damping coefficients.
Thus far, it has been shown that the proposed continuous matching of the linear locally invariant manifolds captures the periodic behavior of the full CPWL system.Two admissible possibilities to continuously match the linear locally invariant slow manifolds have been compared.The results suggest that for ε small enough (0.01), the continuous matching approach gives a good approximation of the full system behavior.For a high value of the perturbation parameter ε = 0.1 and for the chosen parameter set, the CPWL ROM using the max function results in a slightly better approximation of the full system with respect to the peak amplitudes of the slow variables compared with the reduction using the min function, which in turn yields a better reconstruction of the fast variables.Further research is therefore required to better understand the mechanisms behind this effect and determine a priori which slow subspaces should be chosen in the different subregions of the state space in order to obtain the best reduced CPWL model with respect to the slow variables.

Application: forced oscillator with bilinear stiffness
In this section, another variation of the quarter car model with a piecewise linear stiffness is studied.
In contrast to the case with bilinear damping, the convergence property is lost when dealing with a PWL stiffness, since one cannot find a common quadratic Lyapunov function which holds for both subsystems for all possible parameter values [37].Hence, by introducing this type of static PWL nonlinearity in the system, bifurcation scenarios become possible, leading to the appearance of more complex nonlinear phenomena in the system response, such as super-and subharmonic resonances.To illustrate the ability of our reduction approach to capture such nonlinear behaviors, the system shown in Fig. 10 is investigated by means of singular perturbation theory.The system consists of two masses m and εm, connected by a linear spring damper element with stiffness k and damping coefficient c, as well as a unilateral spring with stiffness k, which becomes active only if the relative displacement x − z is negative.The lower mass is connected to the road by a linear spring with the stiffness k ε and the road profile is assumed as a harmonic excitation r (t) = R sin( t).By using the set of with the slow and fast variables defined by x = x 1 x 2 T and y = y 1 y 2 T , respectively, and the continuous lin-ear functions f ± and g ± given by For consistency with Sect.2, we keep the notation f ± (x, y, t; ε) and g ± (x, y, t; ε) even though the functions do not have to depend on all the variables.The switching stiffness coefficient k ± is defined as The switching manifold is denoted by = {x ∈ R 4 : x 1 = 0} with x = x 1 x 2 y 1 y 2 T .Analogous to the previous section, applying the procedure described in Sect. 2 yields the critical manifold M c = {x ∈ R 4 : y = h ± c (x, t)} with the expressions In this case, the switching condition is explicitly given in terms of the slow variable x 1 , which simplifies the following analysis.Again, the continuity and stability of the critical manifold can be easily checked and the reduced dynamics on M c is described by Using Eq. ( 15), one can obtain an explicit expression for the first-order terms of the slow manifolds for both the + and − subsystems, which read as The approximation of the slow manifolds is then obtained by Here, the switching condition depends explicitly on x 1 .
Hence, the problem with regions where both validity conditions are simultaneously valid or invalid, as described in Sect.2, does not occur.The state space can directly be split into two subregions given by x 1 ≥ 0 and x 1 < 0 and the expressions h ± s are substituted correspondingly to obtain the following discontinuous reduced model For the purpose of obtaining a continuous reduced model, the procedure described in Sect. 3 can be applied.Note that the slow dynamics f ± does not depend on the fast variable y 1 .Therefore, only h ± s,2 are needed for the reduced model and h ± s,1 are only used to reconstruct y 1 .To determine a priori which of the alternatives is better suited for the approximations of the fast variables, the desired far field behavior can be investigated similarly to Sect. 4. Since we assume k − = 2k + , the limit ε = 0 implies that min of which only the min function is consistent with the trajectories being attracted to the slow locally invariant manifolds of the corresponding ± subsystems.Therefore, the min function is expected to deliver a better reconstruction of y 1 than max.However, for nonzero ε, the difference between the approximations h ± s,2 reads Therefore, for a large-in-magnitude and positive x 1 , it follows that The opposite results for a large-in-magnitude and negative x 1 .In both cases, the approximation using max agrees with the desired far field behavior and is therefore expected to yield the better reconstruction of the fast variable y 2 .In summary, the resulting continuous reduced dynamics can be obtained in an explicit form, which reads as ẋ = f cont (x, t; ε) and the varying stiffness coefficient given in (46).For a numerical illustration of the accuracy of the different reduced models, we choose the parameter values A brute force frequency sweep-up was performed to obtain the FRC shown in Fig. 11a corresponding to the maximum amplitude of the relative displacement x 1 during one period of oscillation.Figure 11b shows the errors in the approximation of the maximal amplitudes which are measured for every excitation frequency in the sweep-up as follows: where u( ) and ũ( ) are the maximal x 1 amplitudes over one period obtained from the full system and the ROM, respectively, under the harmonic excitation of frequency .The FRC of the full system is shown as a black line.Overall, for ε = 0.01, which can be considered as small enough, both CPWL ROMs capture the nonlinear behavior of the full system.As shown in Fig. 11b, the reduction using the slow half-manifolds yields a more accurate approximation than the reduction along the critical manifold, which is not surprising.The main harmonic resonance (around = 1.09 [rad/s]) and the super-and subharmonic resonance peaks around = 0.58 [rad/s] and = 2.31 [rad/s], respectively, are approximated accurately by the CPWL ROMs using max and min.The small peaks in the errors shown in Fig. 11b at = 2.11 [rad/s] are due to a slightly shifted onset of the subharmonic solution branch in the ROMs compared with the full system.Otherwise, the errors from the CPWL ROMs are at least of magnitude around O(ε).Additionally, the time histories at the peak of the subharmonic resonance are shown in Fig. 12, where the harmonic excitation is illustrated for comparison.

Conclusions
In this paper, a class of slow-fast harmonically forced oscillators with piecewise linear nonlinearities was studied using singular perturbation theory.It was observed that the proposed continuous matching of the linear locally invariant slow manifolds of both subsystems is able to approximate the behavior of the full system with high accuracy for a frequency range around the main harmonic if the perturbation parameter is small enough.Regarding the forced system with bilinear damping, the proposed technique captures the behavior of the full system accurately and circumvents the numerical challenges, which may arise if one takes a switching condition that changes the continuous nature of the original system and yields a discontinuous reduced model.Due to the convergence property, it was shown that the PWL system with positive bilinear damping admits frequency-response curves similar to a linear system.These curves were well approximated by the proposed approach, even for a relatively high nonlinearity resulting from a significant difference between the damping coefficients.For a similar harmonically forced slow-fast oscillator with static PWL nonlinearity instead of the bilinear damping, the system behavior becomes more complex due to the loss of the convergence property and the existence of nonlinear phenomena, such as super-and subharmonic resonances, becomes possible.These nonlinear resonances were accurately captured by the proposed reduction approach for the frequency range around the main resonance.
This study provides insights for the extension of the singular perturbation theory to nonsmooth systems with PWL nonlinearities.The continuous matching of the slow manifolds can be achieved by coupling the fast variables to the slow variables using the min and max functions.A future study is therefore required to establish an a priori criterion as to which continuous matching procedure is better suited for the purpose of model reduction.Moreover, the proposed reduction was shown to give a good approximation of the system behavior around the main resonance frequency.Higher excitation frequencies would therefore interfere with the slow-fast decomposition of the system variables, which could be investigated in future studies.Further work may also focus on extending the proposed reduction method to systems described by (partial) differential inclusions and to compare our approach with existing results such as in [38] and in [39].
Author contributions Both authors contributed to the conceptualization, design and methodology.YK carried out the formal analysis.YK wrote the paper.Both authors commented on previous versions of the manuscript.RL was responsible for funding acquisition and supervision.Both authors read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data Availability Statement
All numerical data in this work have been generated from the considered system equations with the approaches described in this work and the cited works, using MATLAB standard methods.Therefore, it is possible to completely reproduce the data from the information given in this work.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A
Proof of Proposition 1: Since the switching manifold is defined by the condition y 1 = 0, the continuity property implies that the matrices B ± differ only in their first column and the same holds for the matrices D ± .The critical half-manifolds M ± c defined in (28) are obtained from the isolated solutions of the algebraic Eq. (8b), which read as

Fig. 1
Fig. 1 Quarter car model of automotive suspension with bilinear damping

Fig. 2
Fig. 2 Time histories showing the sliding behavior in the discontinuous ROM for = 0.94, c + = 0.01, c − = 1.8, ε = 0.01, m = k = R = 1.The time history of the discontinuous ROM, shown in yellow, does not accurately approximate the

Fig. 6 Fig. 7
Fig.6 Fast convergence of fast coordinates y 1 and y 2 along the full trajectory toward the same coordinates along the continuous PWL reduced systems with ε = 0.01, = 0.94

Fig. 9 Fig. 10
Fig.9 Time histories of the fast variables for = 0.94 and ε = 0.1 (the same legend as in Fig.7)

Fig. 11
Fig.11Frequency-response curve of the harmonically forced 2-DOF oscillator with bilinear stiffness with the parameter setting (59)

Fig. 12
Fig. 12 Time histories at the peak of the subharmonic resonance at = 2.31

Table 2
Choice combinations for h s,i