Surrogate model approach for investigating the stability of a friction-induced oscillator of Duffing's type

TL;DR

This paper introduces a surrogate modeling approach combining Gaussian process regression and adaptive classification to efficiently explore the stability of a Duffing-type oscillator with elasto-plastic friction, enabling rapid detection of instability regions.

Contribution

It presents a novel surrogate strategy integrating kriging and MiVor classification for stability analysis of nonlinear oscillators with complex friction models.

Findings

Efficient detection of instability domains in parametric space.

Use of Lyapunov exponent as an indicator of non-regular behavior.

High proficiency of the MiVor classification in complex response surfaces.

Abstract

Parametric studies for dynamic systems are of high interest to detect instability domains. This prediction can be demanding as it requires a refined exploration of the parametric space due to the disrupted mechanical behavior. In this paper, an efficient surrogate strategy is proposed to investigate the behavior of an oscillator of Duffing's type in combination with an elasto-plastic friction force model. Relevant quantities of interest are discussed. Sticking time is considered using a machine learning technique based on Gaussian processes called kriging. The largest Lyapunov exponent is proposed as an efficient indicator of non-regular behavior. This indicator is estimated using a perturbation method. A dedicated adaptive kriging strategy for classification called MiVor is utilized and appears to be highly proficient in order to detect instabilities over the parametric space and can furthermore be used for complex response surfaces in multi-dimensional parametric domains.

Tables

Table 1. Table 1: Problem parameters

	Problem
	${𝒫}_0$	${𝒫}_1$	${𝒫}_2$	${𝒫}_3$	${𝒫}_4$	${𝒫}_5$
$M$	$1\text{kg}$
$V_0$	$0.1\text{m}.{\text{s}}^{-1}$
$D$	$0.0\text{N}.\text{s}.{\text{m}}^{-1}$
${\mu }_s$	$0.3$
$V_s$	$0.1\text{m}.{\text{s}}^{-1}$
$U_0$	$0.1\text{N}$
$N_0$	$1.0\text{N}$
${\sigma }_0$	$100.0\text{N}.{\text{m}}^{-1}$
${\sigma }_1$	$10.0\text{N}.\text{s}.{\text{m}}^{-1}$
${\sigma }_2$	$0.1\text{N}.\text{s}.{\text{m}}^{-1}$
$\Omega$ (in $\text{rad}.{\text{s}}^{-1}$)	$0.6$	$[0.2,1.0]$	$0.6$	$0.6$	$0.7$	$[0.6,0.9]$
$K_1$ (in $\text{N}.{\text{m}}^{-3}$)	$[0.5,1.0]$	$1.0$	$[0.5,1.0]$	$1.0$	$[0.5,1.0]$	$1.0$
$K_2$ (in $\text{N}.{\text{m}}^{-1}$)	$[0.0,0.6]$	$0.0$	$[0.0,0.5]$	$[0.0,0.5]$	$[0.0,0.6]$	$[0.0,0.5]$
${\mu }_k$	$0.15$	$0.15$	$0.15$	$[0.08,0.18]$	0.15	$[0.10,0.15]$

Table 2. Table 2: Error measures for the sticking time metamodels (Problem 𝒫 0 subscript 𝒫 0 \mathcal{P}_{0} ) with 60 samples.

	Method	MAE [s]	RMAE	RMSE [s]	R2
Metamodel with 20 samples	TPLHD	0.31	1.29	0.46	0.88
Metamodel with 60 samples	TPLHD	0.28	1.27	0.38	0.93
	EIGF	0.16	1.00	0.22	0.97
	MEPE	0.17	1.09	0.27	0.96

Table 3. Table 3: Mean absolute errors in s for the two-dimensional sticking time problem (Problem 𝒫 0 subscript 𝒫 0 \mathcal{P}_{0} ) with four alternative metamodeling techniques and different sample sizes.

	Number of samples
Metamodel 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.25
Radial basis function network	0.46	0.46	0.42	0.25	0.23
Neural network	0.84	0.47	0.32	0.23	0.16

Equations

$$M\ddot{X}(t)=-D\dot{X}(t)-K_1{X}^3(t)-K_{2}X(t)+F_R(V_R)+U_0\sin (\Omega t)\text{,}$$

$$X=z+w\text{.}$$

$$f_R(V_R,z)={\sigma }_{0}z+{\sigma }_1\dot{z}+{\sigma }_2{V}_R,$$

$$\text{sgn}\left(V_R\right)=\{\begin{array}{ll}\frac{V_R}{\text{norm}{V}_R},& \text{if}{V}_R\ne 0,\\ 0,& \text{if}{V}_R=0,\end{array}$$

$$\dot{z}=\left(1-\alpha (z,V_R)\frac{{\sigma }_0}{g(V_R)}z\cdot \text{sgn}\left(V_R\right)\right)V_R\text{,}$$

$$\alpha (z,V_R)=\{\begin{array}{ll}\alpha (z),& \text{if}{V}_R\cdot z\ge 0\text{,}\\ 0,& \text{if}{V}_R\cdot z<0\text{,}\end{array}$$

$$\alpha (z)=\{\begin{array}{ll}0,& \text{if}z\le z_{ba},\\ \frac{1}{2}\left(\sin \left(\pi \frac{\text{norm}z-\frac{z_{max}+z_{ba}}{2}}{z_{max}-z_{ba}}\right)+1\right),& \text{if}{z}_{ba}<\text{norm}z<z_{max},\\ 1,& \text{if}{z}_{max}\le \text{norm}z\text{.}\end{array}$$

$$g(V_R)=N_0\left({\mu }_k+({\mu }_s-{\mu }_k)\exp \left(-\frac{V_R^2}{V_S^2}\right)\right),$$

$$\{\begin{array}{ll}\frac{dX}{dt}& =\dot{X},\\ \frac{d^{2}X}{d{t}^2}& =-\frac{D}{M}\dot{X}-\frac{K_1}{M}{X}^3-\frac{K_2}{M}X+\frac{N_0}{M}{f}_R(V_R)+\frac{U_0}{M}\sin (\Omega t),\\ \frac{dz}{dt}& =\left(1-\alpha (z,V_R)\frac{{\sigma }_0}{g(V_R)}z\cdot \text{sgn}\left(V_R\right)\right)V_R\text{,}\end{array}$$

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})\text{,}$$

$${\mathbf{x}}_{i+1}=\mathbf{G}({\mathbf{x}}_i),$$

$${\mathbf{x}}_{i+1}+{\Delta }_{i+1}$$

$${\mathbf{x}}_{i+1}-{\Delta }_{i+1}$$

$$J\mathbf{G}({\mathbf{x}}_i){\Delta }_i$$

$$J{\mathbf{G}}_j({\mathbf{x}}_i)\approx \frac{\mathbf{G}({\mathbf{x}}_i+{\Delta }_i^j)-\mathbf{G}({\mathbf{x}}_i-{\Delta }_i^j)}{2\delta },$$

$$\{\begin{array}{ll}\dot{x}& =y,\\ \dot{y}& =z,\\ \dot{z}& =-ax-y-4z+y^2+xy.\end{array}$$

$${\mathcal{C}}_{MoB}(x)=\{\begin{array}{ll}1,& \text{if}{\text{LLE}}_{MoB}(\mathbf{x})\ge 0,\\ 0,& \text{if}{\text{LLE}}_{MoB}(\mathbf{x})<0.\end{array}$$

$$Y(\mathbf{x})=\mu +Z(\mathbf{x})\text{.}$$

$$Cov[Z({\mathbf{x}}^{(i)}),Z({\mathbf{x}}^{(j)})]={\sigma }^{2}R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)},\mathbf{\theta })\text{.}$$

$$\hat{\mu }=(1^T{\mathbf{R}}^{-1}1{)}^{-1}{1}^T{\mathbf{R}}^{-1}\mathbf{y}$$

$${\hat{\sigma }}^2=\frac{1}{m}{\left(\mathbf{y}-1\hat{\mu }\right)}^T{\mathbf{R}}^{-1}\left(\mathbf{y}-1\hat{\mu }\right)\text{,}$$

$$R(\mathbf{x}-{\mathbf{x}}^{\prime },\mathbf{l})=\prod _{i=1}^n\left(1+\frac{\sqrt{3}\text{absolutevalue}{x}_i-x_i^{\prime }}{l_i}\right)\exp \left(-\frac{\sqrt{3}\text{absolutevalue}{x}_i-x_i^{\prime }}{l_i}\right)\text{,}$$

$$\hat{\mathbf{\theta }}=\arg \underset{{\mathbf{\theta }}^{\star }}{\min }\psi ({\mathbf{\theta }}^{\star })\text{,}$$

$$\psi (\mathbf{\theta })={\hat{\sigma }}^2(\mathbf{\theta })[\det \mathbf{R}(\mathbf{\theta }){]}^{1/m},$$

$${\mu }_{{\hat{Y}}_0}=\hat{\mu }+{\mathbf{r}}_0^T{\mathbf{R}}^{-1}(\mathbf{y}-1\hat{\mu }).$$

$$r_{0i}=R({\mathbf{x}}^{(0)}-{\mathbf{x}}^{(i)},\mathbf{\theta })i=1,\dots ,m.$$

$${\sigma }_{{\hat{Y}}_0}^2={\sigma }^2\left(1-{\mathbf{r}}_0^T{\mathbf{R}}^{-1}{\mathbf{r}}_0+{\mathbf{u}}_0^T{\left(1^T{\mathbf{R}}^{-1}1\right)}^{-1}{\mathbf{u}}_0\right),$$

$${\mathbf{u}}_0=1^T{\mathbf{R}}^{-1}{\mathbf{r}}_0-1\text{.}$$

$${\mathbf{x}}^{(m+1)}=\arg \underset{{\mathbf{x}}^{\star }\in \mathbb{X}}{\min }\text{Score}\left({\mathbf{x}}^{\star }\right)\text{.}$$

$$EPE(\mathbf{x})=\alpha \hat{e}(\mathbf{x})+(1-\alpha ){\hat{\sigma }}_{\hat{Y}}^2,$$

Full text

1 Introduction

Stick-slip instability can be observed when two rigid bodies in contact slide on each other at low relative velocity. Then, intermittent vibratory effects may occur in form of variations of the frictional force and relative sliding velocity, with a characteristic sawtooth time-displacement curve, corresponding to stick and slip phases rabinowicz1956stick . These dynamic instabilities, also called stick-slip motion result in self-excited vibrations and drastic decreases in performance of some machine parts. Oscillating systems excited by dry friction are frequently encountered in many practical applications, including for instance brake systems (barton2004braking, ; ashraf2017investigation, ), hydraulic cylinders (owen2003reduction, ), gears wu2017comparisons or bearing rubio2007structural . Therefore the extensive analysis of these non-smooth dynamic systems to detect instabilities is a crucial point in engineering design.

Linear spring mass systems placed on a moving belt, commonly referred to as Mass-on-Belt (MoB) systems, have been generally utilized as simple mathematical models to investigate stick-slip system, see e.g. stelter1990stick ; hinrichs1997dynamics ; galvanetto1999dynamics ; jimenez2007two . To provide better accuracy, nonlinear mathematical representations have recently been considered (devarajan2017analytical, ). For instance, the bifurcation behavior of the non-linear MoB system has been investigated in santhosh2016discontinuity ; awrejcewicz2005stick . A Duffing’s type oscillator has been studied analytically in devarajan2017analytical in order to obtain expressions for stick-slip and pure-slip vibration amplitudes and frequencies. An estimation method for the spectrum of Lyapunov Exponents (LE) s proposed by (balcerzak2018spectrum, ) is employed in pikunov2019numerical to analyze the stability of a discontinuous MoB system.

To avoid computing and analyzing the behavior in the time domain, some indicators of instabilities have been proposed in the literature, such as Kolmogorov entropy benettin1976kolmogorov , correlation dimension grassberger1983characterization and sticking time lima2015stick . LE s have been stated by some authors as the most significant indicator for chaotic motion such as in (kocarev2006discrete, ). They quantify the rate of separation of infinitesimally close trajectories and the sensitive dependence on initial conditions exhibited by the system. In practical applications the determination of the Largest Lyapunov Exponent (LLE) is sufficient. If the number is non-positive, then the system is stable. On the contrary, a positive LLE indicates a system with chaotic behavior. Unfortunately, analytic determination of LE s is possible only for simple linear dynamic systems, see balcerzak2018spectrum , and numerical determination of their value is a challenge. In recent years several numerical strategies have been suggested for estimating LEs in case of non-smooth dynamic systems. A thorough review can be found in awrejcewicz2018quantifying . Recently an approximate numerical method based on the estimation of the Jacobian matrix by perturbation method of the initial orthogonal vectors involving an Euler forward scheme has been proposed balcerzak2018spectrum .

In addition, engineers want to prevent instabilities for a large scope of parameters, or scenarios which may be encountered during the working life of the mechanical systems. However the numerical analysis of the LLE in complex systems is computationally expensive, limiting the design space exploration and the optimization of a system at hand. Therefore statistical approximations or surrogate models emulating the LLE over a parametric input space appear very interesting to estimate the response on vast parametric space at low cost. Surrogate model approaches such as polynomial response, see e.g. kleijnen2017regression , kriging also known as Gaussian process regression (kleijnen2009kriging, ; williams1996gaussian, ), support vector regression (clarke2005analysis, ) or radial basis function models (park1991universal, ) have been developed and continuously improved in recent years. Kriging is an accurate interpolative Bayesian surrogate modeling technique, which is stated as the most intensively investigated surrogate approach jiang2017adaptive . Originally proposed by the geostatistician Krige krige1951statistical , it has been utilized for deterministic (sacks1989design, ) and random simulations models (van2003kriging, ). The idea of a kriging model is to provide an estimation of the response of the system on the whole parametric domain from the restricted available information provided by some observations. Two perspectives can be explored. One is generating an accurate surrogate model from the information given a priori. The second is to inform about the localization of crucial areas in the parametric space for the problem of interest, i.e. subdomains in which more information should be gained, if further exploration is possible. In light of recent contributions proposed for dynamical analysis and machine learning tools, the goal of this contribution is to propose surrogate strategies to investigate the oscillatory behavior and efficiently detect the risk of instabilities.

Thus, in this paper surrogate model construction to identify the parametric subdomains where dynamic instabilities may occur is investigated. Pertinent instability indicators are discussed, different metamodeling techniques are contemplated. Among them, kriging appears the most powerful for the problem of interest. Numerical challenges and capacities offered by kriging for the estimation of the sticking time and LLE analysis in a Duffing’s type oscillator with an elasto-plastic friction force model are exposed. More particularly, various adaptive schemes are evaluated in this context.

The paper is organized as follows. In section 2 a dynamic problem non-linear contact behavior is introduced, as well as the accompanying challenges due to the possible unstable behavior. Two possible indicators proposed by the literature to proficiently investigate the system behavior are exposed in section 3. A brief introduction to LE as well as most recent works which allow to estimate the LLE for discontinuous systems in a robust numerical framework is 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 the ability and challenge of the surrogate model approach to investigate unstable dynamic behaviors.

2 Dynamic problem of interest

The mechanical reference problem is a Duffing’s type oscillator with a damping term, as illustrated in Figure 1. A rigid body of mass $M$ is placed on a moving belt with a constant velocity $V_{0}$. The displacement of the body over time $t$ is denoted by $X(t)$. The body’s movement is restricted by a nonlinear spring of stiffness $K_{1}X^{2}+K_{2}$ and a dashpot with damping coefficient $D$ parallel to the spring. The relative velocity between the belt and the body is denoted by $V_{R}(t)$, which equals to $V_{0}-\dot{X}(t)$ with $\dot{X}(t)$ representing the first derivative of the body’s displacement with respect to time, i.e. its velocity. A time-dependent harmonic force $U(t)=U_{0}\sin(\Omega t)$ with amplitude $U_{0}$ and angular frequency $\Omega$ as well as a constant normal load $N_{0}$ is applied to the mass.

The equation of motion of the system reads

[TABLE]

where $\ddot{X}(t)$ is the second derivative of the body’s displacement with respect to time, i.e. the body acceleration and $F_{R}(V_{R})=N_{0}f_{R}(V_{R})$ is the friction force. $f_{R}(V_{R})$ is called the friction function which denotes the friction force per unit of normal load; it is a function chosen by the user. An overview of the common choices proposed in the literature is given in pennestri2016review . Generally, friction force models are classified either as static or dynamic (olsson1998friction, ), where in a dynamic model the friction force does not only depend on the relative velocity between the bodies but also on other state variables. A review of different friction force models for dynamic analysis has been presented in marques2016survey . Here, a dynamic model proposed in dupont2002single is employed. It belongs to the family of elasto-plastic models, which were initially introduced by prandtl1924spannungsverteilung to model the constitutive behavior of deformable bodies. Internal variables are used to represent the non-linear contact behavior. The displacement $X$ of the body is assumed to be the sum of an elastic, i.e. memoryless, contribution denoted by $z$ and a plastic, and so history-dependent, part $w$ such that

[TABLE]

A physical analogy of this model is depicted in Figure 2. During the sticking phase, the plastic displacement is constant, whereas while slipping the elastic component is assumed to be fixed.

The friction model, as illustrated in Figure 3, is based on the idea that any deformation of the asperities is associated with a corresponding friction force $F_{R}$. By analogy with plasticity, the history of the model is characterized by the internal variable $z$, which represents the bristle deflection. The friction function is defined as

[TABLE]

which can be compared to the force model of the classical LuGre model de1995new . The variable $\sigma_{0}$ denotes the bristle stiffness, $\sigma_{1}$ is the average bristle damping coefficient and $\sigma_{2}$ is a viscous component of the friction force.

Consider the following definition

[TABLE]

the velocity of bristle deflection $\dot{z}$ has been defined by dupont2000elasto as

[TABLE]

where $\alpha(z,V_{R})$ is a function, which is utilized to capture stiction effect. It depends on the relative velocity see e.g. townsend1987effect , since it defines the elastic deformation until the breakaway force of the system is exceeded, i.e.

[TABLE]

with

[TABLE]

The maximum bristle deflection is denoted by $z_{max}$, and the bristle deflection for 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 contact velocity increases stribeck1902wesentlichen . It is given by

[TABLE]

where $V_{S}$ is a parameter called the characteristic Stribeck velocity, and $\mu_{s}$ and $\mu_{k}$ denote the static and kinetic friction coefficients, respectively. Therefore, the complete equation of motion of the system in state-space form with the modified elasto-plastic friction model reads

[TABLE]

which can be solved, e.g. via numerical integration.

From the equations of motion, the dynamic behavior of the system can be analyzed. For instance, the time-response behavior plotted in Figure 4a shows the cyclic unstable behavior in time of the mass via the typical cyclic evolution of its displacement and its velocity. The cyclic unstable behavior can also be analyzed via the phase portrait as shown in Figure 4b. Analysis of the stick-slip motion using phase space is detailed in galvanetto1999dynamics .

The temporal analysis has been historically established. However, to automatically scan the parametric space and benefit from modern machine learning tools, it appears interesting to analyze the dynamic behavior via single scalar values, used as quantities of interest for the surrogate model.

3 Quantities of interest to analyze the dynamic system

Among the different scalar values proposed in the literature to investigate the dynamic system two are considered, namely the sticking time and the largest Lyapunov exponent.

3.1 Sticking time

The dynamic cyclic behavior as represented in Figure 4a can be decomposed into sticking and slipping modes. In lima2015stick it has been proposed to investigate the oscillation behavior of dry-friction oscillator using the percent of time in which the body remains in stick mode as Quantity of Interest (QoI). They suggest to improve properties e.g. endurance of mechanical systems by understanding the role of acting parameters on that QoI.

In the following, the sticking time of the MoB is defined as the time during which the relative velocity $V_{R}(t)=V_{0}-\dot{X}(t)$ remains below the threshold of $10^{-4}\,m.s^{-1}$. It is numerically estimated by integrating equation (9) with a six-stage, fifth-order, Runge-Kutta method with variable step-size. In order to avoid localized transient behavior, sticking time is analyzed on a time window between $150$ and $250$ seconds.

The analysis of this QoI is, however, not evident because as shown in Figure 5, the value of the sticking time cannot be directly associated with either stable or unstable behavior. Lima and Sampaio pointed out that both the duration of the sticking mode and the cumulative sticking time lack attention in the literature lima2015stick . However, this analysis is not the scope of this contribution, here the sticking time is simply considered as a candidate QoI to be represented by the metamodel.

3.2 Largest Lyapunov exponent

The second QoI aims at providing an instability indicator to break down the parametric space into stable and unstable domains.

Consider a general equation of motion of an $N$-dimensional time-continuous system recast into first order differential equation system of the form

[TABLE]

with the state-space vector $\bm{x}\in\mathbb{R}^{N}$ and $\bm{f}$ being a set of $N$ functions i.e. $\bm{f}=\left[f_{1}(\bm{x}),\ldots f_{N}(\bm{x})\right]$. Consider an $N$-dimensional sphere of initial conditions. With time the sphere evolves into an ellipsoid whose principal axes contract or expand with rates governed by the spectrum of LE s denoted by $\left\{\lambda_{i}\right\}_{i\in\left[1,N\right]}$. LE s introduced by oseledec1968multiplicative are defined as the average divergence or convergence rates of nearby orbits in state space. A positive exponent indicates local instability in this particular direction and hence characterizes chaotic motion, see rosenstein1993practical . The presence of a positive LLE indicates chaos, whereas negative values are characteristic 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 most useful tools to characterize the stability of a dynamic system in kocarev2006discrete . Therefore, it should turn out as a proficient QoI for building an univariate informative surrogate model.

However, analytic calculation of the LLE is mostly restricted to simple linear systems. For complex systems robust numerical techniques can be employed to estimate the LLE, see e.g. shimada1979numerical ,benettin1980lyapunov , or kantz1994robust .

3.2.1 Numerical estimation of the Lyapunov exponents

When the investigated dynamic system is smooth with an analytically obtainable Jacobian matrix, the LLE can be calculated with e.g. the algorithm given in wolf1986quantifying , in which the system of equations is solved for $N$ nearby initial conditions and correct state-space orientation is maintained by repeatedly orthonormalizing the corresponding set of vectors using the Gram-Schmidt procedure. However, the dry friction problem at hand is a non-smooth system. The Jacobian matrix cannot be determined or is strongly ill-conditioned, which leads to significant numerical problems when employing the algorithm proposed by wolf1986quantifying . To overcome this obstacle, a novel method has recently been presented in balcerzak2018spectrum , which estimates the Jacobian matrix by a truncated Taylor series expansion with small orthogonal perturbations. The spectrum of Lyapunov exponents has been estimated from that numerical method in pikunov2019numerical to analyze the stability of a nonlinear mass-on-belt system. It is shortly summarized as follows:

Consider a discretization in time given by the set $\{t_{i}|i=1,\ldots,n\}$. Let $x(t_{i})$ be denoted by $x_{i}$. The discrete dynamical system of the mapping $\bm{f}$ introduced in equation (10) at time $t_{i}$ reads

[TABLE]

where $\bm{G}(\bm{x}_{i})=\left[G_{1}(\bm{x}_{i}),\ldots,G_{N}(\bm{x}_{i})\right]$. Because the operator $\bm{G}$ in equation (11) may be discontinuous, an analytic expression for the Jacobian can not be given. However, by introducing a perturbation vector $\Delta_{i}=\left[\delta_{1},\ldots\delta_{N}\right]^{T}$, where each element is a scalar value of small magnitude, equation (11) can be rewritten in two ways

[TABLE]

Here $J\bm{G}(\bm{x}_{i})$ is the Jacobian matrix of $\bm{G}(\bm{x}_{i})$, which can be approximated as

[TABLE]

The equations correspond to forward difference and backward difference schemes respectively. Adding these two expressions yields an estimation of each column vector of the Jacobian as

[TABLE]

where $\Delta_{i}^{j}=\delta\bm{e}^{T}_{j}$ and $\bm{e}_{j}$ is the unit vector with unit in the $j$-th element. From this robust numerical estimation of the Jacobian matrix, any algorithm for calculating the LLE can be easily employed even for non-smooth applications. Here, the algorithm mentioned in wolf1986quantifying is considered.

3.2.2 Illustrative examples

The numerical approach is verified using a continuous three-dimensional system for which the analytic Jacobian matrix can be determined. The continuous case suggested by molaie2013simple reads

[TABLE]

Here $a$ is a given parameter. The bifurcation diagram of the problem for the domain $a\in[3.3,3.4]$ is plotted in the upper image of Figure 6. It can be seen that with an increasing value of the parameter the response gets more chaotic. The estimation of LLE using the algorithm proposed in wolf1986quantifying based on the exact Jacobian matrix of the system, is illustrated in the middle picture. The estimation using the same algorithm but based on the Jacobian matrix approximated by the numerical approach previously introduced in section 3.2.1 is plotted in the lowest image in Figure 6.

A perturbation vector $\Delta=10^{-4}\left[1,1,1\right]^{T}$ has been employed. It can be observed 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 bifurcation diagram well.

Let consider the discontinuous system of ordinary differential equations defined by equation (9). Here LLE values can only be evaluated from the numerical approximation. Using a perturbation vector $\Delta=10^{-4}\left[1,1,1\right]^{T}$, the bifurcation diagram and the respective LLE values for this dynamic problem are plotted over the angular frequency in Figure 7. It can be seen that the positive values of the LLE over the parametric domain correspond well with the bifurcation graph. Thus the LLE appears again as a pertinent instability indicator.

3.2.3 Largest Lyapunov exponent as an instability indicator

In order to accurately determine regular or chaotic motion of the system, it is not required to quantitatively estimate the largest Lyapunov exponent, but it seems rather more promising to use its non-negativity or negativity as a binary indicator.

From $\text{LLE}_{MoB}(\bm{x})$ which provides the LLE value for the MoB system with input $\bm{x}$ over the input parametric space, an instability indicator $\mathcal{C}_{MoB}(x)$ can be defined as

[TABLE]

For sake of illustration, the set of LLE values represented in Figure 8a is analyzed through this instability indicator, which is plotted in Figure 8b. The stable behavior associated with values of the indicator $\mathcal{C}_{MoB}=0$ is indicated in gray, whereas the parametric domain corresponding with chaotic behavior i.e. $\mathcal{C}_{MoB}=1$ is highlighted in red. It can be seen that this LLE surface provides a valuable information which could allow to optimize the system by preventing instabilities.

Based on this proficient QoI, the dynamic system could be efficiently analyzed through a surrogate model to evaluate the critical values at low cost. The goal of the surrogate model is to be able to provide an accurate approximation of this surface with only few observations, originating either from physical or numerical experiments. This surrogate model should be the best representation of the surface with respect to the available information, and provide guidance for defining further experiments to be conducted.

4 Construction of surrogate model

The idea behind surrogate models is only summarized here. Given a mapping $\mathcal{M}\,:\mathbb{X}\rightarrow\mathbb{Y}$ between an input $\bm{x}\in\mathbb{X}\subset\mathbb{R}^{n}$ and some uni-variate output $\bm{y}\in\mathbb{Y}\subset\mathbb{R}$. From a set of $m$ observations $\mathcal{D}=\{\left(\bm{x}^{(i)},\,\bm{y}^{(i)}\right),\,i=1,\,\ldots,\,m\}$, global surrogate modeling aims to statistically approximate the output on the whole input space, i.e. to approximate the statistical relationship behind this mapping. The samples in the continuous input space are called design of experiments and are denoted by the set $\mathcal{X}=\{\bm{x}^{(i)},\,i=1,\,\ldots,\,m\}$. For uni-variate output the mapping evaluations are gathered in the vector $\bm{y}=\{y_{i},\,i=1,\,\ldots,\,m\}$. The statistical approximation on the whole input space will be given by a computed surrogate model denoted by $\tilde{\mathcal{M}}$. Among the different existing surrogate modeling approaches, kriging is investigated here.

4.1 Kriging

Kriging was proposed by a mining engineer called Krige for geostatistics study (krige1951statistical, ). The mathematical foundation for the approach has been developed by matheron1963principles , which also established its name. For an overview the reader is referred to kleijnen2009kriging and kleijnen2017regression . It can thus be outlined that the approach can be described from a regression perspective as a statistical extension of deterministic regression stein1991universal , or from a bayesian or machine learning perspective as an update from prior information based on a conditional distribution williams1996gaussian .

4.1.1 Ordinary kriging metamodel

Using ordinary kriging kleijnen2017regression , the mapping between input and output data is approximated by a combination of a global contribution characterized by $\mu$, an unknown constant describing the mean contribution, and localized departures denoted by $Z(\bm{x})$

[TABLE]

In details, $Z(\bm{x})$ is a realization of a stochastic process with zero mean and a variance $\sigma^{2}$ accounting for local deviations. The covariance matrix of $Z(\bm{x})$ is given as

[TABLE]

Here $R(\bm{x}^{(i)},\bm{x}^{(j)},\bm{\theta})$ is the correlation function between two points of the input set. $\bm{\theta}$ is called hyperparameter and contains the unknown correlation parameters between two points. It can thus be seen, as illustrated in Figure 9, that the ordinary kriging model depends on two main contributions, a regression part which is here only a mean contribution, and a stationary Gaussian process based on the assumption of an auto-correlation function, which governs the deviation from that mean.

4.1.2 Prior information

Kriging model is based on some observations, as prior information within a Bayesian approach handcock1993bayesian . It could also be considered as a training stage of the surrogate model. Considered in ordinary kriging as a constant, the prior value of the mean $\hat{\mu}$ can be obtained from the available observation $\bm{y}$ using a generalized least-square estimate as

[TABLE]

with $\bm{1}$ a vector containing only components with scalar value $1$ and $\bm{R}$ the correlation matrix. It can be noticed that this prior evaluation depends on the chosen correlation function. The evaluation of the standard deviation from the prior information is denoted by $\hat{\sigma}$. Then, the variance of the prior information reads

[TABLE]

which, similarly to the mean, depends on the chosen correlation matrix.

4.1.3 Auto-correlation function and hyperparameter evaluation

The realization of the Gaussian process depends on its characteristics, namely its mean, its variance and its auto-correlation.

The auto-correlation function, commonly referred as correlation function for sake of simplicity, is usually assumed by the user. A common choice is the Matérn 3/2 correlation function (matern1960spatial, ), defined as

[TABLE]

in which a different hyperparameter $l_{i}$ can be used for each dimension $n$ of the design of experiments. Due to the lack of knowledge about the observations this correlation structure is arbitrary chosen here.

The hyperparameters need to be determined before constructing the surrogate model. This can for example be done by utilizing the maximum likelihood estimate of the hyperparameter value. This method yields an auxiliary optimization problem dubourg2013metamodel for estimating $\bm{\theta}$ of the form

[TABLE]

in which $\psi$ is the reduced likelihood defined for ordinary kriging as

[TABLE]

with $\det$ denoting the determinant operator.

Since there is no analytic solution for the optimization problem defined by Equation (23), it is necessary to use numerical optimization tools. Bouhlel and Martin mention this step as being the most challenging for the construction of surrogate models with kriging because of the multimodality of the likelihood function bouhlel2019gradient . In this work, a hybridized particle swarm optimization is employed similar to the method suggested by toal2011development .

4.1.4 Surrogate model evaluation

Consider the unobserved value $Y_{0}\equiv Y(\bm{x}^{(0)})$ with $x^{(0)}\in\mathbb{X}$, which needs to be predicted. When utilizing the best linear unbiased predictor for the Gaussian random variate $\hat{Y}_{0}$, the mean of the unobserved value yields

[TABLE]

$\bm{r}_{0}$ describes the cross-correlations between the new point $\bm{x}^{(0)}$ and each available observation with

[TABLE]

Furthermore the variance of the unobserved value can be estimated as

[TABLE]

with

[TABLE]

For proof and further details the reader shall be referred to santner2013design .

4.2 Adaptive sampling

To optimize the interest of the computations, the observation points which define the experimental design are chosen with respect to an adaptive sampling scheme. Starting from a set of initial data $\mathcal{D}_{ini}=\{\left(\bm{x}^{(i)},\,\bm{y}^{(i)}\right),\,i=1,\,\ldots,\,m\}$ a surrogate model $\tilde{\mathcal{M}}$ is constructed with kriging based on this knowledge. By solving some auxiliary optimization problems new points are successively added to the dataset $\mathcal{D}$ until a convergence criterion is reached. The general optimization problem reads

[TABLE]

Here a score function is evaluated which describes generally a tradeoff between an 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 of interest to reduce the prediction error locally. Different techniques have been proposed in the literature such as (jones1998efficient, ; turner2007multidimensional, ; singh2013balanced, ).

An external variance-based technique has recently been proposed in liu2017adaptive . Maximizing Expected Prediction Error (MEPE) is a continuous optimization involving a switch strategy between exploitation and exploration. The bias error is estimated by a continuously approximated leave-one-out cross validation (LOOCV) error $\hat{e}(\bm{x})$. The continuous expected prediction error to be maximized is given by

[TABLE]

where $\hat{e}(\bm{x})$ is the value of $\hat{e}(\bm{x}^{(i)})$ when $\bm{x}$ is located in the Voronoi cell of point $\bm{x}^{(i)}$. Otherwise its value is zero. In order to fasten the computation an approximation of $\hat{e}(\bm{x}^{(i)})$ can be used liu2017adaptive ; sundararajan2000predictive . A balance factor $\alpha$ is introduced to adjust the exploitative bias term and the exploratory variance term by switching adaptively between the two components depending on the estimation quality of the bias term. The computation of the balance factor for the calculation of the sample point $\bm{x}^{(m+q-1)}$ reads

[TABLE]

Here, $e_{true}$ represents the error between the surrogate model and the actual observation for this sample point. The algorithm is summarized in Box 4.2. {kasten}[htbp!] {Algorithmus}

• •

Given a design of experiments $\mathcal{X}=\{\bm{x}^{(1)},\ldots,\bm{x}^{(m)}\}$.

• While the adaptive sampling stopping criterion is not satisfied.

Do:

• Update the balance factor $\alpha$ with equation (30).
• Obtain the LOOCV error at each sample point $\hat{e}(\bm{x}^{(i)})$.
• Obtain a new point by maximizing the $EPE$-criterion of equation (29) over the input
• Update the design of experiments with the new found point and its response.

end

{kasten} Algorithm for MEPE

The second scheme considered here is Expected Improvement for Global Fit (EIGF) lam2008sequential , which is also a continuous optimization scheme, but with a fixed balance between exploitation and exploration. It is based on the maximization of the expected improvement defined as

[TABLE]

where $x_{\ast}$ is the closest sample point to the candidate $x$. The first term, which represents the exploitation contribution, is large when the surrogate model $\hat{Y}(x)$ differs substantially from the response at the nearest point. The second term provides global exploration, it is large when the surrogate model exhibits higher lack of knowledge.

For classification problems, a dedicated adaptive scheme called Monte Carlo-Intersite Voronoi (MiVor) has recently been proposed Class_Krig_sub ; Jan_Master_thesis . That scheme is based on the exploration of the parametric space by Voronoi tessellation. A score is given to all cells to identify the one cell corresponding to the largest lack of knowledge near to a decision region between two classes and so to pertinently add a new observation point. A random switching strategy is used to balance exploration and exploitation during the successive iterations. This allows to reach a compromise between the exploration of the whole parametric space and the addition in some local areas near to the class boundaries which require to be accurately described. The scheme has been described in detail, investigated for various applications and compared with other adaptive schemes in Class_Krig_sub .

5 Surrogate model for dynamic analysis

The surrogate strategy is tested on different applications. The problem parameters are summarized in Table 1. One-, two- and three-dimensional parametric spaces are considered.

The performances for both QoI s are analyzed with respect to a reference solution computed with $m_{ref}$ evaluation points $\{x_{ref,i},\,i=1,\dots,m_{ref}\}$. 5000 points per dimension are considered for the applications. Four indicators are considered to evaluate the error between the metamodel $\hat{Y}$ and the reference response surface $y_{ref}$. The Mean Absolute Error (MAE) is estimated considering the reference points as

[TABLE]

With the same perspective, the Root Mean-Squared Error (RMSE) is estimated as

[TABLE]

The Relative Maximum Absolute Error (RMAE) is given by

[TABLE]

with $\sigma_{y_{ref}}$ the standard deviation of the reference solution. Finally the $\text{R}^{2}$ score reads

[TABLE]

where $\overline{y}_{ref}$ denotes the mean of the reference response values. The parametric input space is normalized in each dimension in order to avoid numerical issues. Let $x_{i}^{u}$ and $x_{i}^{l}$ be the upper and lower limit of the dimension $i$ respectively. The normalized input $x_{i}$ is then given by

[TABLE]

The normalization ensures that $0\leq\bar{x}_{i}\leq 1$ and is denoted by a bar.

5.1 Surrogate model for sticking time: Problem $\mathcal{P}_{0}$

First, the ability of the surrogate model to provide an estimation for the sticking time is investigated. The parametric domain of the problem $\mathcal{P}_{0}$ is given by the two spring stiffnesses with $K_{1}\in[0.5,1.0]\,\text{N}.\text{m}^{-3}$ and $K_{2}\in[0.0,0.6]\,\text{N}.\text{m}^{-1}$. Other parameter values are considered to be deterministic, see Table 1.

The reference response surface for the sticking time over the input domain is plotted in Figure 10a. Initially 20 samples are spread across the domain with Translational Propagation Latin Hypercube Design (TPLHD) viana2010algorithm leading to a metamodel with an absolute error distributed as shown in Figure 10b. It can be noticed that the spread of the sticking time in the parametric domain $\{K_{1},K_{2}\}$ is around 5 seconds and the highest initial absolute error is around 1.6 seconds.

5.1.1 Adaptive kriging strategy

From the initial TPLHD model with 20 samples, surrogate models utilizing 60 samples are built considering the two alternative sampling techniques previously introduced, i.e. Maximizing Expected Prediction Error (MEPE)and Expected Improvement for Global Fit (EIGF). Furthermore a direct TPLHD model with 60 samples is generated. In order to avoid numerical outliers of the stochastic optimization process each strategy is executed twenty times, so that the 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 the initial and final one-shot metamodels based on TPLHD. Among them EIGF appears more proficient than MEPE.

In Figure 11 the convergence of MAE error and RMSE error with the number of samples is plotted for the two adaptive methods. As before, all results are averaged on twenty realizations. It seems that MEPE performs better as shown in Figure 11a. However, by considering the convergence in terms of RMSE error, both methods show a similar convergence behaviour, see Figure 11b.

A set of sample positions of the added points of the two adaptive methods after 150 samples are shown in Figure 12. On these plots the parameters have been normalized to lie between 0 and 1. The initial TPLHD sample points are shown as black circles. The later a point is added the more its color tends towards light red. The contour of the target function is displayed for both cases, it may be observed in Figure 12b that MEPE balances proficiently exploration and exploitation, wheras EIGF focuses on specific domains.

5.1.2 Comparison of kriging and other common metamodel techniques

In this subsection ordinary kriging is compared to three commonly found surrogate modeling techniques, i.e. support-vector machines, radial basis function network and neural networks. Here, the metamodels are here considered without adaptive scheme, and are all based on the same sample set provided by TPLHD. For this comparison the support-vector machine is employed with a radial basis kernel as well as an automatic choice of the scale value for the kernel function as defined in the MATLAB software, see e.g. clarke2005analysis . For the radial basis function network (park1991universal, ), the hidden layer has 50 neurons. The centers are calculated with k-means clustering algorithm and the widths are set to unity. The neural network is employed in form of multilayer perceptron based on 2 hidden-layer with 20 neurons each and sigmoid activation function where Levenberg-Marquardt backpropagation is used for training.

The mean absolute errors as defined by Equation (32) for the obtained metamodels considering different sample sizes are listed in Table 3. It can be noticed that the ranking of the methods depends on the number of samples considered. Support-vector machine performs better than radial basis function network and neural network for a low number of samples whereas it performs worse for number of samples larger than 100. Considering the parameters as detailed before, the results indicate that ordinary kriging yields the best approximation of the target function for all sample sizes for the problem of interest, while simultaneously provding an estimate for the variance of the metamodel over the parametric space.

5.2 Surrogate model for classification of chaotic motion using LLE

The second focus of interest is to predict the parametric subdomains corresponding with chaotic motion of the oscillator at low computational cost using the kriging approach. It has been shown in Jan_Master_thesis ; Class_Krig_sub that issues arise when trying to utilize the common adaptive sampling techniques for a regressive metamodel for that application. Indeed, as previously discussed, the precise value of the LLE is not of particular interest. The main focus of a metamodel for this application is to classify $\hat{\mathcal{M}}_{LLE,MoB}(\bm{x})$ the approximated estimation of $\mathcal{M}_{LLE,MoB}(\bm{x})$ for a point $\bm{x}$ of the parametric domain either above or below 0 to characterize the behavior of the given dynamic system as chaotic or regular motion. Using the adaptive sampling technique MiVor recently proposed in Class_Krig_sub , specifically dedicated to classification based on regression metamodel, is hereafter investigated to classify the LLE values. The required MiVor parameters (see Class_Krig_sub ), which are the initial exploration rate and the decrease factor, are chosen to be $0.4$ and $1.1$ respectively for all the following applications.

5.3 One-dimensional classification for stick-slip instability: Problem $\mathcal{P}_{1}$

The first case of interest is the one-dimensional problem $\mathcal{P}_{1}$, where the goal of the metamodel is to identify the unstable behavior corresponding with the LLE values plotted in Figure 13a. The angular frequency $\Omega$ varies between $0.2$ and $1.0\,\text{rad}.\text{s}^{-1}$, whereas the other input parameters have fixed values as defined in Table 1.

The performance of the metamodel is evaluated with respect to a reference solution obtained with 5000 sample points. The metamodel for classification is judged by how many of the points yielding $\mathcal{C}_{MoB}=0$ on one hand and how many points yielding $\mathcal{C}_{MoB}=1$ on the other hand are predicted. The measure is stated in percent of accurately predicted points. The percentage of correctly predicted points in the unstable regime is denoted by $a^{p}_{\text{LLE}\,\geq\,0}$ and defined as

[TABLE]

with $\hat{n}^{\text{LLE}\,\geq\,0}_{ref,\text{LLE}\,\geq\,0}$ the number of points among the ones predicted with unstable regime by the reference solution which are also predicted with unstable regime using the metamodel, and $\text{n}_{ref,\text{LLE}\,\geq\,0}$ the number of points estimated in the unstable regime by the reference response surface. Furthermore, the percentage of correctly predicted points in the stable regime $a^{p}_{\text{LLE}\,<\,0}$ is given by

[TABLE]

where $\text{n}_{ref,\text{LLE}\,<\,0}$ is the number of points with regular behavior given by the reference response surface. Besides $\hat{n}^{\text{LLE}\,<\,0}_{ref,\text{LLE}\,<\,0}$ denotes the number of points estimated with regular behavior using the kriging approach among the points which are predicted with regular behavior using the reference response surface.

The metamodel is built from the prior information given by five initial sample points set with TPLHD, i.e. $[0.2,0.36,0.52,0.68,0.84]^{T}$. None of them yields a chaotic output. The challenged posed to the adaptive metamodeling strategy is to evaluate if the scheme is able to explore the parametric domain to identify the chaotic subdomain.

The reference solution, the metamodel $\hat{M}$ after 35 samples and the respective sample 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 few observation points. Furthermore as intended with MiVor most of the created samples lie in the area that indicates chaotic behavior.

The evolution of the accuracy with regard to the sample size is shown in Figure 13b. As an arbitrarily chosen stopping criterion the procedure stops after 35 samples, i.e. after adding 30 observation points to the initial samples. Finally, the metamodel correctly evaluates $99.11\,\%$ of the 5000 random points that yield LLE $\geq 0$ and $99.23\,\%$ of the points with LLE $<0$.

5.3.1 Problem $\mathcal{P}_{2}$

Next, the metamodel approach is tested for the test problem $\mathcal{P}_{2}$ which aims at evaluating the ability of the proposed approach to proficiently explore the parametric domain if only a small percentage of the domain has a chaotic motion.

Consider the two-dimensional input domain for the two spring stiffnesses given by $K_{1}\in[0.5,1.0]\,\text{N}.\text{m}^{-3}$ and $K_{2}\in[0.0,0.5]\,\text{N}.\text{m}^{-1}$, see Table 1 for the values of other parameters. The plot of the reference response surface $\mathcal{M}_{LLE,MoB}$ provided by 10000 computations which were spread over the input domain in a space-filling manner using TPLHD is shown in Figure 14a. Large fluctuations on the LLE values can be observed. The analysis of the parametric domain through binary classification corresponding with either stable or chaotic motion is plotted in Figure 14b. As before, regular motion is represented in grey whereas the red area informs about chaotic behavior. It can be seen that only a small fraction of the domain yields chaotic motion. Among the 10000 points computed to build the reference solution only 93 points lead to stick-slip instability, i.e. only less than $1\,\%$ of the points.

The metamodel is built from an initial set of 10 initial experiments created with TPLHD, then adaptive sampling is used until reaching 65 samples. The point locations of the 55 samples created with MiVor over the normalized input space are highlighted in Figure 14c. It can be observed that they are predominantly spread in and around the red area that indicates chaos. This leads to an accurate surrogate model $\hat{\mathcal{M}}_{LLE,MoB}$ for the classification problem as shown in Figure 14d. In details, it yields that $98.93\,\%$ of the points above or 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. It can be seen that $a^{p}_{LLE\geq 0}$ MiVor shows a big jump around 15 sample points. This is due to the fact that the algorithm needs to identify the small area of Figure 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^{p}_{LLE<0}$ appears almost accurate from the beginning with values varying between $99.6-100\,\%$.

5.3.2 Problem $\mathcal{P}_{3}$

Problem $\mathcal{P}_{3}$ aims at evaluating metamodel performance to detect almost disconnected regions of the input domain associated with chaotic behavior.

Consider the input domain for the spring stiffness given by

$K_{2}\in[0.0,0.5]\,\text{N}.\text{m}^{-1}$ and the kinematic friction coefficient $\mu_{k}\in[0.08,0.18]$. Values of all other parameters are listed in Table 1. The LLE reference surface established from 10000 sample points is plotted in Figure 16a. It exhibits large fluctuations. This reference solution yields the classification displayed in Figure 16b. It can be seen that there are two chaotic subdomains which are almost disconnected.

The metamodel is constructed adaptively from 5 points sampled with TPLHD. It is evaluated after adding 100 sample points and the results are averaged over 20 independent realizations. The locations of the 105 MiVor sample points for one realization are shown in Figure 16c. It can be noticed that the method is effective in sampling around the red zones indicating chaos. Furthermore the rest of the domain is sampled with a space-filling design. The resulting metamodel after 105 samples is displayed in Figure 16d and shows an accurate prediction performance. The percentage of correctly identified points for LLE$\geq 0$ is $98.53\%$, whereas $99.19\%$ of the points corresponding to LLE values below 0 are correctly identified. The convergence of the percentage values is displayed in Figure 17. It shows a very proficient sampling performance of MiVor.

5.3.3 Problem $\mathcal{P}_{4}$

The metamodel approach for identifying stick-slip instability is tested on Problem $\mathcal{P}_{4}$, characterized by a highly complex response surface, in which regular and chaotic motion subdomains are embedded into each other. Consider the parametric domain for the two spring stiffnesses given by $K_{1}\in[0.5,1.0]\,\text{N}.\text{m}^{-3}$ and $K_{2}\in[0.0,0.6]\,\text{N}.\text{m}^{-1}$ with other parameters as given in Table 1.

The reference LLE response surface based on the evaluation of 10000 sample points is shown in Figure 18a. The corresponding classification is pictured in Figure 18b. It can be noticed that the chaotic area is complex with some regular subdomains embedded within chaotic behavior.

A set of 5 points sampled using the TPLHD technique is considered as initial observation points. One realization is represented in Figures 18c and 18d. The adaptive sampling technique is evaluated after reaching 115 sample points. Samples needed to generate the metamodel are shown in Figure 18c. The resulting metamodel is represented in Figure 18d. It can be seen that the holes inside the red domain are not detected by the surrogate model. This reduces the percentage of correctly classified points. However an accurate prediction of these holes with kriging would require a significantly larger number of samples. Furthermore identifying a point inside the regular motion subdomain might lead to a perceived reliability of the system dynamics around the observation point which might not be desired from an engineering point of view.

The convergence of the error metrics are displayed in Figure 19, which shows an efficient convergence behavior. The results of the MiVor metamodel after 115 samples are very proficient with values of $a^{p}_{LLE\geq 0}=97.12\,\%$ and $a^{p}_{LLE<0}=97.32\,\%$. Considering the very complex form of the chaotic region including notched boundary and holes as shown in Figure 18b, these percentages of accurately identified points computing only 115 observation points are highly promising for detecting non-stable dynamic parametric domains at low cost.

5.4 Problem $\mathcal{P}_{5}$

Finally the input dimension is extended to the three-dimensional problem $\mathcal{P}_{5}$. Consider the parametric space $\Omega\times K_{2}\times\mu_{K}\in[0.6,0.9]\,\text{rad}.\text{s}^{-1}\times[0.0,0.5]\,\text{N}.\text{m}^{-1}$ $\times[0.1,0.15]$ with $K_{1}=1\,\text{N}.\text{m}^{-3}$. The rest of the parameters remain unchanged with regard to the previous examples as given in Table 1. As reference 15000 observations, which are evenly spread in the parametric space, are evaluated.

The MiVor process is started from an initial set of 5 TPLHD samples, none of which correspond to an unstable dynamic system. The evolution of the averaged error values until including $200$ samples is displayed in Figure 20.

It can be noticed that MiVor achieves proficient value convergence with both metrics indicating an exactitude of above $98\%$ of accurately identified points for both classes with already around 125 observation points. The jittering of the values indicates sharp drops of the LLE response surface. It can be seen that MiVor can be very well expanded to higher dimension on that specific example.

6 Conclusion

Using the latest development proposed in the literature both in terms of numerical approaches allowing to compute the LLE and kriging models for classification, a proficient analysis of the parametric model can be conducted to optimally improve the knowledge of dynamic instabilities. A non-linear oscillator 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 introduced as means to investigate the system behavior. The analysis of this system which appeared highly challenging, as both QoIs have highly fluctuating values in the parametric domain, is very efficiently tackled using an adaptive metamodeling approach. The LLE appears as the most promising QoI in that context. Values equal or above zero indicate chaotic motion accurately. The MiVor sampling technique is able to generate surrogate models able to classify chaotic and regulate motion with only a few observations. It efficiently gains information from the available observation to further guide the investigation and thus allows to identify the instability domain(s) with reasonable number of computations. It therefore provides pertinent information for design optimization, even for complex response surfaces and two- or three-dimensional parametric spaces.

Acknowledgements.

The authors acknowledge the financial support from the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy within the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453).

The results presented in this paper were partially carried out on the cluster system at the Leibniz University of Hannover, Germany.

Conflict of interest The authors declare that they have no conflict of interest.