Robust identification of backbone curves using control-based continuation.

Control-based continuation is a recently-developed approach for testing nonlinear dynamic systems in a controlled manner and exploring their dynamic features as system parameters are varied. In this paper, control-based continuation is adapted to follow the locus where system response and excitation are in quadrature, extracting the backbone curve of the underlying conservative system. The method is applied to a single-degree-of-freedom oscillator under base excitation, and the results are compared with the standard resonant-decay method.


Introduction
The constant drive to improve the performance of engineering structures is increasingly leading to lighter and more flexible designs where nonlinearity is inherent. While nonlinearity is often viewed as detrimental, recent contributions in the literature have shown it can actually be exploited for improving system performance. For instance, nonlinearity was deliberately introduced in the design of sharp acoustic switches and rectifiers [1], vibration absorbers [2] and energy harvesters [3]. The presence of nonlinearity however poses new challenges to engineers because, in contrast with linear systems, nonlinear systems can exhibit a wide variety of complicated dynamic phenomena such as intermittency, quasi-periodic oscillations, chaos, and bifurcations. A problem of particular interest to engineers is the prediction of the response of a system at resonance, where the system is at great risk of failure. Pioneered in the 1960s by Rosenberg [4], the concept of nonlinear normal modes (NNMs) is considered as the natural extension of linear normal modes (LNMs) to nonlinear systems. Oscillations in nonlinear systems are energy-dependent such that the resonance frequency generally varies with the oscillation amplitude. NNMs can trace out this evolution, thereby generating the so-called backbone curve. For systems with moderate linear modal damping, the NNMs of the underlying conservative systems generally describe well the evolution of the resonance frequencies of the damped forced system. NNMs of conservative systems are defined as families of non-necessarily synchronous periodic oscillations [5]. There exist both analytical [6,7] and numerical [8,9] methods to calculate NNMs from a mathematical model. The latter are quite sophisticated, see for example the review [10], and can address complex real-life structures such as a full-scale spacecraft structure [11].
Experimental extraction of modal properties plays a central role in the updating and validation of linear structural models. In the context of nonlinear systems, the system's forced response and backbone curves can be used to estimate model parameters [12][13][14] and apply model updating techniques [15][16][17]. The experimental identification of NNMs was proposed in Refs. [18,19]. Following the principle of linear phase separation techniques, the method isolates a single NNM using an appropriate excitation. The applied force is then stopped and the relation between amplitude and frequency of oscillation is extracted from the free, damped, response of the system, a method often termed resonant decay. The method was successfully applied to several academic systems of moderate complexity as such, a single-degree-of-freedom oscillator [20], a nonlinear beam [19] and a steel frame structure [21]. A phase separation method where multiple NNMs are identified simultaneously from broad-band data was introduced and demonstrated on noisy synthetic data in Ref. [22].
The present paper proposes a novel method for extracting experimentally the underlying NNMs of a forced, damped system. The proposed method is robust to bifurcations and stability changes that arise in the tested system dynamics, and differs from existing contributions, such as resonant decay, in that the backbone curve is no longer post-processed from the measured data but rather directly traced out in the experiment. To this end, the steady-state periodic solutions that describe the backbone curve and a NNM of the underlying conservative system are followed for increasing vibration amplitudes using the control-based continuation (CBC) method. CBC combines stabilizing feedback control and path following techniques to explore the dynamics of a nonlinear system directly during the physical experiment, tracking the evolution of its steady-state response as system parameters are varied. The fundamental ideas underlying CBC were introduced by Sieber and Krauskopf in Ref. [23]. The first experimental demonstration of the method was performed on a parametrically-excited pendulum whose periodic solutions were followed through a saddle-node bifurcation, beyond which point the solutions become unstable [24]. The frequency response of a harmonically excited impact oscillator was studied in Ref. [25] and Barton et al. investigated the periodic solutions of two energy harvesters in Refs. [26,27].
The paper is organized as follows: Section 2 reviews the connection that exists between NNMs and the forced response of a damped system. Section 3 discusses the identification of NNMs within an experiment. The so-called resonant-decay method currently used for extracting NNMs is first presented in Section 3.1. The CBC approach developed in this paper is then introduced and adapted to track the steady-state periodic solutions that define the backbone curve. A single-degree-of-freedom (SDOF) set-up is considered for demonstrating the proposed approach. It is presented in Section 4.1, and experimental results are discussed in Sections 4.2 and 4.3. A mathematical model of the SDOF oscillator is derived in Section 5 in order to further elaborate on the comparison between the proposed CBC approach and resonant decay. Conclusions are drawn in Section 6. Peeters et al. showed that a forced damped system can follow precisely one NNM motion of the underlying conservative (unforced) system provided that an appropriate excitation is applied [18]. Consider a multi-degree-of-freedom (DOF) system with stiffness nonlinearities, the equations of motion can be written as

NNM motions in the presence of damping
where M, C, and K are the linear mass, damping and stiffness matrices, respectively. The displacement, velocity and acceleration vectors are x,ẋ, andẍ. The external loads are f (t). The nonlinear force vector is f nl (x) and it is assumed to contain odd nonlinearities only.
Consider the multi-point multi-harmonic excitation where F k corresponds to the fully-populated amplitude vector of the force for the k th harmonic frequency. The system can be studied when the response is in quadrature with the excitation, i.e. when x(t) is 90 • out of phase with f (t). With this assumption, x(t) can be decomposed as a cosine series where X k are real vectors. The nonlinear force term can be decomposed in a similar cosine series with coefficients F nl,k . Substituting Eqs. (2) and (3) into Eq. (1), each sine and cosine terms in each harmonic can be balanced to give the following relations These equations reveal that if the harmonics of the excitation are all in quadrature with the harmonics of the response, the excitation exactly cancels out the damping force (Eq. (5)) and the periodic response will exactly satisfy the undamped equations of motion, as shown by Eq. (4), which by definition is an NNM motion. If general polynomial nonlinearities in stiffness are considered, Eq. (3) has to include an additional constant term (k = 0) to account for the constant term arising in the Fourier decomposition of the nonlinear force. However, Eqs. (4)(5) remain unaffected because no sine term is introduced, and the damped forced system can therefore still vibrate according to the NNM motion of the underlying conservative system.
Contrary to linear systems that respond only at the forcing frequency, the response of nonlinear systems can comprise many harmonics. The phase quadrature condition should therefore be satisfied for each harmonic. In addition, Equation (5) shows that the excitation should be distributed in space (i.e. applied to each DOF) in order to cancel out damping. This condition is unlikely to be met in practice where the excitation is usually limited to a few number of DOFs, i.e. F k contains only a few non-zero elements. The influence of inadequate forcing distribution was studied in Ref. [28] using the second-order normal form technique. The lack of appropriate forcing was characterized by a phase difference and a net energy transfer between DOFs. This observation shows that the quadrature condition and the underlying NNM motion might never be reached depending on the forces applied to the structure. Force distribution was found particularly critical for closely separated and internally resonant modes. In practice, several studies [19,21] have however shown that much simpler excitations, single-harmonic and single-point, can still isolate the NNMs of a system with well-separated modes. For simple systems, the excitation can be further simplified and conveniently replaced by initial conditions in displacements. Finally, the nature of damping is generally not linear as assumed in Eq. (1). The presence of nonlinear damping can also prevent from reaching the quadrature condition, depending on its form and its importance in the response of the system.

The resonant-decay method
The idea of extracting the amplitude-dependent characteristics of a system from its free (resonant) response is not new. In Ref. [29], Braun and Feldman used the Hilbert transform to study the characteristics of a simulated Duffing oscillator. In Ref. [30], a nonlinear beam was analysed using the wavelet transform. Peeters et al. combined time-frequency analysis and the phase-quadrature result of Section 2 in a two-step methodology for extracting experimentally NNMs [18,19]. The methodology is similar in principle to the linear resonant decay method proposed in Ref. [31]. The first step consists of finding an appropriate excitation force to isolate the NNM of interest. To this end, an harmonic excitation is applied to the system and the forcing frequency is gradually incremented to reach phase quadrature with the response. At quadrature, the system vibrates according to a specific NNM motion of the underlying conservative, unforced system (cf. Eqs. (4)-(5)). The second step of the methodology consists in setting the input force to zero and recording the free response of the system. The oscillations remain close to the NNM of the conservative system as shown in Ref. [18]. The amplitude-(or, equivalently, energy-) dependence of the vibration frequency is then extracted from the free response time series using standard time-frequency analysis tools such as the Hilbert [32], or the Wavelet [30] transform. A difficulty associated with this second step of the methodology is the transient shaker dynamics which prevents the applied force from being actually zero and hence influences the system's response.

The control-based continuation method
CBC is a testing method inspired from numerical continuation techniques [33][34][35] that aims to follow, experimentally, the evolution of the steady-state response of a system (e.g., equilibrium, periodic solution) as parameters are varied. Contrary to numerical simulations, the states of a physical system cannot be set arbitrarily and CBC relies therefore on a controller and its control target x ⋆ as a proxy for the states. The dynamical system of interest is however the uncontrolled system, so the added control system should be non-invasive, i.e. the control input u(t) should vanish for all time. In that case, the steady-state solution x asy observed for the controlled system is also a steady-state solution of the uncontrolled system. In practice, u(t) is only approximately zero. Though the controller does not change the steady-state solution itself, it does change its linearisation thus making unstable orbits become stable if implemented correctly.
The presence of the controller is advantageous because CBC is thus robust to stability changes and bifurcations. However, investigating the stability of the underlying uncontrolled response, and hence detecting bifurcations, is more complicated. In Ref. [25], a number of measures are suggested to overcome this problem but all require turning off the control for a period of time. In many situations this is not desirable as damage could be caused to the experiment or even the experimenter. In Ref. [36], Barton identifies a multi-input multi-output (MIMO) auto-regressive model with exogenous inputs (ARX) from the experiment response perturbed around a steady-state periodic orbit. The model is then exploited to determine the so-called Floquet multipliers and conclude the periodic orbit stability.
In the context of finding steady-state behaviour, choosing x ⋆ such that u(t) ≡ 0 plays the role of the equations of motion of a model. In Refs. [23,25,26,37], the problem is iteratively solved (to experimental accuracy) using a Newton-like algorithm where derivatives are evaluated experimentally using finite differences. Starting from a given steady-state, the search for the next solution is then performed using a pseudo-arclength continuation algorithm (see, for instance, [38]). In this paper, a simplification of this procedure is used because the studied parameter, the total forcing amplitude, and the control signal have the same action on the system. This simplified method is faster than the approach reported in Ref. [23] because no derivative is required. The general aspects of this CBC approach are briefly introduced in Section 3.2.1, and the reader is referred to Ref. [27] for a more detailed description. CBC is then adapted to backbone curve tracking in Section 3.2.2.

Steady-state periodic solutions of the forced system
This section shows how CBC can extract a steady-state periodic solution of the uncontrolled system in response to an excitation f (t). In our experiment, we consider a single-point, single-harmonic forcing of arbitrary phase of the form f (t) = a cos(ωt) + b sin(ωt). The forcing amplitude r = √ a 2 + b 2 is considered as a parameter.
Consider what happens if we pick a specific harmonic forcing f ⋆ (t) defined by the pair of coefficients (a ⋆ , b ⋆ ) and an arbitrary periodic control target signal x ⋆ (t) expanded to m (finite) Fourier modes as A feedback control signal is added to the excitation signal. For simplicity, the particular case of a proportional-plus-derivative (PD) controller as later used in our experimental investigations is considered. The method works however for more general control strategies.
The total input to the system is given by where x(t) is the response of the system. Assuming that the chosen controller is stabilizing, the experiment settles into a periodic steady-state output defined as We assume the experiment has periodic input after the transients have settled. The signal u(t) is generally not equal to zero and the control system is thus invasive. The general solution to this problem is to use a root-finding algorithm to modify the control target However, realizing that f (t) and u(t) have the same action on the experiment, the applied force and control signals can be lumped together such that f (t) = f tot and u(t) ≡ 0 in the first mode. The steady-state response (8) is unchanged as the total input to the system f (t) + u(t) remains unchanged, but (8) is now a steady-state solution of the underlying uncontrolled system of interest. The amplitude of this new forcing at fundamental frequency ω equals r = √ Equations (9) and (10) are not sufficient because the total input force generally does not have the required single-harmonic form due to the presence of nonlinearities. Even if the reference signal x ⋆ (t) is harmonic, the output x asy contains higher-harmonics introduced by the nonlinearities of the experimental system. The control signal u(t) has therefore higher-harmonic Fourier coefficients given by If these coefficients are zero then the forcing f (t) + u(t) is harmonic with amplitude r = √ a 2 + b 2 such that the point (r, x asy ) is a periodic orbit of the uncontrolled system.
The requirement for the coefficients (A u 0 , A u j , B u j ) m (j=2) to be zero is a nonlinear system of 2m − 1 equations in the Fourier coefficients ] of the reference signal x ⋆ . This problem is very similar to the original one, with the notable difference that the first mode (j=1) is no longer included. This first mode usually contains all the instability present in the periodic solution and the root finding problem. Removing this coefficient from the problem allows therefore to use a more effective fixed-point iteration method where derivatives do not need to be evaluated. The k th iteration of the method reads: where In other words, the new control target coefficients X ⋆ k+1 are simply set equal to the Fourier coefficients X of the asymptotic steady-state response x asy reached under the control input defined with coefficients X ⋆ k .
In summary, the overall CBC methodology to trace out the steady-state periodic response of a system in function of the forcing amplitude r is: where (X n , X n−1 ) are the Fourier coefficients of the previous two points along the branch of periodic solutions.
2. Run the experiment with input (7) and x ⋆ defined using the Fourier coefficients X ⋆ .
3. Measure the Fourier coefficients X of the output x(t) after the transients have died out. Although not necessary, the control can be tuned appropriately such that transients die out quickly.

Check if the root-mean-square error
is smaller than the desired tolerance. If not, proceed to the fixed-point iteration. Set X ⋆ = X for all Fourier modes except the first (A ⋆ 1 and B ⋆ 1 are left unchanged) and go to step (3). 5. After fixed-point iterations, the higher-harmonic coefficients of u(t) are below the user-specified tolerance and the total input to the system can be considered as harmonic. The next point X n+1 on the branch is X n+1 = X and the forcing is given by Eqs. (9)-(10).
This method can be regarded as an amplitude sweep carried out at constant forcing frequency. This simplification of the CBC methodology can be applied to other systems than the one presented here, provided that the instability existing in the root-finding problem (A u 0 , A u j , B u j ) m (j=2) = 0 is restricted to the first Fourier mode. This has to be assessed case by case. Note that the methodology described above can be further simplified by omitting f (t). The total input to the system boils down to the sole control signal u(t), which now also plays the role of excitation. The next point X n+1 on the branch is then given by X n+1 = X and (a n+1 , b n+1 ) = (A u 1 , B u 1 ).

Tracking the backbone curve
To track the backbone curve, the oscillation frequency must be updated to maintain the phase quadrature between excitation and response signals. This problem can be formalized as solving the scalar equation where φ out (ω) and φ in (ω) are the phase of the fundamental Fourier modes of the response, x(t), and total input,f (t) + u(t), respectively. The evaluation of q(ω) is performed after the fixed-point iterations and replaces step 5 in the methodology of Section 3.2.1. If q(ω) is below a used-defined tolerance, the point X is recorded as the next point on the backbone curve (with X n+1 = X and forcing given by Eqs. (9)- (10)). If q(ω) is above the prescribed tolerance, the frequency is updated and we go to step 3. Equation (16) can be solved using a Newton-Raphson procedure. A simpler bisection method is also effective.
Equation (16) only accounts for the fundamental harmonic component because the excitation signal considered in the present study is harmonic. If a richer, multi-harmonic, excitation signal was to be considered, the higher-harmonic coefficients would have to be updated to satisfy the quadrature criterion. These extra unknowns could be balanced by extending the quadrature condition (16) to include the phase between the higher harmonics. As these higher-harmonic coefficients would now be considered as part of the excitation, they would not be included in the error (15) and would no longer be updated in the fixed point iteration procedure.

Set-up description
The SDOF oscillator shown in Figure 1 is investigated under base excitation. The oscillator can slide with linear ball bearings along two guide rails. Two linear springs working in the transverse direction are used to attach the oscillator to a metallic supporting structure, which itself is mounted on a perspex base-plate that plays the role of uni-axial shaking  the shaking table oscillations about a central position. The main purpose of this PID controller was to avoid the table drifting -the fine tuning of the control gains was not necessary for CBC. To control the periodic orbits, a PD controller based on the mass displacement was implemented and added to the PID base-displacement control. To determine the controller gains, the system was forced at 3 Hz -a frequency for which two stable steady-state periodic responses coexist, one at high amplitude and one at low amplitude. When the system was at low amplitude, the target coefficients of the controller were set to be equal to the Fourier coefficients of the upper stable periodic response. The first set of gains that was able to realize the transition between the lower and upper responses was selected. The proportional and derivative gains were k p = 0.002 and k d = 0.0005, respectively. The derivative term was computed using a two-point finitedifference approximation after applying a fourth-order IIR Butterworth filter (-3dB cutoff at 25 Hz) to the error signal. The filter is purely for the purpose of control; all other calculations use unfiltered data.
Displacement and force signals are decomposed in real-time into their first seven Fourier modes from the unfiltered data. An effective recursive estimator is used to minimize sampling and noise effects caused by the forcing period not being an integer multiple of the sampling period. The recursive estimator for the k th Fourier coefficients is For convenience, the Fourier decomposition was performed in real time but it does not need to be. The time series can be recorded and post-processed off-line using FFT or least-square techniques.

Experimental results and comparison with resonant decay
The CBC algorithm presented in Section 3.2 was applied to the SDOF shown in Figure 1. The objective function (16) was defined using the phase of the base displacement and of the mass displacement relative to its support, thus focusing on the SDOF oscillator and discarding any influence of the shaking table. The tolerance on (16) was 5 × 10 −3 rad., and the minimum frequency step was 10 −4 Hz. Recorded Fourier coefficients were averaged over 5 samples, and each modification of the target coefficients X ⋆ was followed by a waiting period of 4 seconds to let the transients die out. Starting from rest, the target coefficient A ⋆ 1 was initially increased by 7 mm to overcome friction present in the set-up. The backbone curve was then discretized using a constant amplitude step h = 2.5 mm. The frequency was adapted as described in Section 3.2.2.
The measured backbone curve is presented in black in Figure 2. The total experimental time required to generate the curve was 100 minutes for a total of 41 points. The relative displacement amplitude of the response was computed using the Euclidean norm of the Fourier coefficients. For low amplitudes (below 10 mm), the system presents a softening characteristic which could possibly arise from spring inertial effects [39]. Similar experimental observations were reported in Ref. [20]. At around 10 mm the backbone curve presents a turning point, above which the system presents a hardening characteristic with a resonance frequency increasing from 2.5 Hz to 3.2 Hz in the [10 -70] mm displacement range. At high-amplitude, the fundamental Fourier coefficient is still the main contributor to the response -the relative importance of the largest higher harmonics compared to the fundamental does not exceed 1.5%. This observation confirms our assumption that, in the present case, a single-harmonic excitation was sufficient to accurately reach quadrature and isolate the NNM motion.
The identification of a linear regime of motion was not possible because friction and the very low pretension applied to the springs result in a system where nonlinear effects dominate at all amplitudes.
The backbone curve was also extracted using resonant decay. The free response of the SDOF oscillator was obtained after applying an initial displacement of about 70 mm to the mass. The shaking table was clamped to avoid its dynamics influencing the experimental results. Before extracting the backbone curve, the raw time signal was filtered using a fourth-order IIR Butterworth filter (cutoff frequency at 25 Hz) in order to have a similar filtering action than the seven-mode Fourier decomposition. Whilst there exist many methods for analysing transient signals [30,32], the filtered time series were analysed using zero-crossing and amplitude interpolation as described in Ref. [20]. As such, contrary to CBC, backbone discretization is governed by the damping characteristics of the system. The extracted backbone is reported in blue in Figure 2. Similarly to the CBC results, the resonant-decay backbone presents softening and hardening regions. CBC and free decay agree very well for amplitudes larger than 20 mm. At 20 mm, the curves depart from each other and cease to overlay closely for smaller amplitudes.
The variability of our experimental results is investigated in Figures 3(a) and 3(b) where the results obtained for three different (consecutive) CBC and resonant-decay runs are superimposed, respectively. Overall, the repeatability of both methods is excellent and, arguably, slightly better using CBC than using resonant decay. Interestingly, the backbone curve measured with CBC consistently shows a small amplitude drop at around 20 mm that is not reproduced by resonant decay. However, below 20 mm, the curves appear noisier and discrepancies between curves are noticeable for both methods. This region was found to be very sensitive to experimental conditions such as the day of test and the number of tests performed. We believe that this variability is attributable to friction and complex damping mechanisms present in the mass' ball bearings.
An advantage of the CBC method over resonant decay is to offer means to verify and validate the quality of the experimental results. Besides the convergence of the objective function (16), the assumption of single-harmonic base excitation can be assessed. The invasive character of the (PD) controller can also be measured in order to estimate its influence on the dynamic features of the uncontrolled system. These verifications are performed in Figure 4 where the root-mean-square (RMS) value of three different time series are shown in terms of the Fourier components of the base displacement. The first time series (--) represent the first (fundamental) harmonic component of the base displacement. The second time series (-•-) is made of the sum of the higher-order harmonics, and the third one (-♦-) accounts for the non-harmonic content present in the unfiltered data and that is not decomposed in the seven first Fourier modes. From Figure 4, we can see that the base displacement is essentially mono-harmonic. The relative importance of the higher-harmonics is very low compared to the fundamental component (maximum 2.7%), thus confirming the successful convergence of the fixed-point algorithm used to minimize Equation (15). The non-harmonic part of the signal should ideally be as small as possible. Its relative importance with respect to the fundamental component is about 3%, which is  overall comparable to the importance of the higher-harmonic content. However, several points exceed of the general trend, with the maximum reaching a RMS value equivalent to 17% of the fundamental component. These points are symptomatic of inaccurate measurements because the Fourier decomposition is less representative of the measured signal and the assumption of harmonic forcing no longer holds. These points might be rejected if deemed necessary. Interestingly, the largest value of the non-harmonic content is reached at a response amplitude of 10 mm, which precisely corresponds to the point where the backbone curve changes from softening to hardening. It is also the region where the results present the largest variability (see Figure 3(a)). Finally, the sum of blue and red curves is a direct measure of the invasiveness of the control. This sum is not represented in Figure 4 because it is very close to the blue curve. If the implemented PD control was truly noninvasive, this sum would be identically zero. The predominant source of error in our experiment is the noise amplification through the derivative term in the control signal. Another contribution to the error comes from the important friction effects present in the shaking table, and for which no particular effort in tuning the base displacement controller was made.  Figure 5 shows in a three-dimensional plot (forcing frequency, base amplitude, response amplitude) the data points collected during 17 amplitude sweeps. The sweeps are 0.5 Hz apart and were performed following the procedure described in Section 3.2.1. Extracting a point of Figure 5 is approximately fifteen times faster than extracting a point on the backbone. This important additional cost for backbone continuation is explained by the necessity to solve an extra equation for the phase quadrature and to wait for the transient to die out after each frequency modification. The total experimental time required to generate Figure 5 with 700 data points is about two hours, which is more than extracting the backbone directly.

Forced response of the SDOF oscillator
A Gaussian process (GP) regression [40] was used for constructing a smooth surface out of the experimental data points (grey surfaces in Fig. 5). GP hyper-parameters were calculated by maximizing the marginal likelihood of the hyper-parameters. The obtained interpolant provides a geometrical model of the response surface, thus playing the role of surrogate model for the SDOF dynamical system. Numerical continuation can be used for extracting geometrical features of the forced-response surface such as the curve of saddle-node bifurcations (-). The dark-grey region defined by this curve is a region where periodic solutions of the uncontrolled system are unstable. It would typically be impossible to observe the data points of this region without control.
In Figure 6, the surrogate model is further exploited to trace out the forced response of the oscillator for different constant amplitudes of the base displacement (1.0, 1.5, 2.5, and 3.5 mm). At high amplitudes, the obtained curves are reminiscent of the resonance curves obtained for an idealized Duffing oscillator. The softening and hardening character observed on the measured backbone (also reported in orange) is well reproduced. The proximity of the backbone curve with the saddle-node bifurcation (black dot in Figure 5) is usually challenging for methods like resonant decay because a too-fast variation of the forcing frequency during the force appropriation phase can bring the system in a region where the only stable steady-state solution is at low amplitude, forcing the system to jump down to this solution away from the backbone curve. Moreover, the set of initial conditions (i.e. the basin of attraction) leading to high-amplitude steady-state responses shrinks when the saddle-node bifurcation is approached. Transients can therefore perturb the steady dynamics of the system such that it jumps before the actual bifurcation is reached. When excitation is applied, precisely reaching the phase quadrature condition can be a long and tedious task. CBC does not face such difficulties because the embedded control system guarantees to stay in the vicinity of the prescribed control target.

Model derivation
In this section, a numerical model of the SDOF oscillator presented in Figure 1 is studied in order to further compare and discuss the differences that exist between the CBC and resonant decay methods. We note that no quantitative correspondence between the model and the physical set-up was attempted, rather the emphasis is on reproducing the behaviour qualitatively. The dynamics of the oscillator is modelled as where m = 0.33 Kg is the mass of the oscillator. The linear viscous damping coefficient c was estimated at 0.6 N.s/m using an exponential fitting of the free response envelope. The Coulomb friction force coefficient c nl was estimated around 0.05 following the method described in Ref. [41]. Both viscous damping and friction coefficients will be varied for the purpose of our numerical investigations. The functional form of the nonlinear restoring force was deduced from the shape of the backbone curve observed in Section 4.2. The procedure described in Ref. [14] was followed to estimate the restoring force coefficients (k 1 − k 4 ). Specifically, the response of the system and the applied force are assumed to be of the form x(t) = X cos(ωt) and f (t) = F cos(ωt + φ), with φ an arbitrary phase. These approximations are introduced in Equation (18). The sign(·) terms represent square waves of amplitude 1 and period T = 2π/ω, whose Fourier series expansions are given by sign(sin(ωt)) = 4 π n odd 1 n sin(nωt), sign(cos(ωt)) = 4 π n odd (−1) (n−1)/2 n cos(nωt).
Taking into account Equations (19) and (20), the harmonic coefficients in the equation of motion can be balanced to give, at first order, the following two equations −mω 2 X + k 1 X + 8 3π −cωX − 4c nl π = F sin φ.
In complete analogy with the results presented in Ref. [19] and briefly reviewed in Section 2, Equation (22) shows that if the excitation is in quadrature with the response of the system, i.e. if φ = π/2, it exactly cancels out damping effect. This result not only applies to systems with linear viscous damping, but also to systems including Coulomb friction and any other form of nonlinear damping that can be described using an odd function of the velocity. The left-hand side of Equation (21) is also identically zero and x(t) satisfies the equations of motion of the underlying conservative system, which by definition is a NNM motion.
From Equation (21), an equivalent amplitude-dependent stiffness K eq can be derived as such that the system has a natural frequency The backbone curve obtained in Section 4.2 provides at each point a pair (ω, X), from which the k i 's coefficients can be inferred. The estimated values are k 1 = 79 N/m, k 2 = −165 N/m 2 , k 3 = 3.6 × 10 4 N/m 3 , k 4 = 3 × 10 5 N/m 4 .

Effect of linear viscous damping
The backbone curve and the forced response of the SDOF oscillator were computed using an algorithm combining shooting and pseudo-arclength continuation [42]. The free response of the system was computed using a fifth-order Runge-Kutta algorithm with adaptive time stepping, considering zero applied force and an initial displacement x(0) = 80 mm. The frequency-amplitude dependence of the oscillations was extracted from the noise-free time series using the same method as in Section 4. The forced response of the system is also displayed (-), dashed lines indicate unstable periodic solutions. Overall the backbone curve from continuation is well captured by the resonant-decay method. However small errors are present for amplitudes greater than 50 mm in the larger-damping case. Larger-damping values also decrease the amplitude-frequency discretization of the backbone curve, as indicated by the number of markers along the curve. The CBC approach developed in this paper does not suffer from this limitation because the backbone discretization is governed by the user-defined parameter h.
Although no quantitative comparisons between the experimental and numerical models are attempted, it is interesting to note that the softening effect present in the model is far less important than in the experiment, arguably indicating that the physics underlying the softening effect in the experimental set-up is not properly captured by our model.

Effect of Coulomb friction
The effect of friction is illustrated in Figure 8. Linear viscous damping was set to a low value of 0.06 N.s/m in order to minimize its effect on the response. Contrary to viscous damping, Coulomb friction directly impacts both the quality and discretization of the identified backbone. The softening region is no longer captured for friction coefficients greater than c nl = 0.02, and the oscillator reaches a stick condition when vibration amplitudes are low enough (around 2.4 Hz). These results contrast with the result that would have been obtained with the CBC approach. The backbone curve remains accurately measured as the phase quadrature condition still holds (cf. Equations (21)(22), thus showing that the CBC approach works on systems with specific forms of nonlinear damping such as, for instance, Coulomb friction. Moreover, regions of inaccurate results as observed in Figure 8 are typically difficult to identify based on the sole free response

Conclusions
In this paper, a novel approach for backbone curve identification was developed based on control-based continuation. A SDOF nonlinear oscillator was considered as demonstrator, and the developed method showed a number of advantages over the resonant-decay approach. The control system embedded in CBC overcomes the potential issues arising from the presence of bifurcations and stability changes in the system's response. In this way, the proposed method can be considered as robust with respect to the system's dynamics, and can easily be made automatic. CBC also allows to verify and quantify the quality of experimental results. With resonant decay, it is difficult to remove the influence of the excitation system from the free response dynamics -unless excitation is replaced by initial conditions in displacement (as performed in this paper). CBC naturally overcomes this issue because it focuses on the relative dynamics between the excitation and investigated systems, whilst it guarantees a specific excitation signal form and the non-invasiveness of the control system. Another important facet of the proposed approach is to directly mea-sure the backbone curve in the forced experiment. The method remains valid as long as the phase quadrature condition holds, which is particularly interesting for systems with Coulomb friction as demonstrated in Section 5 using a numerical model. In principle, tolerances on Equations (15) and (16) can be changed to extract the backbone curve with arbitrary accuracy (within experimental limitations). Unlike resonant decay, backbone curve discretization is not governed by the system's damping characteristics but follows a user-defined parameter. Finally, one should note that the method presented here does not only extract the backbone curve but also gives the periodic solutions defining the backbone curve, i.e. the NNMs of the underlying conservative system, which offers new comparison means between theoretical and experimental results.
All these advantages come however at the cost of an increased testing time compared with resonant-decay. Although the proposed methodology is a priori applicable to MDOF systems, it was only tested on a SDOF system for which a single-point, single-harmonic, excitation was sufficient to properly isolate the backbone curve. Future investigations will address more complex systems with closely-spaced (or interacting) modes, for which multi-point multi-harmonic excitations would be necessary.

Data statement
All the experimental data used in this paper has been deposited into the University of Bristol Research Data Repository and is publicly available for download [43].