Large basins of attraction for control-based continuation of unstable periodic states

Numerical continuation tools are nowadays standard to analyse nonlinear dynamical systems by numerical means. These powerful methods are unfortunately not available in real experiments without having access to an accurate mathematical model. Implementing such a concept in real world experiments using control and data processing to track unstable states and their bifurcations, requires robust control techniques with large basins and good global properties. Here we propose design principles for control techniques for periodic states which lead to large basins and which are robust, without the need to have access to a detailed mathematical model. Our analytic considerations for the control design will be based on weakly nonlinear analysis of periodically driven oscillator systems. We then demonstrate by numerical means that in strong


Introduction and experimental context
Dynamical systems theory, in particular the investigation of instabilities and of chaotic motion in nonlinear systems is one of the key themes in theoretical and experimental sciences of the last decades [1,2,3].Numerous tools have been developed to study the complex motion in such systems.From a theoretical perspective continuation of stable as well as unstable solutions and of bifurcations play a prominent role [4] since global structures in phase space and a skeleton of unstable states for the dynamics can be discovered which gives insight into details of the dynamics [5].Equation free analysis has been proposed as a practical tool to derive effective equations of motion from microscopic models by numerical means [6].These effective models can then be studied by analytic and numerical bifurcation analysis to understand the lowdimensional nonlinear structures and pattern formation in complex real world systems.
The idea of a model free approach has taken this concept a step further by proposing an approach where one avoids the intermediate derivation of effective equations of motion, making the idea directly applicable to experiments.In this context the required continuation of stable or unstable dynamical states relies on non-invasive control methods which are directly applicable to real world experiments.Non-invasive control, sometimes called orbit control in the engineering context, has the crucial property that the control forces ultimately tend to zero, so that the stabilised state is a genuine unstable orbit of the original system without control (see as well [7], where this idea has been popularised within the physics community).In addition to numerical implementations of an equation-free approach to obtain bifurcation diagrams at a macroscopic scale, see for instance [8], a successful implementation of this concept has been already demonstrated in mechanical and electro-mechanical hybrid experiments [9,10,11], in electrochemical setups [12], or even in experimental studies of pedestrian flows [13].Key to these experimental implementations of continuation techniques is the availability of suitable control methods which can deal with quite diverse experimental conditions.Control problems have of course a long standing tradition in engineering.They took a large boost during and after world war two, where linear control theory was formalised in a systematic manner.Subsequently these ideas have been extended to nonlinear systems, see e.g.[14,15].Putting the emphasis on non-invasive methods the relevance of control techniques for the purpose of system analysis has been rediscovered in the context of chaotic dynamics [7].Unlike in engineering, control is considered here merely as a kind of spectroscopic tool to identify structures in the phase space of the system.
To illustrate some of the challenges one faces when one aims at implementing control based experimental continuation of unstable states we refer to atomic force microscopy as a paradigm.
A comprehensive theoretical study of control-based continuation in a model of atomic force microscopy can be found in [16].In experimental terms atomic force microscopy is a quite versatile method to inspect surfaces and attached structures or objects thereon, down to the nanoscale [17,18].During the last decades scanning platforms have been substantially augmented, so that frame sizes of 80µm down to 1nm are possible without any drawback in lateral and vertical linearity.During topography acquisition under ambient conditions the relative amplitude drop can serve as setpoint variable while the control variable is the piezo actuator height.To a good degree, thereby the tip-sample separation is kept constant during lateral scanning and a topography is acquired as a discrete data set.Typically, physical lateral resolution is at about 1nm while vertical resolution is at about 1 Å.It is quite common to operate atomic force microscopes in the dynamic mode where the sensor element, a bending micro-cantilever with a tip, is driven into oscillation via a dither or shaker piezo.The nonlinear features in the interaction with the surface turn atomic force microscopy into a complex nonlinear dynamical systems which shows bistability, a variety of bifurcations, and chaotic motion [19,20,21,22].In dynamic force microscopy stable low and high amplitude branches with an unstable branch in between do coexist.Jumping between the two stable branches during acquisition of topography is a common distortion, particularly with molecular species, requiring manual intervention (see figure 1 for experimental images).Especially beyond the resonance frequency of the cantilever a rather broad range of tip sample separations is prone to such imaging instabilities.By applying ideas from time-delayed feedback control to atomic force microscopy noise reduction in the imaging process has been reported [23].
We have mentioned atomic force microscopy just to motivate our plain theoretical considerations.Control challenges which may occur in such an experimental context, such as the lack of a proper mathematical model, fast time scales like those encountered in nanosystems which prevent extensive online data processing, parameter drifts and non-stationary behaviour which preclude a preliminary data based modelling, or the dynamical impact of noise are the topic of our interest.Above all we want to design a control scheme which is robust and comes with a large basin of attraction, and which finally can be set up if just plain measurements of the system are available.In a previous theoretical study [24] we have outlined a scheme to deal in principle with these issues by linear control schemes applied to stroboscopic maps.However, such linear schemes have severe limitations when it comes to global properties of the control in nonlinear setups.In our work we focus on the design of non-invasive control schemes with large basins of attraction.We will demonstrate their success by numerical simulations.Actual implementations in experiments will be addressed elsewhere.
The design of control schemes with good global properties is of course a standard theme in engineering.For instance, the problem of globally stable control design has been solved by exact state-space linearisation for single-input systems, see [14].Having said that, the implementation of such ideas requires some knowledge about the underlying dynamics.In very basic terms the issue of globally stable control has been revisited in the context of control of chaos [25] and has been popularised beyond the remit of engineering problems.To keep our presentation selfcontained section 2 will give a brief sketch of the basic ideas and of the related challenges in a theoretical setting.Design of control schemes with good global properties requires access to some properties of the dynamics.Since we have ultimately applications for atomic force microscopy in mind we focus here on general nonlinear driven oscillators and the corresponding nonlinear resonance behaviour.We will outline in section 3 how weakly nonlinear perturbation expansions will give us analytic access to the dynamics, and in particular to the stroboscopic map of the equations of motion.Based on these analytic estimates we will show in section 4 how to design a globally stable control scheme to stabilise periodic orbits in a non-invasive way.In particular, our design utilises the phase of the driving field as a key component of a globally stable control scheme.While our approach has been based on weakly nonlinear analysis and on the analytic expression of the stroboscopic map, we will show by numerical means that our control scheme also works well beyond the perturbative regime and can thus be applied to generic nonlinear oscillator systems.We demonstrate in section 5 that the control scheme can be used for a data based tracking of unstable orbits and for the generation of a complete bifurcation scenario.While the design of the control scheme has been based on a perturbative treatment of model equations, we will show in section 6 that all the elements of the control scheme can be obtained from data, in particular from scanning a bistable nonlinear resonance curve, as long as higher order harmonic components do not dominate the dynamics of the system.Above all, the control scheme can be implemented without any a priori access to a mathematical model.
Finally, we briefly discuss in the conclusion limitations and merits of the proposed approach, in particular in the context of atomic force microscopy and related experimental setups.

Globally stable control in a nutshell
As we are aiming at controlling periodic states of a dynamical system it seems promising to focus on Poincare maps or stroboscopic maps since periodic states will become fixed points of the time discrete dynamical systems.To illustrate the basic idea how to design a globally stable control scheme consider a time-discrete dynamical system given by a one-dimensional map f µ where the right hand side depends on a parameter µ which will serve as control input.For simplicity and for the purpose of illustration we assume here that the map depends on the parameter µ in an additive way.Such an assumption is by no means essential for the subsequent considerations.We aim at controlling fixed points of the dynamical system, eq. ( 1), where the fixed point manifold, that means the fixed point x * in dependence on the parameter µ, is determined by In order to stabilise such a fixed point we make the parameter µ a dynamical time-dependent quantity.In its simplest instalment the time-dependence of µ is just given by a static relation with the state variable, say where h specifies the control law, i.e., the dependence of the parameter µ n on the state of the system x n .Then the closed loop dynamics reads The actual fixed point to be stabilised is determined by eq. ( 2) together with In geometric terms eqs.( 2) and ( 5) amount to the intersection of two manifolds.The static offset µ R which has been included in the control law, eq. ( 3), serves as an external parameter by which the fixed point to be controlled, x * , and the corresponding actual parameter value µ * can be selected.
We aim to chose the control feedback h(x n ) in such a way so that the fixed point x * becomes a globally stable fixed point of the closed loop dynamics, eq. ( 4).The obvious choice is of course h(x) = −g(x) since the fixed point of the closed loop dynamics becomes superstable and convergence happens in one iteration step.The offset µ R of the control loop becomes in fact the fixed point value to be stabilised and the corresponding parameter value is given by eq. ( 5).
Global convergence can be obtained as well under less stringent conditions.For instance it is sufficient that the closed loop dynamics, eq. ( 4), gives rise to a contraction map.Hence a rough estimate of the full internal dynamics may be sufficient to design a successful control scheme.
As a simple illustration consider the logistic map on the interval [0, 1] which reads To stabilise the fixed point we employ the control scheme µ = µ n with the choice Then the closed loop dynamics yields a contraction on the interval [0, 1] and the dynamics converges globally to a unique fixed point x * .The actual value of the fixed point is determined by the control gain µ R , and the corresponding actual parameter value is given by µ * = µ R h(x * ), see eq. ( 7). Figure 2 shows time traces of successful control of the unstable fixed point for some typical parameter values.
It is worth to mention, that by construction, the control scheme outlined above is non-invasive, even though a control force does not seem to tend to zero.By making the system parameter µ a dynamical quantity µ n , which is adjusted by an instantaneous feedback law such as eq.( 3) or eq.( 7), successful control is signalled by the dynamical parameter tending to a limit value µ ∞ (see figure 2).In such a case the resulting dynamics for x, in our case a fixed point (see figure 2), coincides by construction with an orbit of the system without control at fixed parameter value µ ∞ .Hence, we have a non-invasive control scheme which stabilises a proper unstable orbit of the original equations of motion.In fact, the fixed point which is stabilised in figure 2 (red symbols), coincides with the unstable fixed point of the logistic map for the limit parameter value µ ∞ (figure 2, cyan symbols).We also note that the reference value µ R in the feedback law tunes the limit value µ ∞ and the control target, see eq.( 5), but is by no means per se the target value of the control scheme.To design a globally stable control law one needs the full knowledge or at least a sufficiently accurate estimate of the underlying equations of motion.That is of course far from surprising as a globally stable solution requires the control law to remove all potentially occurring basin boundaries of the stabilised state.Such boundaries are normally caused by unstable saddles and their stable manifolds which, as said, have to be removed from the equations of motion.
The reasoning summarised in this section is far from novel.Such considerations are at the heart of any sensible control design and they are normally covered in any basic course in control engineering.The main idea is also well established for, say, almost a century and simplified versions have also popped up in the physics context, see e.g.[25].Nevertheless we felt it useful to recall these basic considerations as we will use such reasoning to design globally stable control of periodic states in oscillator systems in the next sections.
3 Stroboscopic map for weakly nonlinear oscillators We will mainly deal with control of unstable periodic solutions in two-dimensional driven oscillator systems, which are described by the set of differential equations Here γ denotes the viscous damping constant, U (x) the anharmonic part of the potential, and h (s/c) the amplitudes of a harmonic periodic driving force.As pointed out in the previous section we will require some degree of access to the full dynamics, hence we adopt a system with a small parameter ε which can be employed to perform an analytic perturbation expansion.
The purpose of the model, eq.( 9), is twofold.For ε = 1 the model constitutes a general nonlinear driven oscillator which will be used to test our control design (see sections 5 and 6 ), even to the level where we base control only on a recorded time series of the model.On the other hand we will use the regime of small ε in eq.( 9) to develop, inspired by analytic perturbation expansion, a globally stable control design in section 4. In particular, we will show that access to the amplitude and the phase of the driving field is sufficient to set up the control scheme.
While this design will work by construction in the perturbative regime we will also demonstrate that good global properties persist for typical oscillators.
As far as dynamical phenomena of eq.( 9) are concerned we can simply resort to numerical computations of stroboscopic maps, that means computing and study the properties of the map (x n , v n ) → (x n+1 , v n+1 ) by numerical means.While in our analytical studies we will keep a general potential we use for numerical purposes a simple Duffing oscillator with potential given by and we chose in numerical simulations a standard set of parameter values A dominant dynamical signature of such oscillator systems is a nonlinear resonance which occurs even beyond the perturbative regime of small ε values, where a bistability between a large and a small amplitude branch appears.That can be vividly illustrated by a numerical computation of time traces of the stroboscopic map which tend toward fixed points, and where the bistability occurs when the amplitude of the driving field is changed in a quasi-stationary manner, see figure 3.This bistability will be at our centre of interest as we aim to stabilise the unstable branch which separates the two stable states in the bistable region, with a view towards performing  One of the two amplitudes h (c) or h (s) in eq. ( 9) seem to be redundant because Here θ denotes the phase of the driving field relative to the cross section of the stroboscopic map, i.e., relative to the times n2π/ω of the observation.The phase θ can be easily eliminated by choosing an appropriate cross section for the stroboscopic map, i.e., by shifting the time of observation by a constant amount.However, we will see soon that having access to the phase of the driving field will be a key element to construct a globally stable control scheme.In fact using the phase of the driving field for the purpose of control is by no means novel, see e.g.[26].
Such feedback has been proposed in experimental contexts, for instance, for noise reduction at the microscale [27,28].However, neither non-invasive control, nor control with good global properties has been at the centre of interest in these studies.
As we have seen in the previous section a successful design of a suitable control feedback requires some basic understanding of the dynamics of the underlying system.That means some analytic access to the stroboscopic map is beneficial.While a closed analytic expression for the stroboscopic map is never available in non-trivial situations, we can get some insight if we constrain to the perturbative regime of small ε values.If x(t, x n , v n , ε) and v(t, x n , v n , ε) denote the time dependent solution of eq. ( 9) with initial condition (x n , v n ) the exact stroboscopic map is given by ( straightforward series expansion of the solution of eq. ( 9) in terms of ε we are able to derive an analytic approximation of the stroboscopic map.At lowest order O(ε 0 ), eq. ( 9) reduces to a harmonic oscillator and the solution with initial condition (x n , v n ) reads At first order O(ε) we obtain the linear inhomogeneous system ẋ(1 with initial condition x (1) (0) = 0, v (1) (0) = 0.The solution can be easily computed and we finally obtain for the solution after a period 2π/ω, using the expression eq. ( 13) x (1) To evaluate the remaining integrals we introduce the amplitude r n and the phase φ n of the solution eq. ( 13) by Then we have A simple substitution shows that the last integral is an odd function of r n , that means the integral can be written as r n times an even function.Therefore we can define where the quantity w can be viewed as an effective suitably averaged force.Then eq. ( 17) reads if we take the abbreviations ( 16) and ( 18) into account.The real and imaginary parts of this expression yield the remaining integrals contained in eq. ( 15).Eqs. ( 13) and ( 15) result in the stroboscopic map at first order in the expansion parameter This result of the perturbation expansion is valid for any type of potential U (x), as long as the integral in eq.( 18) does not vanish and defines a meaningful effective force w.In particular, the result, eq.( 20), does not rely on any symmetry properties of the potential1 .For the particular case of the Duffing oscillator, eq. ( 10), eq. ( 18) readily gives the explicit expression It is rather straightforward to compute the fixed points of the map, eq. ( 20) in closed analytic form.In fact, for the fixed point (x * , v * ) eq. ( 20) results in the linear system of equations Solving this system for x * and v * we get and using the definition of the stationary amplitude, r 2 * = x 2 * + v 2 * /ω 2 (see eq. ( 16)) we finally obtain an implicit equation for r * To demonstrate the accuracy of the perturbation expansion we compare bifurcation diagrams obtained directly from the numerical integration of the equations of motion with those computed from the first order perturbation result, eq. ( 20) or eq.( 24). Figure 4 shows the stationary amplitude in dependence of the driving amplitude for quasi-stationary parameter upsweeps and downsweeps.The results obtained from the analytic first order expression are surprisingly accurate when compared with numerical simulations for small values of ε.Actually deviations turn out to be so small that at the scale used in figure 4 no difference between the simulation of the differential equation, the iteration of the analytic map eqs.( 20) and ( 21), and the analytic expression eq. ( 24) is discernible.For large values of ε the first order truncation in eq. ( 20) fails when iterations of this map are considered, as the sequence of iterates tends to diverge.Nevertheless, the bifurcation diagram of stationary states for large values of ε, see figure 3, is still in qualitative and to some extent even in quantitative agreement with the data obtained analytically for small ε, cf.figures 3 and 4. While time scales and transients for small and large values of ε vastly differ there seems to be little change in the stationary states, i.e., in the location of fixed points.This coincidence is surprising but not totally unexpected as our perturbation scheme is essentially equivalent to the averaging principle or the principle of harmonic balance.While those approaches are formally first order perturbation schemes they sometimes perform well even beyond the perturbative regime as these approaches can also be viewed as a non-systematic mean field type expansion.

Design of globally stable control
Our main goal is to construct a control feedback so that stabilisation of an unstable state occurs for a large set of initial conditions, if possible even globally in the entire phase space.Following the basic reasoning outlined in section 2 we would need some information about the underlying stroboscopic map.At least in a perturbative regime of small values of ε such information is available as demonstrated in section 3. Hence we will base our control design on the expression eq. ( 20) for the stroboscopic map.Thanks to two driving field amplitudes, i.e. thanks to the fact that we have taken explicitly the phase of the driving field into account, a field amplitude occurs as an additive part in each of the components of the stroboscopic map, and we can use these amplitudes for control purposes.Following the reasoning outlined in section 2 it looks tempting to implement a control feedback which removes the nonlinear part, so that the remaining dissipative linear contribution ensures global convergence towards a stabilised fixed point.Such a reasoning leads us to the design where r 2 n = x 2 n + v 2 n /ω 2 , see eq. ( 16).The reasoning which lead us to the design, eq. ( 25), follows the idea of feedback linearisation which is well established in the engineering context, see e g. [30].However, we want to emphasise that the control scheme can be implemented easily by experimentally accessible parameters such as the amplitude and the phase of the driving field.
Evaluating the fixed point condition under control (see eq. ( 20) with h (c/s) replaced by h so that the two offset values h ∞ which correspond to the stabilised fixed point (x * , v * ), as already illustrated in figure 2 for our toy model.Having access to the phase of the driving field, i.e., having thereby access to both terms of the driving field has turned out to be crucial for our design.
While the setup defined in eq. ( 25) will by definition work for the stroboscopic map in first order of ε we still need to confirm whether the scheme has good convergence properties when applied to the full equations of motion.For that purpose let us first consider numerical simulations for small values of ε where first order perturbation theory has turned out to be quite accurate even at a quantitative level, see section 3. Implementing the control means that at the beginning of each period of the drive we readjust the amplitudes of the driving field using eq.( 25), i.e., using the current values of x n = x(n2π/ω) and v n = v(n2π/ω).We could use in principle any values for the offsets h R .However we are aiming at stabilising unstable states in the bistable regime and for driving fields which have no non-trivial phase (i.e.driving fields with ultimate value h (s) ∞ = 0).In practice one could tune such offsets when stabilisation has been achieved, but there is also a way to determine suitable offsets a priori, which we will address in the next section.For the moment we just make up "suitable" values out of thin air.
Figure 5 shows time traces obtained for the driven Duffing oscillator with the control scheme along the lines of eq. ( 25).The time traces of the phase space coordinates prove successful stabilisation with a limiting value of r 2 ∞ = 0.448 . .., while the time traces of the control forces h Hence, even if the basin for successful control is finite, its size is so large that the scheme can be considered as sufficiently robust.It is not so surprising that the control design works quite successfully for small values of ε since the analysis of the previous section has shown that this case is quite well covered by the perturbative analysis.As a test for our approach we address parameter settings beyond the perturbative regime and apply the control defined by eq. ( 25) to the Duffing oscillator, eqs.( 9) and (10), for larger values of ε, say ε = 1.Time traces of the oscillator subjected to control with the parameter setup used in figure 3 are shown in figure 6.The performance in the strongly nonlinear regime, ε = 1, is in fact comparable to the perturbative regime, cf.figure 5, even though the transients are now much shorter thanks to the larger dissipation εγ.For large times the amplitude tends towards the limit r 2 ∞ = 0.411 . . .and the amplitude of the driving field has the limiting value h (c) ∞ = 0.140 . .., so that the asymptotic state is in fact on the unstable branch of the strongly nonlinear Duffing oscillator, see figure 3.In addition there is as well a small component h (s) ∞ = −0.0066 . . .giving rise to a non-vanishing phase in the drive, see eq. ( 12), which could be removed by a slight adjustment of the offsets h R .Overall our control design performs quite well even beyond the perturbative regime.To judge the overall performance we investigate by numerical means the basin of the control with parameter settings used in figure 6, i.e., h (c) R = 0.03515 and h (s) R = 0.05713.We are no longer within the perturbative regime and the basin is finite, see figure 7.However, the basin is still quite large and covers all the states which occur in plain parameter sweeps, cf.figure 3. Hence the control design can be considered as a success even from an experimental point of view.The boundary of the basin resembles the fractal structure caused by homoclinic tangles.
Hence, the basin boundary in our case may be caused by the stable manifold of a saddle.

Control parameter setting and tracking
The control offsets h R determine the final amplitudes of the driving field and the fixed point to be stabilised.Since the control design has good global properties one could start with any values for the offsets.After successful control a continuation of unstable states could be done, as usual, by small quasi-stationary changes of the offsets.x(0) Figure 7: Basin for stroboscopic control, eq. ( 25), with control parameter setting h (c) R = 0.03515 and h (s) R = 0.05713 (see figure 6) applied to the driven Duffing oscillator, eqs.( 9) and ( 10) with parameter values as in eq. ( 11) and ε = 1.0.The filled circle (cyan) indicates the stabilised orbit.
If one employs the knowledge from the weakly nonlinear analysis one can do better and determine a priori estimates for suitable offsets to stabilise an unstable state with an approximate amplitude r 2 * .Given a value for r 2 * , eq. ( 24) tells us within the limits of first order perturbation theory that the corresponding amplitudes of the driving field read when we impose a constraint on the phase of the driving field.Using these values to estimate the actual fixed point coordinates via eq.( 23) we have Finally eq. ( 26) yields the estimates for suitable control offsets as We have in fact used eq.( 29) to determine offsets for the control in the previous section (with r * = 0.4).The actual fixed point controlled, see figure 5, then differs slightly from the estimate r 2 * , in particular, if the parameters of the equations of motion are not within the range of the first order perturbative treatment.
We use eq.( 29) for tracking of stable and unstable states without resorting to a quasistationary change of the control offsets h (c) R and h (s) R .Figure 8 shows the result for the Duffing oscillator in a strongly nonlinear regime, ε = 1, with the parameter setup as in figure 3.Even though we are well beyond the validity of the perturbative regime our control design performs reasonably well, allowing to track the unstable state within the region of bistability.9) and eq.( 10), for ε = 1, with parameter settings as in eq. ( 11), and with the control scheme defined in eq. ( 25) and (21).Control offsets for the tracking have been taken from eq. ( 29) (with r 2 * serving as the curve parameter).Left: Successful tracking of the stationary state r 2 ∞ (blue, full symbols) in dependence on the amplitude of the driving field, see eq. (12).For comparison the corresponding data of the oscillator without control are shown as well (open symbols, cyan), see figure 3. Right: resulting amplitudes of the driving field, eq. ( 25), when control has been successful.Cyan h There are visible deviations in figure 8 between the fixed point without control and the fixed point subjected to control, in particular, close to the fold instabilities which bound the region of bistability.In addition, the unstable branch shows a seemingly subtle structure.These deviations are caused by the control action resulting in a non-vanishing value for h (s) , see figure 8, that means in a non-vanishing phase for the driving field, see eq. ( 12).Hence our observation for the amplitude r 2 * under control does not correspond to the setup used in the system without control, see figure 3, where we have h (s) = 0.If the orbit were a plain harmonic such an effect would not matter.Since the orbit in phase space is not a perfect ellipse the measured amplitude depends on the actual phase of the driving field.
There are a couple of ways to compensate for this effect.On the one hand we can slowly tune h ∞ becomes zero, which essentially amounts to an experimental root finding problem.One has to keep in mind, however, that beyond the perturbative regime the rotational symmetry shared by the low order perturbation expansion eq. ( 20) is not valid any longer.
Therefore the stability of the unstable state may get lost during this adjustment.It may in fact happen that within our design successful control for a non-vanishing phase of the driving field does not per se translate into successful control for vanishing phase of the driving field when the offsets are readjusted accordingly.On the other hand, there is no need to perform such an adjustment.The phase of the driving field, eq. ( 12), can be compensated for if we consider the stroboscopic map at a different time, i.e., if we formally chose a different cross section.If we record data at times n2π/ω + θ/ω, that means if we compute a renormalised amplitude r 2 n based on x(n2π/ω + θ/ω) and v(n2π/ω + θ/ω) then such a value of r 2 n is identical to the value one obtains from an ordinary stroboscopic map in a Duffing oscillator where the driving field has vanishing phase θ = 0.That means, while we still base control on discrete time points n2π/ω we use a suitable time lag in the data recording to compensate for the non-vanishing final value of h (s) ∞ .By this trivial change we obtain a perfect match of the controlled and the uncontrolled data, see figure 9, where even subtle details of the unstable branch are detected with considerable accuracy.It is again worth to stress that the control, the data processing, and the control based continuation can be easily implemented in experiments as only a plain time series is required.12), the estimate of the stationary amplitude r 2 * has been computed from time series data x(n2π/ω+θ/ω) and v(n2π/ω+θ/ω) with a time lag given by the phase of the driving field.Stationary amplitude as a function of the amplitude of the driving field with control (full symbols, blue).For comparison the corresponding data without control, see figure 3, are shown as well (open symbols cyan).

Data-driven control design
Our control design, eq. ( 25), was inspired by the first order perturbation expansion of the stroboscopic map, eq. ( 20).This design has been quite successful even beyond the perturbative regime, where eq. ( 20) fails to model the dynamics properly, as exemplified by the results of  24) and (30) to the data shown in figure 3, i.e., least square fit for the nonlinear resonance line of the driven Duffing oscillator, eqs.( 9) and (10) with parameters as in eqs.( 11), h (s) = 0, and ε = 1.Line: eq. ( 24) with eq. ( 30) and w 0 = −0.2338,w 1 = 0.3290.Symbols: data obtained from a parameter upsweep and parameter downsweep (cf. figure 3).
driven Duffing oscillator in the parameter setting specified above.Even though there is no a priori guarantee that the approach will work, as the perturbation expansion eq. ( 20) definitely does not capture the dynamics any longer, the result shown in figure 11 demonstrates that control works quite robust even for an initial condition which is not close to the target state.Even though the parameter setups used to produce the data in figures 6 and 11 coincide, the stabilised fixed points differ slightly as the effective force w(r 2 ), i.e., the actual control feedback is not the same, so that the final driving amplitudes h (s/c) ∞ in both cases differ slightly.
In addition to time traces we have as well computed the basin of the control, obtained with the data-based control design.The results shown figure 12 indicate in fact a considerable improvement as compared to figure 7 where the design was based on the first order perturbation scheme.Hence, the approach outlined in this section, which is based on the data obtained from a simple parameter upsweep and downsweep provides a promising strategy for a wider class of driven oscillator systems.One may even exploit the dependence of the parameter sweeps on the driving frequency ω in conjunction with the analytic expression eq. ( 24) to explore the damping mechanism in more detail.Details will be published elsewhere.

Conclusion
We have succeeded with our main aim to design a globally robust non-invasive control scheme for the stabilisation of periodic orbits, to enable control-based continuation.The control is based on the measurement of a single data point per period of the drive, so that such a scheme is applicable in fast systems like atomic force microscopy, when no extensive data processing can be done during control, and where no accurate mathematical model is available for preprocessing.
The design of our control scheme was initially inspired by the weakly nonlinear analysis of oscillator systems.Having access to the phase of the driving field has turned out to be a key for the success of the control, since thereby we have been able to ensure for large basins of attraction for the stabilised state.The control feedback contains an effective averaged force, so that the design already gives vital information about the underlying dynamics.Numerical simulations indicate that the control scheme works quite well beyond the perturbative regime, and further systematic numerical studies look promising.Above all, we have shown that the required details of the control design can be obtained from a simple scan of the nonlinear resonance curve, so that in fact no underlying mathematical model is needed from the outset.
We have based our analysis on a one degree of freedom mechanical oscillator model with fairly general potential.For the theoretical analysis we have used the assumption that higher order harmonics play a limited role, even though our numerical studies show that such a constraint can be relaxed to some extent.However, it still needs to be investigated how the present approach can be generalised to higher dimensional driven dynamical systems.Again weakly nonlinear analysis could provide a hint how to proceed in this context.
We have here considered a model with a simple viscous damping.In real world experiments, such as atomic force microscopy, the treatment of all losses by a single viscous damping constant γ may be regarded as undercomplex, even in the seemingly simple situation of a stiff, hydrophobic sample in a vacuum environment.While the dynamics of an oscillating microcantilever beam itself behaves largely linear, the complexity of the system arises from the largely unknown dynamics of the multitude of effects taking part in the nanoscopic junction between probe and sample.These unknowns and nonlinearities do pose the challenges.Currently, it is hard to say what physics or (bio)chemistry happens at the turning point of the tip in the vicinity of the sample, the point of strongest interaction which shapes contrast and stability of the measurement.While force-distance spectroscopy can help to unravel force contributions from molecule layering, visco-elasticity, plastic deformation, electrical polarisation and attraction, electro-chemical reactions, rupturing molecular bonds, entropic interaction, depletion forces, oxidation, receptor-ligand, or other (see, for instance [31,32]), this is of limited value to judge the situation of topography acquisition in the dynamic mode at cantilever frequencies in the 10 to 500 kHz regime.At a phenomenological level the complex dissipation processes can be modelled by a state dependent damping, and its impact can be analysed within a weakly nonlinear perturbation expansion as well, resulting in an additional effective damping term which along the side of the effective averaged force enters the shape of the nonlinear resonance curve.Both contributions, the effective damping and the effective force can be disentangled by monitoring in addition the frequency dependence of the nonlinear resonance.Hence, by having access to stable as well as unstable branches we can quite accurately determine potential and damping from measured data and thus contribute, for instance, to the outstanding challenge to understand the dissipative mechanisms in atomic force microscopy.
A robust control scheme which enables large basins of attraction is one of the keys to implement model free continuation of bifurcations in experiments, to make the power of continuation tools in the study of mathematical models available as a 21st century data-based spectroscopic tool.Our control design meets all these constraints, so that we can track unstable states even when no quasi-stationary parameter sweep can be implemented.Having developed a suitably robust control scheme we have taken the next step towards finally implementing control based continuation in real world complex experiments.Details in that direction will be reported elsewhere.

Figure 1 :
Figure 1: Dynamic force microscopy of a conductive polymer blend (PEDOT:PSS) thin film deposited by spin coating on a glass substrate (rotational frequency is 5000 revolutions per minute, film thickness ≈ 32 nm).Cantilever type: Nanosensors SSS-NCHR with a force constant in the regime of 10-13 N/m and a nominal tip radius of 2 nm.AFM instrument: Park Systems NX-20 (in ambient air).AM amplitude drop A/A f ree ≈ 0.75.Left: topography, middle: phase between excitation and cantilever bending, right: amplitude.The phase exhibits a clear bimodal distribution peaking at -33 • and +14 • (not shown).

Figure 3 :
Figure 3: Bifurcation diagram of eqs.(9) and(10) with ε = 1, h (s) = 0, and other parameters as in eq.(11).Data have been obtained from a numerical simulation of the stroboscopic map, i.e., from time traces evaluated at integer multiples of the period of the driving field.Skipping a transient of 50 iterations the orbit settles on a stable fixed point (x * , v * ).The dependence of r 2 * = x 2 * + v 2 * /ω 2 on h (c) is shown for a quasi-stationary parameter upsweep (cyan, full symbols) and a parameter downsweep (red, open symbols).The black line shows the analytic result obtained from a first order perturbation theory, see eq. (24).
finally a control based data driven continuation of bifurcations in experiments such as atomic force microscopy.

Figure 4 : 2 * = x 2 * + v 2 *
Figure 4: Bifurcation diagram of eq.(9) for ε = 0.05 (other parameters are as in figure3).Data are obtained from a numerical simulation of the stroboscopic map (red, open symbols), or using iterates of the analytic map derived by first order perturbation theory eqs.(20) and (21) (cyan, full symbols).Skipping a transient of 800 iterations the orbit settles on a stable fixed point (x * , v * ).The dependence of r 2 * = x 2 * + v 2 * /ω 2 on h (c) is shown for a quasi-stationary parameter upsweep (left) and a parameter downsweep (right).The black line shows the analytic expression obtained for the fixed point at first order perturbation theory, see eq. (24).
design select the components of the fixed point to be stabilised.The amplitudes h (c/s) n settle on values of the driving amplitude h (c) ∞ = 0.128 . . .and h (s) ∞ = −0.0002 . .., which correspond to the unstable branch right in the middle of the bistable region, cf.figure 4. We have also checked how the control performs for different initial conditions.For initial conditions in the range −5 ≤ x(0) ≤ 5, −5 ≤ v(0) ≤ 5 we always find successful stabilisation of the unstable periodic state, in line with what we expect from the first order perturbation treatment.For larger values of the initial condition occasionally solutions seem to diverge, which is far from surprising given that the equations of motion and the control scheme contain cubic terms.

Figure 8 :
Figure 8: Control based continuation of the stationary state in the Duffing oscillator, eq.(9) and

Figure 9 :
Figure 9: Control based tracking of the stationary state in the Duffing oscillator with parameters and control design as in figure 8. To compensate for the non-vanishing phase of the driving field, see eq. (12), the estimate of the stationary amplitude r 2 * has been computed from time series data x(n2π/ω+θ/ω) and v(n2π/ω+θ/ω) with a time lag given by the phase of the driving field.Stationary

Figure 10 :
Figure 10: Least square fit of the analytic expression eqs.(24) and(30) to the data shown in figure3,