Surrogate model approach for investigating the stability of a friction-induced oscillator of Duffing's type
NNoname manuscript No. (will be inserted by the editor)
Surrogate model approach for investigating thestability of a friction-induced oscillator of Duffing’stype
Jan N. Fuhg · Am´elie Fau
Received: date / Accepted: date
Abstract
Parametric studies for dynamic systems are of high interest to de-tect instability domains. This prediction can be demanding as it requires arefined exploration of the parametric space due to the disrupted mechanicalbehavior. In this paper, an efficient surrogate strategy is proposed to inves-tigate the behavior of an oscillator of Duffing’s type in combination with anelasto-plastic friction force model. Relevant quantities of interest are discussed.Sticking time is considered using a machine learning technique based on Gaus-sian processes called kriging. The largest Lyapunov exponent is proposed asan efficient indicator of non-regular behavior. This indicator is estimated usinga perturbation method. A dedicated adaptive kriging strategy for classifica-tion called MiVor is utilized and appears to be highly proficient in order todetect instabilities over the parametric space and can furthermore be used forcomplex response surfaces in multi-dimensional parametric domains.
Keywords
Dry friction · Lyapunov exponents · Non-smooth system · Surrogate model · Adaptive sampling · Machine learning · Stick-slipinstability
Jan N. Fuhg (Corresponding author)Institute of continuum mechanicsLeibniz Universit¨at HannoverAppelstraße 1130167 Hannover, GermanyE-mail: [email protected].: +49 (0)511.7 62 - 22 85Fax: +49 (0)511.7 62 - 54 96Am´elie FauInstitute of mechanics and computational mechanicsLeibniz Universit¨at HannoverAppelstraße 9A30167 Hannover, Germany a r X i v : . [ ee ss . S Y ] J u l Jan N. Fuhg, Am´elie Fau
Stick-slip instability can be observed when two rigid bodies in contact slide oneach other at low relative velocity. Then, intermittent vibratory effects mayoccur in form of variations of the frictional force and relative sliding velocity,with a characteristic sawtooth time-displacement curve, corresponding to stickand slip phases [1]. These dynamic instabilities, also called stick-slip motionresult in self-excited vibrations and drastic decreases in performance of somemachine parts. Oscillating systems excited by dry friction are frequently en-countered in many practical applications, including for instance brake systems[2, 3], hydraulic cylinders [4], gears [5] or bearing [6]. Therefore the exten-sive analysis of these non-smooth dynamic systems to detect instabilities is acrucial point in engineering design.Linear spring mass systems placed on a moving belt, commonly referred to asMass-on-Belt (MoB) systems, have been generally utilized as simple mathe-matical models to investigate stick-slip system, see e.g. [7, 8, 9, 10]. To providebetter accuracy, nonlinear mathematical representations have recently beenconsidered [11]. For instance, the bifurcation behavior of the non-linear MoBsystem has been investigated in [12, 13]. A Duffing’s type oscillator has beenstudied analytically in [11] in order to obtain expressions for stick-slip andpure-slip vibration amplitudes and frequencies. An estimation method for thespectrum of Lyapunov Exponents (LE)s proposed by [14] is employed in [15]to analyze the stability of a discontinuous MoB system.To avoid computing and analyzing the behavior in the time domain, some indi-cators of instabilities have been proposed in the literature, such as Kolmogoroventropy [16], correlation dimension [17] and sticking time [18]. LEs have beenstated by some authors as the most significant indicator for chaotic motionsuch as in [19]. They quantify the rate of separation of infinitesimally close tra-jectories and the sensitive dependence on initial conditions exhibited by thesystem. In practical applications the determination of the Largest LyapunovExponent (LLE) is sufficient. If the number is non-positive, then the system isstable. On the contrary, a positive LLE indicates a system with chaotic behav-ior. Unfortunately, analytic determination of LEs is possible only for simplelinear dynamic systems, see [14], and numerical determination of their value isa challenge. In recent years several numerical strategies have been suggestedfor estimating LEs in case of non-smooth dynamic systems. A thorough re-view can be found in [20]. Recently an approximate numerical method basedon the estimation of the Jacobian matrix by perturbation method of the initialorthogonal vectors involving an Euler forward scheme has been proposed [14].In addition, engineers want to prevent instabilities for a large scope of param-eters, or scenarios which may be encountered during the working life of themechanical systems. However the numerical analysis of the LLE in complexsystems is computationally expensive, limiting the design space explorationand the optimization of a system at hand. Therefore statistical approximationsor surrogate models emulating the LLE over a parametric input space appearvery interesting to estimate the response on vast parametric space at low cost. riging approach to investigate stick-slip instability 3
Surrogate model approaches such as polynomial response, see e.g. [21], krigingalso known as Gaussian process regression [22, 23], support vector regression[24] or radial basis function models [25] have been developed and continuouslyimproved in recent years. Kriging is an accurate interpolative Bayesian surro-gate modeling technique, which is stated as the most intensively investigatedsurrogate approach [26]. Originally proposed by the geostatistician Krige [27],it has been utilized for deterministic [28] and random simulations models [29].The idea of a kriging model is to provide an estimation of the response ofthe system on the whole parametric domain from the restricted available in-formation provided by some observations. Two perspectives can be explored.One is generating an accurate surrogate model from the information given apriori . The second is to inform about the localization of crucial areas in theparametric space for the problem of interest, i.e. subdomains in which more in-formation should be gained, if further exploration is possible. In light of recentcontributions proposed for dynamical analysis and machine learning tools, thegoal of this contribution is to propose surrogate strategies to investigate theoscillatory behavior and efficiently detect the risk of instabilities.Thus, in this paper surrogate model construction to identify the parametricsubdomains where dynamic instabilities may occur is investigated. Pertinentinstability indicators are discussed, different metamodeling techniques are con-templated. Among them, kriging appears the most powerful for the problem ofinterest. Numerical challenges and capacities offered by kriging for the estima-tion of the sticking time and LLE analysis in a Duffing’s type oscillator withan elasto-plastic friction force model are exposed. More particularly, variousadaptive schemes are evaluated in this context.The paper is organized as follows. In section 2 a dynamic problem non-linearcontact behavior is introduced, as well as the accompanying challenges dueto the possible unstable behavior. Two possible indicators proposed by theliterature to proficiently investigate the system behavior are exposed in section3. A brief introduction to LE as well as most recent works which allow toestimate the LLE for discontinuous systems in a robust numerical frameworkis also presented. In section 4, an overview of the kriging approach is given.Section 5 represents the core of this contribution, with a detailed review of theability and challenge of the surrogate model approach to investigate unstabledynamic behaviors.
The mechanical reference problem is a Duffing’s type oscillator with a dampingterm, as illustrated in Figure 1. A rigid body of mass M is placed on a movingbelt with a constant velocity V . The displacement of the body over time t isdenoted by X ( t ). The body’s movement is restricted by a nonlinear spring ofstiffness K X + K and a dashpot with damping coefficient D parallel to thespring. The relative velocity between the belt and the body is denoted by V R ( t ),which equals to V − ˙ X ( t ) with ˙ X ( t ) representing the first derivative of the Jan N. Fuhg, Am´elie Fau body’s displacement with respect to time, i.e. its velocity. A time-dependentharmonic force U ( t ) = U sin( Ωt ) with amplitude U and angular frequency Ω as well as a constant normal load N is applied to the mass. Fig. 1 – Scheme of the analyzed nonlinear mass-on-belt system.The equation of motion of the system reads M ¨ X ( t ) = − D ˙ X ( t ) − K X ( t ) − K X ( t ) + F R ( V R ) + U sin( Ωt ) , (1)where ¨ X ( t ) is the second derivative of the body’s displacement with respectto time, i.e. the body acceleration and F R ( V R ) = N f R ( V R ) is the frictionforce. f R ( V R ) is called the friction function which denotes the friction forceper unit of normal load; it is a function chosen by the user. An overview of thecommon choices proposed in the literature is given in [30]. Generally, frictionforce models are classified either as static or dynamic [31], where in a dynamicmodel the friction force does not only depend on the relative velocity betweenthe bodies but also on other state variables. A review of different friction forcemodels for dynamic analysis has been presented in [32]. Here, a dynamic modelproposed in [33] is employed. It belongs to the family of elasto-plastic models,which were initially introduced by [34] to model the constitutive behaviorof deformable bodies. Internal variables are used to represent the non-linearcontact behavior. The displacement X of the body is assumed to be the sumof an elastic, i.e. memoryless, contribution denoted by z and a plastic, and sohistory-dependent, part w such that X = z + w . (2)A physical analogy of this model is depicted in Figure 2. During the stickingphase, the plastic displacement is constant, whereas while slipping the elasticcomponent is assumed to be fixed.The friction model, as illustrated in Figure 3, is based on the idea that anydeformation of the asperities is associated with a corresponding friction force F R . By analogy with plasticity, the history of the model is characterized by riging approach to investigate stick-slip instability 5 Fig. 2 – Physical analogy of the elasto-plastic model representing thecontact behavior - the block displacement X over the surface is broken downinto an elastic component z and a plastic part w , modified from [35]the internal variable z , which represents the bristle deflection. The frictionfunction is defined as f R ( V R , z ) = σ z + σ ˙ z + σ V R , (3)which can be compared to the force model of the classical LuGre model [36].The variable σ denotes the bristle stiffness, σ is the average bristle dampingcoefficient and σ is a viscous component of the friction force. Fig. 3 – Basic concept of friction force formulation. Asperities in the contactsurface are modeled as elastic bristles with damping. For simplicity theLuGre model considers an average bristle deflection z for the determinationof the friction force. Jan N. Fuhg, Am´elie Fau
Consider the following definitionsgn ( V R ) = (cid:40) V R (cid:107) V R (cid:107) , if V R (cid:54) = 0 , , if V R = 0 , (4)the velocity of bristle deflection ˙ z has been defined by [35] as˙ z = (cid:18) − α ( z, V R ) σ g ( V R ) z · sgn ( V R ) (cid:19) V R , (5)where α ( z, V R ) is a function, which is utilized to capture stiction effect. Itdepends on the relative velocity see e.g. [37], since it defines the elastic defor-mation until the breakaway force of the system is exceeded, i.e. α ( z, V R ) = (cid:40) α ( z ) , if V R · z ≥ , if V R · z < α ( z ) = , if z ≤ z ba , (cid:32) sin (cid:32) π (cid:107) z (cid:107) − z max + z ba z max − z ba (cid:33) + 1 (cid:33) , if z ba < (cid:107) z (cid:107) < z max , , if z max ≤ (cid:107) z (cid:107) . (7)The maximum bristle deflection is denoted by z max , and the bristle deflectionfor breakaway condition as z ba . The function g ( V R ) models the Stribeck effect,which refers to an increase of the friction coefficient when the relative contactvelocity increases [38]. It is given by g ( V R ) = N (cid:18) µ k + ( µ s − µ k ) exp (cid:18) − V R V S (cid:19)(cid:19) , (8)where V S is a parameter called the characteristic Stribeck velocity, and µ s and µ k denote the static and kinetic friction coefficients, respectively. Therefore,the complete equation of motion of the system in state-space form with themodified elasto-plastic friction model reads dXdt = ˙ X, d Xdt = − DM ˙ X − K M X − K M X + N M f R ( V R ) + U M sin( Ωt ) , dzdt = (cid:18) − α ( z, V R ) σ g ( V R ) z · sgn ( V R ) (cid:19) V R , (9)which can be solved, e.g. via numerical integration.From the equations of motion, the dynamic behavior of the system can beanalyzed. For instance, the time-response behavior plotted in Figure 4a showsthe cyclic unstable behavior in time of the mass via the typical cyclic evolutionof its displacement and its velocity. The cyclic unstable behavior can also be riging approach to investigate stick-slip instability 7 (a) Time response (b)
Phase portrait
Fig. 4 – Time response and phase portrait for the nonlinear mass-on-beltsystem with X in m, ˙ X in m . s − , t in s, Ω = 0 . . s − , M = 1 kg, V = 0 . . s − , D = 0 . . m − , K = 1 N . m − , K = 0 . . m − , µ s = 0 . µ k = 0 . V s = 0 . . s − , U = 0 . N = 1 . σ = 100 . . m − , σ = 10 . . m − and σ = 0 . . m − .analyzed via the phase portrait as shown in Figure 4b. Analysis of the stick-slipmotion using phase space is detailed in [9].The temporal analysis has been historically established. However, to automat-ically scan the parametric space and benefit from modern machine learningtools, it appears interesting to analyze the dynamic behavior via single scalarvalues, used as quantities of interest for the surrogate model. Among the different scalar values proposed in the literature to investigate thedynamic system two are considered, namely the sticking time and the largestLyapunov exponent.3.1 Sticking timeThe dynamic cyclic behavior as represented in Figure 4a can be decomposedinto sticking and slipping modes. In [18] it has been proposed to investigate theoscillation behavior of dry-friction oscillator using the percent of time in whichthe body remains in stick mode as Quantity of Interest (QoI). They suggestto improve properties e.g. endurance of mechanical systems by understandingthe role of acting parameters on that QoI.In the following, the sticking time of the MoB is defined as the time duringwhich the relative velocity V R ( t ) = V − ˙ X ( t ) remains below the thresholdof 10 − m.s − . It is numerically estimated by integrating equation (9) with asix-stage, fifth-order, Runge-Kutta method with variable step-size. In order toavoid localized transient behavior, sticking time is analyzed on a time window Jan N. Fuhg, Am´elie Fau
Fig. 5 – Sticking time as indicator for chaotic motion. Upper image:Bifurcation diagram over angular frequency. Lower image: Sticking time overangular frequency.between 150 and 250 seconds. The analysis of this QoI is, however, not evidentbecause as shown in Figure 5, the value of the sticking time cannot be directlyassociated with either stable or unstable behavior. Lima and Sampaio pointedout that both the duration of the sticking mode and the cumulative stickingtime lack attention in the literature [18]. However, this analysis is not the scopeof this contribution, here the sticking time is simply considered as a candidateQoI to be represented by the metamodel.3.2 Largest Lyapunov exponentThe second QoI aims at providing an instability indicator to break down theparametric space into stable and unstable domains.Consider a general equation of motion of an N -dimensional time-continuoussystem recast into first order differential equation system of the form˙ x = f ( x ) , (10)with the state-space vector x ∈ R N and f being a set of N functions i.e. f = [ f ( x ) , . . . f N ( x )]. Consider an N -dimensional sphere of initial conditions.With time the sphere evolves into an ellipsoid whose principal axes contract orexpand with rates governed by the spectrum of LEs denoted by { λ i } i ∈ [1 ,N ] . LEsintroduced by [39] are defined as the average divergence or convergence ratesof nearby orbits in state space. A positive exponent indicates local instabilityin this particular direction and hence characterizes chaotic motion, see [40]. riging approach to investigate stick-slip instability 9 The presence of a positive LLE indicates chaos, whereas negative values arecharacteristic of regular motion. Therefore to ensure the stability of the system,computing only the LLE is sufficient. The LLE is considered as one of the mostuseful tools to characterize the stability of a dynamic system in [19]. Therefore,it should turn out as a proficient QoI for building an univariate informativesurrogate model.However, analytic calculation of the LLE is mostly restricted to simple linearsystems. For complex systems robust numerical techniques can be employedto estimate the LLE, see e.g. [41],[42], or [43].
When the investigated dynamic system is smooth with an analytically ob-tainable Jacobian matrix, the LLE can be calculated with e.g. the algorithmgiven in [44], in which the system of equations is solved for N nearby ini-tial conditions and correct state-space orientation is maintained by repeatedlyorthonormalizing the corresponding set of vectors using the Gram-Schmidtprocedure. However, the dry friction problem at hand is a non-smooth sys-tem. The Jacobian matrix cannot be determined or is strongly ill-conditioned,which leads to significant numerical problems when employing the algorithmproposed by [44]. To overcome this obstacle, a novel method has recently beenpresented in [14], which estimates the Jacobian matrix by a truncated Taylorseries expansion with small orthogonal perturbations. The spectrum of Lya-punov exponents has been estimated from that numerical method in [15] toanalyze the stability of a nonlinear mass-on-belt system. It is shortly summa-rized as follows:Consider a discretization in time given by the set { t i | i = 1 , . . . , n } . Let x ( t i )be denoted by x i . The discrete dynamical system of the mapping f introducedin equation (10) at time t i reads x i +1 = G ( x i ) , (11)where G ( x i ) = [ G ( x i ) , . . . , G N ( x i )]. Because the operator G in equation(11) may be discontinuous, an analytic expression for the Jacobian can notbe given. However, by introducing a perturbation vector ∆ i = [ δ , . . . δ N ] T ,where each element is a scalar value of small magnitude, equation (11) can berewritten in two ways x i +1 + ∆ i +1 = G ( x i + ∆ i ) ≈ G ( x i ) + J G ( x i ) ∆ i , x i +1 − ∆ i +1 = G ( x i − ∆ i ) ≈ G ( x i ) − J G ( x i ) ∆ i . (12)Here J G ( x i ) is the Jacobian matrix of G ( x i ), which can be approximated as J G ( x i ) ∆ i ≈ G ( x i + ∆ i ) − G ( x i ) ,J G ( x i ) ∆ i ≈ − G ( x i − ∆ i ) + G ( x i ) . (13) The equations correspond to forward difference and backward difference sche-mes respectively. Adding these two expressions yields an estimation of eachcolumn vector of the Jacobian as J G j ( x i ) ≈ G ( x i + ∆ ji ) − G ( x i − ∆ ji )2 δ , (14)where ∆ ji = δ e Tj and e j is the unit vector with unit in the j -th element. Fromthis robust numerical estimation of the Jacobian matrix, any algorithm forcalculating the LLE can be easily employed even for non-smooth applications.Here, the algorithm mentioned in [44] is considered. The numerical approach is verified using a continuous three-dimensional sys-tem for which the analytic Jacobian matrix can be determined. The continuouscase suggested by [45] reads ˙ x = y, ˙ y = z, ˙ z = − ax − y − z + y + xy. (15)Here a is a given parameter. The bifurcation diagram of the problem for thedomain a ∈ [3 . , .
4] is plotted in the upper image of Figure 6. It can beseen that with an increasing value of the parameter the response gets morechaotic. The estimation of LLE using the algorithm proposed in [44] based onthe exact Jacobian matrix of the system, is illustrated in the middle picture.The estimation using the same algorithm but based on the Jacobian matrixapproximated by the numerical approach previously introduced in section 3.2.1is plotted in the lowest image in Figure 6.A perturbation vector ∆ = 10 − [1 , , T has been employed. It can be ob-served that the numerical approximation of the Jacobian matrix performs well.The results provided by the analytical expression are accurately reproduced.It can also be seen that LLE mirrors the chaotic behavior of the bifurcationdiagram well.Let consider the discontinuous system of ordinary differential equations definedby equation (9). Here LLE values can only be evaluated from the numericalapproximation. Using a perturbation vector ∆ = 10 − [1 , , T , the bifurcationdiagram and the respective LLE values for this dynamic problem are plottedover the angular frequency in Figure 7. It can be seen that the positive valuesof the LLE over the parametric domain correspond well with the bifurcationgraph. Thus the LLE appears again as a pertinent instability indicator. In order to accurately determine regular or chaotic motion of the system, itis not required to quantitatively estimate the largest Lyapunov exponent, but riging approach to investigate stick-slip instability 11
Fig. 6 – Bifurcation diagram of the problem defined by Eq. (15) over theparameter a (upper image) analysed with the perspective of the LLEestimated from the analytic expression of the Jacobian matrix (middle image)and from the numerical approximation of the Jacobian matrix (lower image). Fig. 7 – Description of bifurcations by positive values of the LLE it seems rather more promising to use its non-negativity or negativity as abinary indicator.From LLE
MoB ( x ) which provides the LLE value for the MoB system withinput x over the input parametric space, an instability indicator C MoB ( x ) canbe defined as C MoB ( x ) = (cid:40) , if LLE MoB ( x ) ≥ , , if LLE MoB ( x ) < . (16)For sake of illustration, the set of LLE values represented in Figure 8a isanalyzed through this instability indicator, which is plotted in Figure 8b. Thestable behavior associated with values of the indicator C MoB = 0 is indicatedin gray, whereas the parametric domain corresponding with chaotic behaviori.e. C MoB = 1 is highlighted in red. It can be seen that this LLE surfaceprovides a valuable information which could allow to optimize the system bypreventing instabilities. (a)
LLE values over parametric space (b)
Indicator C MoB map (red: C MoB = 1,gray: C MoB = 0)
Fig. 8 – Example of LLE analysis over the spring stiffnesses ( K in N . m − , K in N . m − )Based on this proficient QoI, the dynamic system could be efficiently analyzedthrough a surrogate model to evaluate the critical values at low cost. The goalof the surrogate model is to be able to provide an accurate approximationof this surface with only few observations, originating either from physical ornumerical experiments. This surrogate model should be the best representationof the surface with respect to the available information, and provide guidancefor defining further experiments to be conducted. The idea behind surrogate models is only summarized here. Given a mapping M : X → Y between an input x ∈ X ⊂ R n and some uni-variate output riging approach to investigate stick-slip instability 13 y ∈ Y ⊂ R . From a set of m observations D = { (cid:0) x ( i ) , y ( i ) (cid:1) , i = 1 , . . . , m } ,global surrogate modeling aims to statistically approximate the output onthe whole input space, i.e. to approximate the statistical relationship behindthis mapping. The samples in the continuous input space are called designof experiments and are denoted by the set X = { x ( i ) , i = 1 , . . . , m } . Foruni-variate output the mapping evaluations are gathered in the vector y = { y i , i = 1 , . . . , m } . The statistical approximation on the whole input spacewill be given by a computed surrogate model denoted by ˜ M . Among thedifferent existing surrogate modeling approaches, kriging is investigated here.4.1 KrigingKriging was proposed by a mining engineer called Krige for geostatistics study[27]. The mathematical foundation for the approach has been developed by[46], which also established its name. For an overview the reader is referred to[22] and [21]. It can thus be outlined that the approach can be described froma regression perspective as a statistical extension of deterministic regression[47], or from a bayesian or machine learning perspective as an update fromprior information based on a conditional distribution [23]. Using ordinary kriging [21], the mapping between input and output data isapproximated by a combination of a global contribution characterized by µ , anunknown constant describing the mean contribution, and localized departuresdenoted by Z ( x ) Y ( x ) = µ + Z ( x ). (17)In details, Z ( x ) is a realization of a stochastic process with zero mean and avariance σ accounting for local deviations. The covariance matrix of Z ( x ) isgiven as Cov [ Z ( x ( i ) ) , Z ( x ( j ) )] = σ R ( x ( i ) , x ( j ) , θ ) . (18)Here R ( x ( i ) , x ( j ) , θ ) is the correlation function between two points of the inputset. θ is called hyperparameter and contains the unknown correlation param-eters between two points. It can thus be seen, as illustrated in Figure 9, thatthe ordinary kriging model depends on two main contributions, a regressionpart which is here only a mean contribution, and a stationary Gaussian pro-cess based on the assumption of an auto-correlation function, which governsthe deviation from that mean. Kriging model is based on some observations, as prior information within aBayesian approach [48]. It could also be considered as a training stage ofthe surrogate model. Considered in ordinary kriging as a constant, the prior
Fig. 9 – Scheme of an ordinary kriging metamodel Y ( x ) as the superpositionof a constant mean µ of the samples, shown by black circles, and anuncertain deviation described by the Gaussian process Z ( x ).value of the mean ˆ µ can be obtained from the available observation y using ageneralized least-square estimate asˆ µ = ( T R − ) − T R − y (19)with a vector containing only components with scalar value 1 and R thecorrelation matrix. It can be noticed that this prior evaluation depends on thechosen correlation function. The evaluation of the standard deviation from theprior information is denoted by ˆ σ . Then, the variance of the prior informationreads ˆ σ = 1 m ( y − ˆ µ ) T R − ( y − ˆ µ ) , (20)which, similarly to the mean, depends on the chosen correlation matrix. The realization of the Gaussian process depends on its characteristics, namelyits mean, its variance and its auto-correlation.The auto-correlation function, commonly referred as correlation function forsake of simplicity, is usually assumed by the user. A common choice is theMat´ern 3/2 correlation function [49], defined as R ( x − x (cid:48) , l ) = n (cid:89) i =1 (cid:32) √ | x i − x (cid:48) i | l i (cid:33) exp (cid:32) − √ | x i − x (cid:48) i | l i (cid:33) , (21)in which a different hyperparameter l i can be used for each dimension n of thedesign of experiments. Due to the lack of knowledge about the observationsthis correlation structure is arbitrary chosen here.The hyperparameters need to be determined before constructing the surro-gate model. This can for example be done by utilizing the maximum likeli-hood estimate of the hyperparameter value. This method yields an auxiliaryoptimization problem [50] for estimating θ of the formˆ θ = arg min θ (cid:63) ψ ( θ (cid:63) ) , (22) riging approach to investigate stick-slip instability 15 in which ψ is the reduced likelihood defined for ordinary kriging as ψ ( θ ) = ˆ σ ( θ )[det R ( θ )] /m , (23)with det denoting the determinant operator.Since there is no analytic solution for the optimization problem defined byEquation (23), it is necessary to use numerical optimization tools. Bouhlel andMartin mention this step as being the most challenging for the constructionof surrogate models with kriging because of the multimodality of the likeli-hood function [51]. In this work, a hybridized particle swarm optimization isemployed similar to the method suggested by [52]. Consider the unobserved value Y ≡ Y ( x (0) ) with x (0) ∈ X , which needs to bepredicted. When utilizing the best linear unbiased predictor for the Gaussianrandom variate ˆ Y , the mean of the unobserved value yields µ ˆ Y = ˆ µ + r T R − ( y − ˆ µ ) . (24) r describes the cross-correlations between the new point x (0) and each avail-able observation with r i = R ( x (0) − x ( i ) , θ ) i = 1 , . . . , m. (25)Furthermore the variance of the unobserved value can be estimated as σ Y = σ (cid:16) − r T R − r + u T (cid:0) T R − (cid:1) − u (cid:17) , (26)with u = T R − r − . (27)For proof and further details the reader shall be referred to [53].4.2 Adaptive samplingTo optimize the interest of the computations, the observation points whichdefine the experimental design are chosen with respect to an adaptive samplingscheme. Starting from a set of initial data D ini = { (cid:0) x ( i ) , y ( i ) (cid:1) , i = 1 , . . . , m } a surrogate model ˜ M is constructed with kriging based on this knowledge.By solving some auxiliary optimization problems new points are successivelyadded to the dataset D until a convergence criterion is reached. The generaloptimization problem reads x ( m +1) = arg min x (cid:63) ∈ X Score ( x (cid:63) ) . (28)Here a score function is evaluated which describes generally a tradeoff betweenan exploration and an exploitation component. Using the exploration term, the domain is globally examined in order to detect unknown regions of interest.The exploitation aims to generate data points in the pre-identified regions ofinterest to reduce the prediction error locally. Different techniques have beenproposed in the literature such as [54, 55, 56].An external variance-based technique has recently been proposed in [57]. Max-imizing Expected Prediction Error (MEPE) is a continuous optimization in-volving a switch strategy between exploitation and exploration. The bias erroris estimated by a continuously approximated leave-one-out cross validation(LOOCV) error ˆ e ( x ). The continuous expected prediction error to be maxi-mized is given by EP E ( x ) = α ˆ e ( x ) + (1 − α )ˆ σ Y , (29)where ˆ e ( x ) is the value of ˆ e ( x ( i ) ) when x is located in the Voronoi cell ofpoint x ( i ) . Otherwise its value is zero. In order to fasten the computation anapproximation of ˆ e ( x ( i ) ) can be used [57, 58]. A balance factor α is intro-duced to adjust the exploitative bias term and the exploratory variance termby switching adaptively between the two components depending on the esti-mation quality of the bias term. The computation of the balance factor for thecalculation of the sample point x ( m + q − reads α = . , if q = 10 .
99 min (cid:20) . e true ( x ( m + q − )ˆ e ( x ( m + q − ) , (cid:21) if q > . (30)Here, e true represents the error between the surrogate model and the actualobservation for this sample point. The algorithm is summarized in Box 1. – Given a design of experiments X = { x (1) , . . . , x ( m ) } .While the adaptive sampling stopping criterion is not satisfied.Do: Update the balance factor α with equation (30).Obtain the LOOCV error at each sample point ˆ e ( x ( i ) ).Obtain a new point by maximizing the EP E -criterion of equation (29) overthe inputUpdate the design of experiments with the new found point and its response.end
Box 1 – Algorithm for MEPEThe second scheme considered here is Expected Improvement for Global Fit(EIGF) [59], which is also a continuous optimization scheme, but with a fixedbalance between exploitation and exploration. It is based on the maximizationof the expected improvement defined as EI ( x ) = ( µ ˆ Y ( x ) − y ( x ∗ )) + σ Y ( x ) (31) riging approach to investigate stick-slip instability 17 where x ∗ is the closest sample point to the candidate x . The first term, whichrepresents the exploitation contribution, is large when the surrogate modelˆ Y ( x ) differs substantially from the response at the nearest point. The secondterm provides global exploration, it is large when the surrogate model exhibitshigher lack of knowledge.For classification problems, a dedicated adaptive scheme called Monte Carlo-Intersite Voronoi (MiVor) has recently been proposed [60, 61]. That scheme isbased on the exploration of the parametric space by Voronoi tessellation. Ascore is given to all cells to identify the one cell corresponding to the largest lackof knowledge near to a decision region between two classes and so to pertinentlyadd a new observation point. A random switching strategy is used to balanceexploration and exploitation during the successive iterations. This allows toreach a compromise between the exploration of the whole parametric space andthe addition in some local areas near to the class boundaries which require tobe accurately described. The scheme has been described in detail, investigatedfor various applications and compared with other adaptive schemes in [60]. The surrogate strategy is tested on different applications. The problem param-eters are summarized in Table 1. One-, two- and three-dimensional parametricspaces are considered.
Problem P P P P P P M V . . s − D . . s . m − µ s . V s . . s − U . N . σ . . m − σ . . s . m − σ . . s . m − Ω (in rad . s − ) 0 . . , .
0] 0 . . . . , . K (in N . m − ) [0 . , .
0] 1 . . , .
0] 1 . . , .
0] 1 . K (in N . m − ) [0 . , .
6] 0 . . , .
5] [0 . , .
5] [0 . , .
6] [0 . , . µ k .
15 0 .
15 0 .
15 [0 . , .
18] 0.15 [0 . , . Table 1 – Problem parametersThe performances for both QoIs are analyzed with respect to a reference so-lution computed with m ref evaluation points { x ref,i , i = 1 , . . . , m ref } . 5000points per dimension are considered for the applications. Four indicators areconsidered to evaluate the error between the metamodel ˆ Y and the referenceresponse surface y ref . The Mean Absolute Error (MAE) is estimated consid- ering the reference points asMAE = 1 m ref m ref (cid:88) i =1 (cid:12)(cid:12)(cid:12) y ref ( x ref,i ) − ˆ Y ( x ref,i ) (cid:12)(cid:12)(cid:12) . (32)With the same perspective, the Root Mean-Squared Error (RMSE) is esti-mated as RMSE = (cid:118)(cid:117)(cid:117)(cid:116) m ref m ref (cid:88) i =1 (cid:16) y ref ( x ref,i ) − ˆ Y ( x ref,i ) (cid:17) . (33)The Relative Maximum Absolute Error (RMAE) is given byRMAE = max i ∈ [1 ,m ref ] (cid:16)(cid:12)(cid:12)(cid:12) y ref ( x ref,i ) − ˆ Y ( x ref,i ) (cid:12)(cid:12)(cid:12)(cid:17) σ y ref , (34)with σ y ref the standard deviation of the reference solution. Finally the R score reads R = 1 − (cid:80) m ref i =1 (cid:16) y ref ( x ref,i ) − ˆ Y ( x ref,i ) (cid:17) (cid:80) m ref i =1 (cid:0) y ref ( x ref,i ) − y ref (cid:1) , (35)where y ref denotes the mean of the reference response values. The parametricinput space is normalized in each dimension in order to avoid numerical issues.Let x ui and x li be the upper and lower limit of the dimension i respectively.The normalized input x i is then given by¯ x i = x i − x li x ui − x li . (36)The normalization ensures that 0 ≤ ¯ x i ≤ P First, the ability of the surrogate model to provide an estimation for the stick-ing time is investigated. The parametric domain of the problem P is given bythe two spring stiffnesses with K ∈ [0 . , .
0] N . m − and K ∈ [0 . , .
6] N . m − .Other parameter values are considered to be deterministic, see Table 1.The reference response surface for the sticking time over the input domain isplotted in Figure 10a. Initially 20 samples are spread across the domain withTranslational Propagation Latin Hypercube Design (TPLHD) [62] leading toa metamodel with an absolute error distributed as shown in Figure 10b. Itcan be noticed that the spread of the sticking time in the parametric domain { K , K } is around 5 seconds and the highest initial absolute error is around1.6 seconds. riging approach to investigate stick-slip instability 19 (a) Reference solution (b)
Contour of absolute error (in s) forthe initial surrogate model (20 sampleswith TPLHD)
Fig. 10 – Sticking time full response surface (Problem P ) and initialabsolute error of the metamodel. K is given in N . m − , K in N . m − . From the initial TPLHD model with 20 samples, surrogate models utilizing60 samples are built considering the two alternative sampling techniques pre-viously introduced, i.e. Maximizing Expected Prediction Error (MEPE) andExpected Improvement for Global Fit (EIGF). Furthermore a direct TPLHDmodel with 60 samples is generated. In order to avoid numerical outliers of thestochastic optimization process each strategy is executed twenty times, so thatthe performances can be examined as mean performance of these realizations.The error measures for the different surrogate models are listed in Table 2.The two adaptive schemes provide much more proficient metamodels than theinitial and final one-shot metamodels based on TPLHD. Among them EIGFappears more proficient than MEPE.
Method MAE [s] RMAE RMSE [s] R Metamodelwith 20 samples TPLHD 0.31 1.29 0.46 0.88Metamodelwith 60 samples TPLHD 0.28 1.27 0.38 0.93EIGF 0.16 1.00 0.22 0.97MEPE 0.17 1.09 0.27 0.96
Table 2 – Error measures for the sticking time metamodels (Problem P )with 60 samples. In Figure 11 the convergence of MAE error and RMSE error with the numberof samples is plotted for the two adaptive methods. As before, all resultsare averaged on twenty realizations. It seems that MEPE performs better asshown in Figure 11a. However, by considering the convergence in terms ofRMSE error, both methods show a similar convergence behaviour, see Figure11b. (a)
MAE error in s (b)
RMSE error in s
Fig. 11 – Evolution of the errors with the sample size for the two adaptivesampling techniques for the sticking time metamodels (Problem P ) until150 samples.A set of sample positions of the added points of the two adaptive methods after150 samples are shown in Figure 12. On these plots the parameters have beennormalized to lie between 0 and 1. The initial TPLHD sample points are shownas black circles. The later a point is added the more its color tends towardslight red. The contour of the target function is displayed for both cases, it maybe observed in Figure 12b that MEPE balances proficiently exploration andexploitation, wheras EIGF focuses on specific domains. (a) EIGF - 150 samples (b)
MEPE - 150 samples
Fig. 12 – Sample positions after 150 samples for the sticking time problem(Problem P ). riging approach to investigate stick-slip instability 21 In this subsection ordinary kriging is compared to three commonly found sur-rogate modeling techniques, i.e. support-vector machines, radial basis functionnetwork and neural networks. Here, the metamodels are here considered with-out adaptive scheme, and are all based on the same sample set provided byTPLHD. For this comparison the support-vector machine is employed with aradial basis kernel as well as an automatic choice of the scale value for thekernel function as defined in the MATLAB software, see e.g. [24]. For the ra-dial basis function network [25], the hidden layer has 50 neurons. The centersare calculated with k-means clustering algorithm and the widths are set tounity. The neural network is employed in form of multilayer perceptron basedon 2 hidden-layer with 20 neurons each and sigmoid activation function whereLevenberg-Marquardt backpropagation is used for training.The mean absolute errors as defined by Equation (32) for the obtained meta-models considering different sample sizes are listed in Table 3. It can be noticedthat the ranking of the methods depends on the number of samples considered.Support-vector machine performs better than radial basis function networkand neural network for a low number of samples whereas it performs worse fornumber of samples larger than 100. Considering the parameters as detailed be-fore, the results indicate that ordinary kriging yields the best approximationof the target function for all sample sizes for the problem of interest, whilesimultaneously provding an estimate for the variance of the metamodel overthe parametric space.
Number of samplesMetamodel method 20 50 100 150 200
Ordinary kriging 0.31 0.18 0.15 0.10 0.08
Support vector machine 0.37 0.32 0.28 0.26 0.25Radial basis function network 0.46 0.46 0.42 0.25 0.23Neural network 0.84 0.47 0.32 0.23 0.16
Table 3 – Mean absolute errors in s for the two-dimensional sticking timeproblem (Problem P ) with four alternative metamodeling techniques anddifferent sample sizes.5.2 Surrogate model for classification of chaotic motion using LLEThe second focus of interest is to predict the parametric subdomains corre-sponding with chaotic motion of the oscillator at low computational cost usingthe kriging approach. It has been shown in [61, 60] that issues arise when try-ing to utilize the common adaptive sampling techniques for a regressive meta-model for that application. Indeed, as previously discussed, the precise value of the LLE is not of particular interest. The main focus of a metamodel forthis application is to classify ˆ M LLE,MoB ( x ) the approximated estimation of M LLE,MoB ( x ) for a point x of the parametric domain either above or below 0to characterize the behavior of the given dynamic system as chaotic or regularmotion. Using the adaptive sampling technique MiVor recently proposed in[60], specifically dedicated to classification based on regression metamodel, ishereafter investigated to classify the LLE values. The required MiVor param-eters (see [60]), which are the initial exploration rate and the decrease factor,are chosen to be 0 . . P The first case of interest is the one-dimensional problem P , where the goalof the metamodel is to identify the unstable behavior corresponding with theLLE values plotted in Figure 13a. The angular frequency Ω varies between0 . . . s − , whereas the other input parameters have fixed values asdefined in Table 1. (a) LLE value and metamodel forclassification (b)
Evolution of MiVor performance withnumber of samples
Fig. 13 – One-dimensional LLE example (Problem P ). Ω is given in rad . s − The performance of the metamodel is evaluated with respect to a referencesolution obtained with 5000 sample points. The metamodel for classificationis judged by how many of the points yielding C MoB = 0 on one hand and howmany points yielding C MoB = 1 on the other hand are predicted. The measureis stated in percent of accurately predicted points. The percentage of correctlypredicted points in the unstable regime is denoted by a p LLE ≥ and defined as a p LLE ≥ = ˆ n LLE ≥ ref, LLE ≥ n ref, LLE ≥ (37)with ˆ n LLE ≥ ref, LLE ≥ the number of points among the ones predicted with unstableregime by the reference solution which are also predicted with unstable regime riging approach to investigate stick-slip instability 23 using the metamodel, and n ref, LLE ≥ the number of points estimated in theunstable regime by the reference response surface. Furthermore, the percentageof correctly predicted points in the stable regime a p LLE < is given by a p LLE < = ˆ n LLE < ref, LLE < n ref, LLE < , (38)where n ref, LLE < is the number of points with regular behavior given by thereference response surface. Besides ˆ n LLE < ref, LLE < denotes the number of pointsestimated with regular behavior using the kriging approach among the pointswhich are predicted with regular behavior using the reference response surface.The metamodel is built from the prior information given by five initial samplepoints set with TPLHD, i.e. [0 . , . , . , . , . T . None of them yields achaotic output. The challenged posed to the adaptive metamodeling strategyis to evaluate if the scheme is able to explore the parametric domain to identifythe chaotic subdomain.The reference solution, the metamodel ˆ M after 35 samples and the respectivesample positions at the end of the adaptive process are shown in Figure 13a.It can be seen that a proficient metamodel has been created with only fewobservation points. Furthermore as intended with MiVor most of the createdsamples lie in the area that indicates chaotic behavior.The evolution of the accuracy with regard to the sample size is shown in Figure13b. As an arbitrarily chosen stopping criterion the procedure stops after 35samples, i.e. after adding 30 observation points to the initial samples. Finally,the metamodel correctly evaluates 99 .
11 % of the 5000 random points thatyield LLE ≥ .
23 % of the points with LLE < P Next, the metamodel approach is tested for the test problem P which aimsat evaluating the ability of the proposed approach to proficiently explore theparametric domain if only a small percentage of the domain has a chaoticmotion.Consider the two-dimensional input domain for the two spring stiffnesses givenby K ∈ [0 . , .
0] N . m − and K ∈ [0 . , .
5] N . m − , see Table 1 for the valuesof other parameters. The plot of the reference response surface M LLE,MoB provided by 10000 computations which were spread over the input domain ina space-filling manner using TPLHD is shown in Figure 14a. Large fluctuationson the LLE values can be observed. The analysis of the parametric domainthrough binary classification corresponding with either stable or chaotic mo-tion is plotted in Figure 14b. As before, regular motion is represented in greywhereas the red area informs about chaotic behavior. It can be seen that onlya small fraction of the domain yields chaotic motion. Among the 10000 pointscomputed to build the reference solution only 93 points lead to stick-slip in-stability, i.e. only less than 1 % of the points.The metamodel is built from an initial set of 10 initial experiments createdwith TPLHD, then adaptive sampling is used until reaching 65 samples. The (a) M LLE,MoB (b)
Classification provided by M LLE,MoB (c)
65 MiVor samples (d)
Classification provided by ˆ M LLE,MoB with 60 samples
Fig. 14 – LLE response surface, LLE classes and metamodel approximationand observations for the 2D LLE case P . ¯ K and ¯ K are the normalizedspring stiffnesses.point locations of the 55 samples created with MiVor over the normalizedinput space are highlighted in Figure 14c. It can be observed that they arepredominantly spread in and around the red area that indicates chaos. Thisleads to an accurate surrogate model ˆ M LLE,MoB for the classification problemas shown in Figure 14d. In details, it yields that 98 .
93 % of the points aboveor equal 0 are correctly classified and 99 .
91 % of the points below 0.The convergence of the error measures for MiVor is depicted in Figure 15. Itcan be seen that a pLLE ≥ MiVor shows a big jump around 15 sample points.This is due to the fact that the algorithm needs to identify the small area ofFigure 14b. After the jump the metric quickly converges to the optimal state.As a large part of the LLE domain is below the threshold value a pLLE< appearsalmost accurate from the beginning with values varying between 99 . −
100 %. P Problem P aims at evaluating metamodel performance to detect almost dis-connected regions of the input domain associated with chaotic behavior.Consider the input domain for the spring stiffness given by K ∈ [0 . , .
5] N . m − and the kinematic friction coefficient µ k ∈ [0 . , . riging approach to investigate stick-slip instability 25 Fig. 15 – Convergence of correctly classified points for 2D LLE case P .fluctuations. This reference solution yields the classification displayed in Figure16b. It can be seen that there are two chaotic subdomains which are almostdisconnected. (a) M LLE,MoB (b)
Classification provided by M LLE,MoB (c)
105 MiVor samples (d)
Classification provided by ˆ M LLE,MoB with 105 samples
Fig. 16 – LLE response surface, LLE classes, metamodel approximation andobservations for the 2D LLE case P . The metamodel is constructed adaptively from 5 points sampled with TPLHD.It is evaluated after adding 100 sample points and the results are averaged over20 independent realizations. The locations of the 105 MiVor sample points forone realization are shown in Figure 16c. It can be noticed that the methodis effective in sampling around the red zones indicating chaos. Furthermorethe rest of the domain is sampled with a space-filling design. The resultingmetamodel after 105 samples is displayed in Figure 16d and shows an accurateprediction performance. The percentage of correctly identified points for LLE ≥ . .
19% of the points corresponding to LLE values below 0are correctly identified. The convergence of the percentage values is displayedin Figure 17. It shows a very proficient sampling performance of MiVor.
Fig. 17 – Convergence of correctly classified points for 2D LLE problem P . P The metamodel approach for identifying stick-slip instability is tested on Prob-lem P , characterized by a highly complex response surface, in which regularand chaotic motion subdomains are embedded into each other. Consider theparametric domain for the two spring stiffnesses given by K ∈ [0 . , .
0] N . m − and K ∈ [0 . , .
6] N . m − with other parameters as given in Table 1.The reference LLE response surface based on the evaluation of 10000 samplepoints is shown in Figure 18a. The corresponding classification is picturedin Figure 18b. It can be noticed that the chaotic area is complex with someregular subdomains embedded within chaotic behavior.A set of 5 points sampled using the TPLHD technique is considered as initialobservation points. One realization is represented in Figures 18c and 18d. Theadaptive sampling technique is evaluated after reaching 115 sample points.Samples needed to generate the metamodel are shown in Figure 18c. Theresulting metamodel is represented in Figure 18d. It can be seen that the holes riging approach to investigate stick-slip instability 27 inside the red domain are not detected by the surrogate model. This reducesthe percentage of correctly classified points. However an accurate prediction ofthese holes with kriging would require a significantly larger number of samples.Furthermore identifying a point inside the regular motion subdomain mightlead to a perceived reliability of the system dynamics around the observationpoint which might not be desired from an engineering point of view. (a) M LLE,MoB (b)
Classification provided by M LLE,MoB (c)
115 MiVor samples (d)
Classification provided by ˆ M LLE,MoB with 115 samples
Fig. 18 – LLE response surface, LLE classes, metamodel approximation andobservations for the 2D LLE case P . ¯ K and ¯ K are the normalized springstiffnesses. LLE is unitless.The convergence of the error metrics are displayed in Figure 19, which shows anefficient convergence behavior. The results of the MiVor metamodel after 115samples are very proficient with values of a pLLE ≥ = 97 .
12 % and a pLLE< =97 .
32 %. Considering the very complex form of the chaotic region includingnotched boundary and holes as shown in Figure 18b, these percentages ofaccurately identified points computing only 115 observation points are highlypromising for detecting non-stable dynamic parametric domains at low cost.
Fig. 19 – Convergence of correctly classified points for 2D LLE problem P .5.4 Problem P Finally the input dimension is extended to the three-dimensional problem P .Consider the parametric space Ω × K × µ K ∈ [0 . , .
9] rad . s − × [0 . , .
5] N . m − × [0 . , .
15] with K = 1 N . m − . The rest of the parameters remain unchangedwith regard to the previous examples as given in Table 1. As reference 15000observations, which are evenly spread in the parametric space, are evaluated.The MiVor process is started from an initial set of 5 TPLHD samples, noneof which correspond to an unstable dynamic system. The evolution of theaveraged error values until including 200 samples is displayed in Figure 20.
50 100 150 200020406080100
Fig. 20 – Convergence of correctly classified points for 3D LLE problem P .It can be noticed that MiVor achieves proficient value convergence with bothmetrics indicating an exactitude of above 98% of accurately identified pointsfor both classes with already around 125 observation points. The jittering of riging approach to investigate stick-slip instability 29 the values indicates sharp drops of the LLE response surface. It can be seenthat MiVor can be very well expanded to higher dimension on that specificexample. Using the latest development proposed in the literature both in terms of nu-merical approaches allowing to compute the LLE and kriging models for clas-sification, a proficient analysis of the parametric model can be conducted tooptimally improve the knowledge of dynamic instabilities. A non-linear oscilla-tor of Duffing’s type with an elasto-plastic friction model has been investigated.The largest Lyapunov exponent as well as the sticking time have been intro-duced as means to investigate the system behavior. The analysis of this systemwhich appeared highly challenging, as both QoIs have highly fluctuating val-ues in the parametric domain, is very efficiently tackled using an adaptivemetamodeling approach. The LLE appears as the most promising QoI in thatcontext. Values equal or above zero indicate chaotic motion accurately. TheMiVor sampling technique is able to generate surrogate models able to classifychaotic and regulate motion with only a few observations. It efficiently gainsinformation from the available observation to further guide the investigationand thus allows to identify the instability domain(s) with reasonable numberof computations. It therefore provides pertinent information for design opti-mization, even for complex response surfaces and two- or three-dimensionalparametric spaces.
Acknowledgements
The authors acknowledge the financial support from the DeutscheForschungsgemeinschaft under Germanys Excellence Strategy within the Cluster of Excel-lence PhoenixD (EXC 2122, Project ID 390833453).The results presented in this paper were partially carried out on the cluster system at theLeibniz University of Hannover, Germany.
Conflict of interest
The authors declare that they have no conflict of interest.
References
1. E Rabinowicz. Stick and slip.
Scientific American , 194(5):109–119, 1956.2. D Barton and A Blackwood.
Braking 2004: Vehicle Braking and ChassisControl , volume 6. John Wiley & Sons, 2004.3. N Ashraf, D Bryant, and JD Fieldhouse. Investigation of stick-slip vi-bration in a commercial vehicle brake assembly.
International Journal ofAcoustics & Vibration , 22(3), 2017.4. William S Owen and Elizabeth A Croft. The reduction of stick-slip fric-tion in hydraulic actuators.
IEEE/ASME transactions on mechatronics ,8(3):362–371, 2003.5. Q Wu, S Luo, Qu T, and X Yang. Comparisons of draft gear dampingmechanisms.
Vehicle System Dynamics , 55(4):501–516, 2017.
6. D Rubio and L San Andres. Structural stiffness, dry friction coefficient,and equivalent viscous damping in a bump-type foil gas bearing.
Journalof engineering for gas turbines and power , 129(2):494–502, 2007.7. P Stelter. Stick-slip vibrations and chaos.
Phil. Trans. R. Soc. Lond. A ,332(1624):89–105, 1990.8. N Hinrichs, M Oestreich, and K Popp. Dynamics of oscillators with impactand friction.
Chaos, Solitons & Fractals , 8(4):535–558, 1997.9. U Galvanetto and S Bishop. Dynamics of a simple damped oscillatorundergoing stick-slip vibrations.
Meccanica , 34(5):337–347, 1999.10. M Jim´enez, J Bielsa, R Rodr´ıguez, and C Bernad. Two FEM approachesfor the prediction and quantification of stick-slip phenomena on rubber-metal sliding contacts. In
IUTAM Symposium on Computational Methodsin Contact Mechanics , pages 291–309. Springer, 2007.11. K Devarajan and B Balaram. Analytical approximations for stick–slipamplitudes and frequency of duffing oscillator.
Journal of Computationaland Nonlinear Dynamics , 12(4):044501, 2017.12. B Santhosh, S Narayanan, and C Padmanabhan. Discontinuity inducedbifurcations in nonlinear systems.
Procedia IUTAM , 19:219–227, 2016.13. J Awrejcewicz and D Sendkowski. Stick-slip chaos detection in coupledoscillators with friction.
International Journal of Solids and Structures ,42(21-22):5669–5682, 2005.14. M Balcerzak, A Dabrowski, A Stefa´nski, and J Wojewoda. Spectrum ofLyapunov exponents in non-smooth systems evaluated using orthogonalperturbation vectors. In
MATEC Web of Conferences , volume 148, page10003. EDP Sciences, 2018.15. D Pikunov and A Stefanski. Numerical analysis of the friction-inducedoscillator of Duffing’s type with modified LuGre friction model.
Journalof Sound and Vibration , 440:23–33, 2019.16. G Benetti, L Galgani, and J Strelcyn. Kolmogorov entropy and numericalexperiments.
Physical Review A , 14(6):2338, 1976.17. P Grassberger and I Procaccia. Characterization of strange attractors.
Physical Review Letters , 50(5):346, 1983.18. R Lima and R Sampaio. Stick-mode duration of a dry-friction oscillatorwith an uncertain model.
Journal of Sound and Vibration , 353:259–271,2015.19. L Kocarev, J Szczepanski, J Amig´o, and I Tomovski. Discrete chaos-i:Theory.
IEEE Transactions on Circuits and Systems I: Regular Papers ,53(6):1300–1309, 2006.20. J Awrejcewicz, A Krysko, N Erofeev, V Dobriyan, M Barulina, andV Krysko. Quantifying chaos by various computational methods. part1: Simple systems.
Entropy , 20(3):175, 2018.21. J Kleijnen. Regression and kriging metamodels with their experimentaldesigns in simulation: a review.
European Journal of Operational Research ,256(1):1–16, 2017.22. J Kleijnen. Kriging metamodeling in simulation: A review.
Europeanjournal of operational research , 192(3):707–716, 2009. riging approach to investigate stick-slip instability 31
23. C Williams and C Rasmussen. Gaussian processes for regression. In
Ad-vances in neural information processing systems , pages 514–520, 1996.24. S Clarke, J Griebsch, and T Simpson. Analysis of support vector regressionfor approximation of complex engineering analyses.
Journal of mechanicaldesign , 127(6):1077–1087, 2005.25. J Park and I Sandberg. Universal approximation using radial-basis-function networks.
Neural computation , 3(2):246–257, 1991.26. P Jiang, Y Zhang, Q Zhou, X Shao, J Hu, and L Shu. An adaptivesampling strategy for kriging metamodel based on Delaunay triangulationand topsis.
Applied Intelligence , pages 1–13, 2017.27. D Krige. A statistical approach to some basic mine valuation problems onthe witwatersrand.
Journal of the Southern African Institute of Miningand Metallurgy , 52(6):119–139, 1951.28. J Sacks, W Welch, T Mitchell, and H Wynn. Design and analysis ofcomputer experiments.
Statistical science , pages 409–423, 1989.29. W Van Beers and J Kleijnen. Kriging for interpolation in random simu-lation.
Journal of the Operational Research Society , 54(3):255–262, 2003.30. E Pennestr`ı, V Rossi, P Salvini, and P Valentini. Review and comparisonof dry friction force models.
Nonlinear dynamics , 83(4):1785–1801, 2016.31. H Olsson, K ˚Astr¨om, C De Wit, M G¨afvert, and P Lischinsky. Frictionmodels and friction compensation.
Eur. J. Control , 4(3):176–195, 1998.32. F Marques, P Flores, JC Claro, and H Lankarani. A survey and com-parison of several friction force models for dynamic analysis of multibodymechanical systems.
Nonlinear Dynamics , 86(3):1407–1443, 2016.33. P Dupont, V Hayward, B Armstrong, and F Altpeter. Single stateelastoplastic friction models.
IEEE Transactions on automatic control ,47(5):787–792, 2002.34. L Prandtl. Spannungsverteilung in plastischen K¨orpern. In
Proceedings ofthe 1st International Congress on Applied Mechanics , pages 43–54, 1924.35. P Dupont, B Armstrong, and V Hayward. Elasto-plastic friction model:contact compliance and stiction. In
American Control Conference, 2000.Proceedings of the 2000 , volume 2, pages 1072–1077. IEEE, 2000.36. C De Wit, H Olsson, K Astrom, and P Lischinsky. A new model forcontrol of systems with friction.
IEEE Transactions on automatic control ,40(3):419–425, 1995.37. W Townsend and Jr Salisbury. The effect of Coulomb friction and stictionon force control. In
Proceedings. 1987 IEEE International Conference onRobotics and Automation , volume 4, pages 883–889. IEEE, 1987.38. R Stribeck. Die wesentlichen Eigenschaften der Gleit-und Rollenlager.
Zeitschrift des Vereines Deutscher Ingenieure , 46:1341–1348, 1902.39. V Oseledec. A multiplicative ergodic theorem. Liapunov characteristicnumber for dynamical systems.
Trans. Moscow Math. Soc. , 19:197–231,1968.40. M Rosenstein, J Collins, and C De Luca. A practical method for cal-culating largest Lyapunov exponents from small data sets.
Physica D:Nonlinear Phenomena , 65(1-2):117–134, 1993.
41. I Shimada and T Nagashima. A numerical approach to ergodic problem ofdissipative dynamical systems.
Progress of theoretical physics , 61(6):1605–1616, 1979.42. G Benettin, L Galgani, A Giorgilli, and J Strelcyn. Lyapunov characteris-tic exponents for smooth dynamical systems and for Hamiltonian systems;a method for computing all of them. part 1: Theory.
Meccanica , 15(1):9–20, 1980.43. H Kantz. A robust method to estimate the maximal Lyapunov exponentof a time series.
Physics letters A , 185(1):77–87, 1994.44. A Wolf. Quantifying chaos with Lyapunov exponents.
Chaos , 16:285–317,1986.45. M Molaie, S Jafari, J Sprott, and S Golpayegani. Simple chaotic flows withone stable equilibrium.
International Journal of Bifurcation and Chaos ,23(11):1350188, 2013.46. G Matheron. Principles of geostatistics.
Economic geology , 58(8):1246–1266, 1963.47. A Stein and LCA Corsten. Universal kriging and cokriging as a regressionprocedure.
Biometrics , pages 575–587, 1991.48. M Handcock and M Stein. A Bayesian analysis of kriging.
Technometrics ,35(4):403–410, 1993.49. B Mat´ern. Spatial variation: Meddelanden fran statens skogsforskningsin-stitut.
Lecture Notes in Statistics , 36:21, 1960.50. V Dubourg, B Sudret, and F Deheeger. Metamodel-based importancesampling for structural reliability analysis.
Probabilistic Engineering Me-chanics , 33:47–57, 2013.51. M Bouhlel and J Martins. Gradient-enhanced kriging for high-dimensionalproblems.
Engineering with Computers , 35(1):157–173, 2019.52. D Toal, N Bressloff, A Keane, and C Holden. The development of ahybridized particle swarm for kriging hyperparameter tuning.
Engineeringoptimization , 43(6):675–699, 2011.53. T Santner, B Williams, and W Notz.
The design and analysis of computerexperiments . Springer Science & Business Media, 2013.54. D Jones, M Schonlau, and W Welch. Efficient global optimization ofexpensive black-box functions.
Journal of Global optimization , 13(4):455–492, 1998.55. Cameron J Turner, Richard H Crawford, and Matthew I Campbell. Mul-tidimensional sequential sampling for NURBs-based metamodel develop-ment.
Engineering with Computers , 23(3):155–174, 2007.56. P Singh, D Deschrijver, and T Dhaene. A balanced sequential designstrategy for global surrogate modeling. In
Simulation Conference (WSC),2013 Winter , pages 2172–2179. IEEE, 2013.57. H Liu, J Cai, and Y Ong. An adaptive sampling approach for krigingmetamodeling by maximizing expected prediction error.
Computers &Chemical Engineering , 106:171–182, 2017.58. S Sundararajan and S Keerthi. Predictive approaches for choosing hy-perparameters in Gaussian processes. In
Advances in neural information riging approach to investigate stick-slip instability 33 processing systems , pages 631–637, 2000.59. C Lam.
Sequential adaptive designs in computer experiments for responsesurface model fit . PhD thesis, The Ohio State University, 2008.60. JN Fuhg and A Fau. An innovative adaptive kriging approach for effi-cient binary classification of mechanical problems. manuscript submittedto publication , in 2019.61. JN Fuhg. Adaptive surrogate models for parametric studies.Master’s thesis, Leibniz Universit¨at Hannover, Arxiv platformhttps://arxiv.org/abs/1905.05345, 2019.62. F Viana, G Venter, and V Balabanov. An algorithm for fast optimallatin hypercube design of experiments.