The basin stability of bi-stable friction-excited oscillators

: Stability considerations play a central role in structural dynamics to determine states that are robust against perturbations during the operation. Linear stability concepts, such as the complex eigenvalue analysis, constitute the core of analysis approaches in engineering reality. However, most stability concepts are limited to local perturbations, i.e. they can only measure a state’s stability against small perturbations. Recently, the concept of basin stability has been proposed as a global stability concept for multi-stable systems. As multi-stability is a well-known property of a range of nonlinear dynamical systems, this work studies the basin stability of bi-stable mechanical oscillators that are affected and self-excited by dry friction. The results indicate how the basin stability complements the classical binary stability concepts for quantifying how stable a state is given a set of permissible perturbations.


Introduction
The dynamics of systems affected by friction are most often studied in the context of friction-excited vibrations (FIV). Prominent examples for FIV in mechanical structures and machines range from brake systems [1][2][3][4], clutches [5], drill strings [6] to artificial hip joints [7] and others. FIV often arise through positive energy feedback from a friction interface with the structure, i.e. through self-excitation [8][9][10]. Sub-critical Hopf bifurcations [11,12] and isolated solution branches [13][14][15] are a common observation in those systems, such that bi-and multi-stable systems have been reported numerously [14,16,17]. The computation of those nonlinear responses (periodic, quasi-period orbits, chaotic trajectories) is a well-established field of research [18][19][20][21], mostly resulting in the identification of complicated bifurcation diagrams [11,13,[22][23][24]. The stability of the solutions is usually assessed by local Lyapunov-type stability metrics [25,26]. Hence, the stability statement is often a binary one that measures the state's robustness against small perturbations. However, the actual size of permissible perturbations, i.e. those for which the trajectory would still return back to the state, is not given. In a multi-stable system configuration, the long term steady-state behavior thus depends on the choice of initial conditions or the size of instantaneous perturbations. Once the system enters another basin of attraction, severe jumping phenomena may occur. Typically, such jumps are related to a change from a stable steady sliding state to high-intensity periodic vibrations or stick-slip cycles [27,28], or from one periodic solution to another periodic solution [11].
This work investigates a rather novel technique denoted as basin stability to estimate the size of the system's basins of attraction in a subset of the state space. The basins' size estimation can then be considered a global stability metric, i.e. indicating how likely the system is to end up on one of the co-existing stable steady-states. Therefore, those probabilities add new insights to the rather binary stability statements derived from local perturbation-based approaches. We study a friction oscillator excited by a falling friction slope and a second oscillator excited through binary flutter instability.
Our results indicate that the basin stability analysis is a robust and easily applicable model-agnostic technique. It can reveal the actual picture of the long-term behavior for a given set of perturbations, thus augmenting classical bifurcation and stability studies. Using the basin stability analysis, some solutions can even be ruled out if one can guarantee strict control over the instantaneous perturbations to system trajectories or operating conditions.

The concept of basin stability
We study nonlinear dynamic systemṡ with the states x (t) in the N-dimensional state space. The long-term asymptotic behavior is denoted as attractor A [29] denotes the subset of states that converge to the same attracting set A. The basin boundaries are related to unstable solutions of the system which represent separatrices of the basins in state space. Depending on the size and shape of its basin, an attractor can be more or less robust against non-small perturbations. There are multiple ways to compute the basins of attraction, e.g. through Lyapunov functions [30]. These methods are known for some canonical, low-dimensional, and well-studied systems. However, they are not readily available, or straight-forward to compute, for any generic and high-dimensional nonlinear dynamical system, such as frictional oscillators which are studied in this work. The basin stability proposed by Menck et al. [31,32] is a global stability concept for complex systems that aims at measuring stability against non-small perturbations by a volume-based probabilistic approach. Conceptually, the basin stability measures the volumetric share of all basins of attraction in a hypervolume of the state space. For a computationally feasible solution, a distribution ρ (x) of perturbations is drawn from a reference subset Q ⊂ R N of the state space, representing a set of states to which the system may be pushed to through non-small perturbations with Q ρ (x) dx = 1. For each perturbation, the steady-state behavior of the dynamical system is obtained through time-marching integrations. Then, the fraction of perturbed states that converged to the specific attractor A denotes an estimate for the basin stability S B (A), i.e. [31,33] Here, κ B denotes an indicator function that classifies a steady state solution x (t) to belong to the attractor A. Therefore, S B (A) is an estimate for the volume share of the basin of attraction B A given the reference subset Q ⊂ R N sampled by ρ (x) [31,34]. Naturally, for a k-multi-stable system, the basin stability values of all k attractors add up to unity ∑ k i S B (A i ) = 1. As the basin stability computation can be considered a repeated Bernoulli experiment [31], the standard error of the basin stability estimate is which can be used to find a subset Q that ensures a low standard error. Recently, systems with fractal, riddled, and intermingled basin boundaries were studied [32] indicating the robustness of the basin stability concept. All basin stability computations in this work were obtained from the open-source package bSTAB [35] available at https://github.com/TUHH-DYN/bSTAB/tree/v1.0.
x LC x EP x 1 n randomly chosen states x  Figure 1 displays a schematic for illustrating the practical computation of basin stability values. A nonlinear dynamical system with two states x = [x 1 , x 2 ] is studied 1 . The system exhibits three solutions: A stable equilibrium position (x EP ), an unstable periodic orbit (x UPO ), and a stable limit cycle (x LC ). The distribution of perturbations ρ (x) is chosen such that all solutions are contained in Q and n = 100 samples are drawn uniformly at random. The trajectories starting from n EP = 37 states in the basin B EP converge towards the equilibrium position, while n LC = 63 states are located in the basin B LC and thus converge to the stable limit cycle. As a result, the basin stability estimates are S B (EP) = 0.37 and S B (LC) = 0.63, respectively. Because the separatrix, that is the unstable periodic orbit, is explicitly known for the system, the basin volumes can be determined analytically. The exact volumetric fractions of B EP and B LC in Q are 0.3275 and 0.6725, respectively. Therefore, the basin stability computed from n = 100 samples is a good approximation for the system at hand 2 . 1 In fact, the system is the single-degree-of-freedom frictional oscillator to be discussed in Section 3 2 Appendix A.2 indicates that n ≈ 300 samples are required for a very close approximation of the analytical results

Bi-stable oscillator with falling friction slope
As a first system, we study the dynamics and the stability of a single-degree-of-freedom oscillator mẍ + cẋ + kx = F, see Figure 2 (a), with velocity-dependent friction as proposed by Papangelo et al. [12]. Specifically, the friction characteristic µ (v rel ) is a velocity-dependent weakening function featuring the static friction coefficient the reference velocity v 0 and the contact normal load N. The non-dimensional form of the equations of motion is obtained through normalization (·) of the quantities accordingly to the work of Papangelo et. al [12]. The velocity-dependence introduces a dynamic instability that gives rise to friction-induced vibrations (FIV) for 0 ≤ṽ d ≤ 1.84. Moreover, the friction nonlinearity enables the system to exhibit a bi-stable behavior, such that a stable steady sliding state and a stable stick-slip cycle co-exist for a range of belt velocities 1.11 ≤ṽ d ≤ 1.84, see Figure 2 (b). Atṽ d = 1.15, the steady sliding state loses stability at a subcritical Hopf bifurcation point. In the bi-stability regime, and depending on the initial condition or instantaneous perturbations, the system will either end up in the low-energy steady sliding state, or on the high-intensity stick-slip cycle. Both solutions are locally stable and attractive, i.e. robust against small perturbations.

Bi-stable oscillator with mode-coupling
As a second system, we study a frictional oscillator [8,11], which (in contrast to the first system) experiences FIV through a mode-coupling instability. The system features a main oscillating mass that is in dry Coulomb-type frictional contact with a conveyer belt. A second mass is connected to the main mass through a nonlinear joint element in diagonal direction, thereby geometrically coupling the vertical and horizontal movement of the main mass. The relative sliding velocity is assumed to always be positive, such that no stick-slip cycles can arise. For the nonlinear joint element, a cubic stiffness nonlinearity k nl is chosen [11]. The equations of motion and parameter values are given in Appendix B and the model is displayed in Figure 4. Previous studies have revealed the complicated bifurcation behavior of this system, including super-and sub-critical Hopf bifurcations as well as isolated solution  branches [11,13,14]. In this study, a variation of the horizontal stiffness k x is performed. A sub-critical Hopf bifurcation point is found at k x = 32.3, see Figure 5 (a). Below, a stable limit cycle and the unstable fixed point exist. Above this value there is a bi-stable range up to k x = 33.0 with a co-existing stable limit cycle and the stable fixed point. The eigenvalues' real parts in Figure 5 (b) exhibit the classical forking behavior that is related to the mode-coupling instability mechanism in this system. At the point of instability, one eigenvalue crosses into the positive plane. The basin stability S B of both stable solutions is computed for n = 500 random initial conditions drawn from Q (x,ẋ) : [0, 0.5] × [0, 0.25] (all other initial conditions are fixed to 0). Figure 5 (c) depicts the basin stability as a function of the horizontal stiffness. In the bi-stability range 32.3 ≤ k x ≤ 33.0 the basin stability values indicate that the limit cycle solution is the dominating one for lower stiffness values. For larger stiffness values the fixed point solution is the most probable for our choice of Q. Hence, within this rather short bi-stability range, a minor variation of the horizontal stiffness value would crucially affect the probability of arriving either on the low-energy steady-sliding state, or on the high-energy limit cycle, which may cause increased wear, audible vibrations and other effects in realistic systems. Such kind of statement about the global stability w.r.t. non-small perturbations would not have been easily accessible through the bifurcation diagram or the local stability analysis. These results are clearly related to the shape of the unstable periodic orbit, i.e. the separatrix of both basins of attraction. While the qualitative basin stability values for a variation in the initial conditions for x would have been readable from the bifurcation diagram in Figure 5 (a), this task quickly becomes complex once more degrees-of-freedom (DOFs) are considered. For example, let us consider a reference subset Q that captures certain variations for multiple DOFs, instead of variations for a single DOF as shown before. Figure 6 displays the basin stability values for three different choices of Q, i.e. different variations of the range of possible initial conditions: In the first case, some initial variations in the horizontal position and large variations in the vertical displacement of the main mass are allowed. In the second case, variations in the initial velocities are studied, and in the third case also variations in the secondary mass' initial conditions are considered. Such scenarios would quickly become impractical for studying permissible perturbations, i.e. the global stability of each solution, using bifurcation diagrams and subdividing the state space by the unstable solutions. The concept of basin stability automates this process through the Monte Carlo sampling, allowing for a easy-scaling and consistent estimation of the relevant basin volumes. In fact, even though the three reference sets are very different in their value ranges, the resulting basin stability analysis displayed in Figure 6 does not change qualitatively. The turning point, i.e. the point after which the FP solution dominates over the LC solution for increasing values of k x , changes only slightly: For Q 1 this point is found at k x = 32.9, while it is k x = 32.7 and k x = 32.55 for Q 2 and Q 3 , respectively. Hence, the basin stability is not very sensitive to the choice of Q for this system. In a situation in which the overall qualitative behavior of the basin stability values may have seem obvious, the quantitative evaluation would have become difficult to obtain from the bifurcation diagrams. Especially for higher-dimensional systems and specific subset choices the basin stability analysis represents a highly robust approach to estimate the probability of arriving on either of the competing solutions, which we will illustrate in the next section.

Bi-stable oscillator with isolated periodic solution
The third dynamical system studied in this work is a weakly-damped variant of the system proposed in the previous section and sketched in Figure 4. Here, the damping parameters are reduced by a factor of 10 to d x = d y = d z = d lin = 0.002. This system configuration has already been studied in [13,14] where the authors found an isolated solution branch resulting from the damping variation. Figure 7 (a) displays the bifurcation diagram for the horizontal stiffness k x . The fixed point solution loses stability through a sub-critical Hopf bifurcation at k x = 32.24 to a limit cycle solution, hereafter denotes as LC1. Interestingly, a second stable limit cycle solution is born for k x < 29.9, which is found to be an isolated branch [14], hereafter denoted as LC2. That is, this solution is not connected to Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 November 2020 doi:10.20944/preprints202011.0234.v1 any other solution path. As a result, the system may jump from the fixed point solution to the first limit cycle for 32.24 ≤ k x ≤ 33.0, and then from the limit cycle to the isolated branch for k x < 29.9. Hence, within a rather narrow parameter range, two jumping phenomena between different solutions may occur. It is therefore of great interest to investigate the probability of arriving on either of those solutions for some prescribed set of initial conditions. For for the bi-stability range featuring the two periodic solutions LC1 and LC2 (k x < 29.9) the basin stability analysis reveals that LC1 is the by far most probable solution. A maximum of 21% of the trajectories converge to the isolated solution branch, while the remaining trajectories converge to the first periodic orbit. Particularly interesting is the parameter regime 27.4 ≤ k x ≤ 29.9. Here, the basin stability indicates that LC1 is globally stable, even though the stable isola still co-exists. However, due to the choice of Q, no initial conditions were drawn for the basin related to LC2. Hence, if the range of initial conditions and perturbations can be quantified or limited for some specific system, the basin stability analysis can also help to rule out jumping phenomena between co-existing solutions.
Another interesting observation is the following: the basin stability values in this specific setting do not follow the qualitative trend of the respective amplitudes reported in Figure 7 (a). S B (LC1) keeps increasing along the stiffness parameter, while the corresponding amplitude of the horizontal vibration amplitude shows a different behavior. Theoretically, it is clear that the vibration amplitudes do not relate to the size of the basins of attraction. However, on the first sight classical bifurcation diagrams may suggest that one solution is more attractive if it has a larger vibration amplitude. At this point, the basin stability represents a technique to quantify the attractiveness in a highly consistent manner.
Lastly, we discuss our previous thought on the benefits of having a robust methodical approach to estimating the basin volumes through Monte Carlo sampling irrespective of the dynamical system at hand (so-called model-agnostic techniques). Especially for such low-dimensional systems as shown before, one might raise the issue of using computation-heavy sampling methods, even though the basins of attraction are readily available once the bifurcation diagram is known. Figure 8 displays the state space of each DOF at k x = 27, hence in a configuration where the two periodic orbits co-exists. It becomes clear that even for this 3 DOF oscillator (6 states), the analytical calculation of the basin volumes can quickly become a challenge. There is no straight-forward way to computing the volumes in the six-dimensional space from the intertwined basins separated by the unstable orbits, especially looking at the z coordinate. Therefore, the basin stability analysis is not only relevant for systems featuring larger number of states, but also for rather low-dimensional systems.

Conclusions
This work proposed to augment the classical local stability analysis of friction-excited oscillators by their basin stability. The concept of basin stability allows to assign global stability metrics to multi-stable solutions in a highly automated manner including error estimates. For three different friction-excited systems we show that the knowledge about global stability with respect to a specific set of initial conditions can provide important insights into the long-term dynamics. Particularly for well-controlled perturbations, this approach allows to estimate the probabilities of arriving on either of multiple stable solutions, and even to rule out some steady-state behavior. As a result, we suggest to include the basin stability analysis into the toolbox of techniques that are applied to study the nonlinear dynamics of multi-stable systems, especially when operating conditions are well-known.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: where (·) indicates a non-dimensional quantity.

Appendix A.2. Convergence of basin stability values
The number of samples n is varied to answer the question of how many samples from Q are required for a robust approximation of the basin stability values. Figure A1 displays the convergence of the basin stability values for the single-DOF oscillator case and the corresponding analytical values. where u is the relative displacement in the joint between the main mass and the secondary mass, given by u = − √ 2 2 x − √ 2 2 y + z. The parameter values for the reference configuration are given by m = m 1 = 1, k x = 32.5, k y = 20, k z = 100, k lin = 10, k nl = 5, d x = d y = d z = d lin = 0.02, µ = 0.65.