Controlling coexisting attractors of a class of non-autonomous dynamical systems

This paper studies a control method for switching stable coexisting attractors of a class of non-autonomous dynamical systems. The central idea is to introduce a continuous path for the system's trajectory to transition from its original undesired stable attractor to a desired one by varying one of the system parameters according to the information of the desired attractor. The behaviour of the control is studied numerically for both non-smooth and smooth dynamical systems, using a sof-impact and a Duffing oscillator as examples. Special attention is given to identify the regions where the proposed control strategy is applicable by using the path-following methods implemented via the continuation platform COCO. It is shown that the proposed control concept can be implemented through either using an external control input or varying a system parameter. Finally, extensive numerical results are presented to validate the proposed control methods.


Introduction
Non-autonomous dynamical systems exhibit a rich variety of different long-term behaviours coexisting for a given set of parameters, which is referred as multistability or coexisting attractors. It is a common phenomenon in science and nature that relies crucially on the initial conditions of the system. For example, multistability can be found in many engineering systems, such as the vibro-impact capsule [1,2], the electronic circuits [3,4], the gas laser [5] and the drilling system [6,7]. It is also a fundamental property of biological systems, such as the spiking neurons [8], tumor progression [9], cell fate transitions [10,11] and cell cycle control [12,13]. Control of multistability has been an active research field in the past few decades, see [14]. These coexisting attractors are extremely sensitive to noise due to their fractally interwoven basins of attraction. Switching the coexisting attractors while protecting against perturbation-induced basin hopping is vital in practice, particularly for engineering applications. Since the engineering systems may present different performance at different coexisting attractors (e.g. energy efficiency), control the switching among these stable attractors could offer them more flexibility. For example, stick-slip motion and constant rotation coexist in the drilling system [15,16], and prevention of the stick-slip motion can maintain the drilling at a high penetration rate so reducing the operation cost. For this reason, numerous methods have been developed to control multistability, and a control strategy for avoiding perturbation-induced attractor switching is always desirable.
In the present work, we will develop a novel control method for controlling the coexisting attractors of a class of non-autonomous dynamical systems, including both smooth and non-smooth dynamical systems. It is well known that the OGY method [17] was initially developed for stabilizing an unstable periodic orbit embedded in a chaotic motion via adjusting the system parameter in a small neighbourhood. For the same purpose, Pyragas [18] developed a delayed feedback controller to achieve the stabilisation of unstable periodic orbits in a chaotic system. In [19], Pyragas and Tamaševičius used the delayed feedback control method to stabilise an analogue circuit. In [20,21], Pyragas and Pyragas introduced the act-and-wait concept to reduce the dimension of phase space of the systems with the delayed feedback control. Although the above studies focused on stabilising unstable periodic attractors embedded in the chaotic attractor, those control concepts also can be applied to control coexisting attractors. For example, Lai [22] constructed a hierarchy of paths for targeting the desired attractor by introducing a feedback perturbation to the system with fractal basins. Wang et al. [23] adopted the similar concept to control an undesired attractor to a desired one for a low dimensional network system. In [24], Zhang et al. used the delayed feedback control to switch the unwanted attractors to a period-1 attractor for a soft impacting oscillator, and in [25], Páez Chávez et al. developed a numerical approach for analysing the dynamical properties of the non-smooth systems with time delays. On the other hand, there are also many other methods for controlling the multistability of dynamical systems. Arecchi et al. [26] discovered that the external noise can bridge the coexisting states of a forced Duffing equation. In [27], Pisarchik and Goswami employed a slow external periodic perturbation to the system parameter of a bistable system to annihilate its coexisting attractors. In [28], an intermittent control was designed to provide an impulsive force for non-autonomous dynamical systems to switch between coexisting attractors, which was verified by both numerical and experimental results. Liu and Páez Chávez [29] developed a linear augmentation control law to control the multistability of a soft impacting oscillator, and analysed the dynamical properties of the control law by using the path-following (continuation) techniques for non-smooth dynamical systems. Nevertheless, it should be noted that most of these control methods for attractor's switching of non-autonomous dynamical systems relied on external input, while few works have concerned about the switching by utilising system's properties, e.g. modulation of system parameter. It is obvious that the latter approach is less invasive, so the original system could be maintained, which is advantageous to some systems whose external input is hard to access, such as the machining process [30][31][32]. Therefore, in the present work, we will focus on the control method that relies on the original properties of the system only to achieve the switching among coexisting attractors.
According to [14], Pisarchik and Feudel classified the existing methods for controlling multistability into three categories: feedback, non-feedback and stochastic controls. They have suggested that the most efficient way of ensuring a predefined behaviour for the system is to annihilate all the other coexisting attractors. However, annihilation of undesired attractors in some dynamical system could change the existing structure of solutions leading to the emergence of new complex basins of attraction. One of the simplest ways to achieve this could be to apply an impulsive external perturbation to direct system's trajectory from one basin to another one [14]. In [33], Kaneko's study suggested that a short pulse can cause system state to jump from one attractor to another one if the pulse's amplitude is sufficiently large. Chizhevsky et al. proposed to apply a short-pulsed perturbation in the form that switching the system off and on for a very short time, so the system can run from a different initial condition. The problem of these methods is that the short pulse was applied in the form of non-feedback. If the amplitude of the pulse is small, the system regains the same attractor after a few iterations. In addition, these methods are only effective for the systems with few coexisting attractors, and the short-pulsed control may become uncertain if more attractors coexist with fractal basins. Hence, continuous control by utilising system's feedback signals is more preferable for non-smooth dynamical systems that have considerable numbers of coexisting attractors, e.g. in the near-grazing dynamics [34].
In the present work, the primary focus is to address the continuous switching between two of coexisting stable attractors by varying a system parameter without affecting their original dynamics. In order to achieve this, we propose a continuous control method that can adjust a system parameter based on the information of trajectory of the desired attractor. This control method is applied to the controlled system continuously until its trajectory is sufficiently close to the desired one. We develop two control strategies based on this concept, the so-called linear and nonlinear control strategies, where the former is implemented through an external control input and the latter is applied via a system parameter. The advantage of the nonlinear control strategy is that it depends only on the original properties of system parameter and does not rely on any external input. To demonstrate the applicability of the proposed control method to both non-smooth and smooth systems, an impact and a Duffing oscillators were employed in the present paper.
The rest of this paper is organized as follows. Section 2 introduces the mathematical model of the periodically forced mechanical oscillator subjected to a one-sided soft constraint and some mathematical preliminaries for the proposed control method. In Section 3, the effectiveness of this control method on a single-degree-of-freedom oscillator with piecewise-smooth nonlinearity is studied. In Section 4, the proposed control method is adopted to the Duffing oscillator. Finally, concluding remarks are drawn in Section 5.

Design of feedback control strategy
We consider periodically forced systems with a single control input of the forṁ where the input u(τ ) is a scalar function of time, andẎ denotes differentiation with respect to time τ . We assume that the uncontrolled system, (1) with u(τ ) = 0, has an attractor Y d (·) (where the subscript in Y d stands for "desired"), but that the initial condition Y u,0 is away from Y d (0), and possibly outside the basin of attraction of Y d (·). We also assume that this "desired" attractor Y d (·) is (internally) stable in the sense that it has no positive Lyapunov exponents. In practice our method is intended to be applied to periodic orbits. Our test examples will be single-degree-of-freedom oscillators with periodic forcing and multiple coexisting attractors. The typical scenario we envisage is that Y u,0 is on one of the other ("undesirable") attractors of the uncontrolled system.

Single-degree-of-freedom oscillator with piecewise-smooth nonlinearity
In this section, we will use the single-degree-of-freedom oscillator with piecewise-smooth nonlinearity shown in Fig. 1 as an example to study the proposed control method. Soft impacts occur in mechanical systems when an object hits an obstacle with a negligible mass but with a non-negligible stiffness, see e.g. [35][36][37][38][39]. As can be seen from Fig. 1, it is assumed that the discontinuity boundary is fixed at x = e, with e > 0 being the nondimensional gap. The equations of motion of the oscillator are in form (1), where Y (τ ) := (x(τ ), v(τ )) T . We will consider three different cases. Defining for scalarũ the three cases for control input are The first case, F lin has linear control input and the second case, F a , has parametric control input, varying the forcing amplitude, while the third case adjusts the gap e. The function H(·) stands for the Heaviside step function. In the right-hand side F frc , defined in (2), the variables (including time) and parameters of the system are nondimensionalised according to where y 0 > 0 is an arbitrary reference length, ω n is the natural frequency, ω is the frequency ratio, β is the stiffness ratio, ζ is the damping ratio, e is the nondimensional gap between the mass and the secondary spring and a is the nondimensionalised amplitude of the external excitation. The right-hand side given in Eq.
(2) is a typical non-autonomous dynamical system with the one-sided elastic constraint considered as the non-smoothness that can lead to complex phenomena, such as the grazing bifurcation [34,40], the coexistence of multistable attractors [41,42] and the chaotic motions [37]. Here, we consider the situation that system (1) with u(τ ) = 0 has many coexisting attractors within some specific ranges of the system parameters. In particular, we consider a stable attractor Y d with nonpositive Lyapunov exponents exists, which means this attractor is periodic or quasi-periodic. Usually, this Y d (·), our desired attractor, is a stable periodic orbit. For Y d , small changes in system parameters do not affect the response of the system significantly, such that the attractor Y d will persist.

Distance-reducing feedback control
Let us introduce the following definition.
Define the difference between the desired and the current states as is the distance at time τ from the desired attractor Y d . We will investigate a simple feedback control strategy that aims to reduce this distance ∆ over time, such that [d/dτ ]∆(τ ) < 0. Requiring that u(τ ) has permissible values between −M 1 and M 1 and Lipschitz constant M 2 , we adjust u at every time τ ≥ τ 0 according tou starting from u(τ 0 ) = 0, for as long as ∆(τ ) > ǫ with some tolerance ǫ ≪ 1. In Eq. The term [d/dτ ]∆(τ ) is a function of (τ, Y u (τ ), Y d (τ ), u(τ )), making the derivative with respect to u non-zero: In the practical algorithms described in Sections 2.3 and 2.4, we will apply Eq. (6) only when [d/dτ ]∆(τ ) > 0, otherwise, we will setu(τ ) = 0.

Implementation with finite sampling step -linear case
Let us assume the sampling time step, h > 0, such that τ i = τ 0 + ih. If the control input u enters the right-hand side linearly with a constant coefficient vector b, such that F (τ, Y, u) = F (τ, Y ) + bu (as in example F lin given in (2), where b = (0, 1) T ), the notation in the definition (6) of the control u and the resulting expression (7) simplifies. Thus, we formulate Algorithm 1 for this common case separately. In our formulation of Algorithm 1 with the sampling step h, the control input u(τ ) is of type "zero-order hold". That is, u(τ ) is a constant u i on the sampling interval [τ i , τ i+1 ], and the Lipschitz constant M 2 applies to changes per time step: |u i+1 − u i | ≤ hM 2 .
After introducing the algorithm, the following main theorem can be obtained.
Theorem 2.1. Let the dynamical systemẎ (τ ) = F (τ, Y ) + (0, u(τ )) T with u = 0 have two stable coexisting attractors Y c and Y d , and let M 1,2 be bounded intervals in R. We assume that there exists a time τ * > τ 0 , such that for all sufficiently small sampling steps h the control U and trajectory Y u

Algorithm 1 Linear control
Step 0: Choose M 1 and M 2 , where M 1 is the boundary of |u(τ )|, M 2 is the boundary of |u(τ )|, andu(τ ) := du(τ ) dτ . Take the initial control u(τ 0 ) = 0 ∈ M 1 , and set the iteration index i = 0 and the time step h. Below Step 1: Compute the range ofu(τ i ) that satisfies the following criterion calling this range the feasible range foru(τ i ).
Step 2: If this feasible range is greater than M 2 and ensure that u i+1 satisfies M 1 , take the minimum value of |u(τ i )| and go to Step 3. Otherwise, choose the valueu(τ i ) such that both M 1 and M 2 are satisfied and is the closest to the feasible range ofu(τ i ), and go to Step 3.
Step 3: . . , n * ). Then these two stable attractors are controllable by this external controller.
According to Theorem 2.1, the following lemma can be obtained.

Implementation with finite sampling step -nonlinear case
By adjusting the accessible parameter of non-autonomous dynamical systems, such as the amplitude of excitation, system energy is altered, so the dynamical property of the system can be controlled (see e.g. [17,22,23]). In this subsection, we will apply the control method in Algorithm 1 to an accessible system parameter to achieve the switching between two stable coexisting attractors. The detailed algorithm of the new control method (Algorithm 2) is given as below, and the following theorem is introduced. Theorem 2.3. Let the dynamical systemẎ (τ ) = F (τ, Y, u p (τ )) with u p = 0 have two stable coexisting attractors Y c and Y d , and let M p,κ , κ = 1, 2 be bounded intervals in R. We assume that there exists a time τ * > τ 0 , such that for all sufficiently small sampling steps h the control u p and trajectory Y u

Algorithm 2 Nonlinear control
Step 0: Choose M p,1 and M p,2 , which are the boundaries of u p (τ ) andu p (τ ), respectively. Take the initial control u p (τ 0 ) = 0 ∈ M p,1 , and set the iteration i := 0 and the time step h. while the termination criterion Step 1: Compute the feasible range ofu p (τ i ) to satisfy the following criterion Step 2: If this range is greater than M 2 and ensure that u p,i satisfied M p,1 , take the minimum value of |u p (τ i )| and go to Step 3. Otherwise, choose the valueu p (τ i ) such that both M p,1 and M p,2 are satisfied and is closest to the feasible range ofu p (τ i ). Then go to Step 3.
Step 3: (1) for τ > τ i . Increment i by 1 and return to Step 1. end while obtained by Algorithm 2 with bounds M p,κ satisfy . . , n * ). Then these two stable attractors are controllable by this external controller.
According to Theorem 2.3, the following lemma can be obtained.

Control of non-smooth dynamical systems
This section will show the effectiveness of the proposed control methods by using the impact oscillator shown in Fig. 1, which is a typical non-smooth dynamical system exhibiting many coexisting attractors at its near-grazing dynamics [42]. The following parameters were used for which two stable attractors, a period-2 and a period-5 responses, coexist as shown in Fig. 2. ζ = 0.01, e = 1.26, a = 0.7, β = 28 and ω = 0.85.

Linear control
Firstly, the dynamical response of the impacting system of type (2) with right-hand side F lin given in (3) is presented in Fig. 3 at where the original period-5 attractor was switched to the period-2 attractor by using the external control strategy. The time step of the simulations was fixed at h = 0.002, and the control strategy was implemented at τ = 591.358. As can be seen from Figs. 3(b) and (c), the period-5 response experienced a transition and was settled down to the period-2 response around τ = 637. Fig. 3(d) shows the transition by grey line on the phase plane and indicates the steady-state response by red line. Figs. 3(e) and (f) also demonstrate the effectiveness of the control where the distance between the present and the desired attractors was reduced once the control sequence was applied. It can be seen that the overall trend of the distance was decreased, and according to the simulation, it was about 0.026 at τ = 637 and was about nil after τ = 650. In addition, it can be seen from Fig. 3(f) that, during the control process, no control was applied when the distance between the two trajectories was decreasing. Then a continuous increase in the control signal to 0.832 at τ = 593.01 and a continuous decrease to 0 at τ = 594.65 were recorded. Thereafter, the control experienced intermittent control actions, and the amplitudes of the control actions decreased as the two trajectories were closer. Finally, the control was turned off when the control target was achieved. Fig. 4 shows the control result from the period-2 to the period-5 attractor by using the external control strategy. According to the simulation, the control was switched on at τ = 591.358, and the  distance between the two trajectories was decreased to 0.01 at τ = 643.5. The control signal became zero gradually when the distance was sufficiently small at τ = 651.9.

Nonlinear control
Next, we will test the effectiveness of the nonlinear control strategy (Algorithm 2) by varying the amplitude of excitation a and the gap e, as given in (4) and (5). Again, the control target here is to switch the response of system (1) with right-hand side (4) between the period-5 and period-2 attractors by varying its amplitude of excitation.   5 shows the control of the impacting system from the period-5 attractor to the period-2 attractor by varying its amplitude of excitation a. Based on our calculation, when the nonlinear control strategy was applied at τ = 591.358, the controlled trajectory of the system experienced a transition as observed in Figs. 5(b), (c) and (d). As can be seen from Fig. 5(e), the distance between the desired and the controlled trajectories in 2-norm was decreased indicating the controlled trajectory approached to the desired one, and this distance was reduced to 0.001 at τ ≈ 640.64. Fig. 5(f) presents the time history of the control signal u p (τ ) that initiated from τ = 591.358 and terminated at τ = 640.854. Thereafter, u p (τ ) = 0, and the amplitude of excitation was back to its original value, a = 0.7.
The control from the period-2 attractor to the period-5 attractor by varying the amplitude of excitation a is presented in Fig. 6. The control strategy was applied at τ = 591.358, but the distance between the two trajectories was not decreased continuously. Therefore, compared to the external control strategy, it took a longer time for the system to settle down to the period-5 attractor. According to the simulation, the control was switched off at τ = 689.1 when the distance was reduced to 0.001. When applying the control input u p (τ ) to the gap e, as given in (5), the control aims again to switch the response of the soft-impact system (1) with right-hand side F e (τ, Y, u p ) given in (5), between the period-5 and period-2 attractors by varying its gap e. Fig. 7 presents the control result from the period-5 to the period-2 attractor by varying system's gap e. The control strategy was applied at τ = 591.358 and was switched off at τ = 659.6. During the control period, the trajectory of the system experienced a transition, and the distance between the controlled and the desired trajectories was decreased from 1.03 to 0.001. Within the same time duration, the control signal u p (τ ) reached a maximum value u p (τ ) = 0.21 and decreased to −0.02 at τ = 630.326. Thereafter, the control signal did not change significantly and reduced to nil gradually after τ = 659.49.
To demonstrate the switching from the period-2 to the period-5 attractor by varying system's gap, Fig. 8 presents the control result. Based on the calculation, the control strategy was applied at τ = 591.358, and the distance between the two trajectories was reduced from 1.03 to 0.001 at τ = 703.794. Thereafter, the control signal u p (τ ) decreased to nil, and the control target was achieved. Compared to the control result shown in Fig. 7, the transition from the period-2 to the period-5 attractor took a longer time and had a more complex transient response, which was due to the complexity of the period-5 response.

Bifurcation analysis of the coexisting attractors
In this section we will study in detail the effect of the control parameters (excitation amplitude a and mass-spring gap e) on the period-2 and period-5 coexisting attractors studied in the previous section. For this purpose we will employ path-following methods for piecewise-smooth dynamical systems, using  the continuation platform COCO [43]. The precise COCO-implementation for the impact oscillator (2) can be found in a previous publication by the authors [29], which will be adopted in the present work. Contact time As can be seen from the previous section, in order to apply the proposed control mechanism it is essential to identify parameter regimes where the considered period-2 and period-5 attractors maintain both their stability properties and orbit structure. For this purpose we will carry out a one-parameter continuation of the underlying attractors with respect to the excitation amplitude a, using the contact time as solution measure, i.e. the time the impacting mass stays in contact with the secondary spring per orbital period.
The result of the process described above is depicted in Fig. 9, which shows yellow and green curves corresponding to the numerical continuation of the period-2 and period-5 attractors detected in Fig. 2, respectively. In both cases, the parameter window where the corresponding periodic solutions remain stable and maintain their orbit structure is determined by period-doubling and grazing bifurcations of limit cycles. Specifically, the period-2 solution traced along the yellow branch loses stability when the excitation amplitude decreases below a ≈ 0.63111, where the solution undergoes a period-doubling bifurcation (PD1). Here, the original period-2 orbit (with one impact per orbital period) becomes unstable and a family of period-4 orbits is born, see for instance the test solution plotted in Fig. 9(f), right after the bifurcation occurs. On the other hand, when the parameter a increases, a grazing bifurcation is found at a ≈ 1.54462, where the solution makes tangential contact with the impact boundary x = e, see Fig.  9(g). After this point, a small window of period-2 solutions with two impacts per orbital period exists, and they lose stability via a fold bifurcation at a ≈ 1.54486 (not shown in the diagram). An analogous ag replacements scenario is found for the period-5 attractor (with three impacts per orbital period) depicted in Fig. 9(c). As before, the window of stability (and orbit structure preservation) for this solution is determined by the period-doubling bifurcation PD2 (a ≈ 0.64564) and the grazing point GR2 (a ≈ 0.71258), which results in a significantly smaller window than the one obtained for the period-2 attractor.
With the results of the one-parameter continuation we are now in position to determine a parameter region in the a-e plane where the considered period-2 and period-5 attractors maintain both their stability properties and orbit structure. To this end, we will perform a two-parameter continuation of the codimension-one bifurcations detected above. Fig. 10(a) shows the locus of the period-doubling points PD1 (blue curve), PD2 (black curve) and grazing bifurcations GR1 (green curve), GR2 (red curve) encountered in Fig. 9(a). In this figure, two regions are highlighted, in grey and yellow colors. The grey area represents the parameter region in which the stable period-2 solution (with one impact per orbital period, see panel (b)) exists. The yellow region corresponds to the coexistence of the latter solution with the stable period-5 orbit (with three impacts per orbital period, see panel (c)). Furthermore, several test points have been selected in order to illustrate the validity of the highlighted yellow area in the a-e plane. Specifically, pairs of coexisting attractors have been computed at the test points P1 (a = 0.68, e = 1.26), P2 (a = 1.1, e = 2.05) and P3 (a = 1.5, e = 2.8), see panels (b)-(g) in Fig. 10. In this way, the yellow area can be used as a reference for the applicability of the proposed control scheme, so as to guarantee that the parametric perturbations do not bring the system to a regime where either of the considered attractors lose stability or the intended orbit structure.

Controlling multiple coexisting attractors
Here, we consider achieving the switch among three coexisting attractors, as shown in Fig. 11, by the nonlinear control strategy (Algorithm 2) through varying the amplitude. At the beginning, the control of the impacting system from the initial attractor (period-7 attractor with large amplitude) to desired attractor (period-7 attractor with small amplitude) is considered, as shown in Fig. 12(a)-(c). In details, in Fig. 12(a), the distance between the desired and the controlled trajectories in 2-norm was decreased and finally was reduced to nil. Fig. 12(c) presents the time history of the control signal u p (τ ) that initiated from τ = 589.417 and terminated at τ = 763.2234. The transition on the displacement can be observed from Fig. 12(b). Secondly, the control from the period-7 attractor with small amplitude to the period-3 attractor by varying the amplitude of excitation a is presented in Fig. 12(d)-(f). Specifically, the control strategy was applied at τ = 589.417 and was switched off at τ = 657.5665. During this period, the trajectory experienced a transition, as shown in Fig. 12(e), and the distance between the controlled and the desired trajectories was decreased from 0.182 to 0.001, as shown in Fig. 12(d). In the meanwhile, the control signal did not change significantly and reduced to nil gradually. After that, the controlled trajectory approached the desired trajectory spontaneously. Finally, Fig. 12(g)-(i) present the control result from the period-3 attractor to period-7 attractor with large amplitude by varying the amplitude. The control strategy was also applied at τ = 589.417, and switched off at τ = 666.3881 as shown in Fig. 12(i). Within the same time duration, the trajectory of the controlled system approached to the desired trajectory in Fig. 12(h), and the distance between the controlled and the desired trajectories was decreased from 0.403 to 0.001 as shown in Fig. 12(g). The above processes show that the proposed method can present good performances on controlling more coexisting attractors.

Nonlinear control
In this section, the Duffing oscillator representing smooth dynamical systems is employed to test the versatility of the proposed control method. The Duffing system, which is known to have many coexisting attractors without control, can be described by where Y (τ ) := (x(τ ), v(τ )) T , and The following parameters: Γ = 1.9, ω = 1.2, p 1 = 0.9 and p 2 = 1 were considered in this study. At these parameter values, the system without control (u p = 0) has two coexisting attractors, depicted in Fig. 13, which are a period-1 small and a large amplitude attractors with their Poincaré sections denoted by black and violet dots, respectively. The control aims for system (12) to switch the two stable attractors shown in Fig. 13 by varying the stiffness of the nonlinear spring p 2 . The control result for the switching from the large to the small amplitude attractor is shown in Fig. 14, where the control strategy was applied to the original attractor at τ = 418.879, and the controlled trajectory experienced a transition until τ = 439. During this time the 2-norm distance between the control and the desired trajectories was reduced from 1.338 to 0.001. The control signal reached the maximum u p (τ ) = 0.3 and the minimum u p (τ ) = −0.3 for several times before it was switched off.
The control from the small to the large amplitude attractor is presented in Fig. 15, where the control strategy was switched on at τ = 418.879 and was switched off at τ = 440.246 when the distance between the two trajectories was decreased to 0.001. Compared to the switching in Fig. 14, the transition from the small to the large amplitude attractor took a longer time.

Bifurcation analysis of the coexisting attractors
Analogous to Section 3.3, in this section our main concern will be to study in detail the effect of the control parameters p 1 , p 2 on the small-and high-amplitude oscillations of the Duffing system (12), see Fig. 13. To this end, we will employ path-following methods for limit cycles, implemented via the continuation platform COCO [43], along with its routines for bifurcation detection and two-parameter continuation of codimension-1 bifurcations.  The starting point for our study is the high-amplitude periodic solution shown in Fig. 16(c), computed for p 1 = 0.8. Panel (a) presents the result of the numerical continuation of this orbit with respect to the control parameter p 1 . In this diagram, changes of stability are detected, which are marked with solid (for stable solutions) and dashed (unstable solutions) lines. The window of stability of the the high-amplitude orbit is bounded from above by the fold bifurcation F2 (p 1 ≈ 0.90352). At this point, a branch of unstable periodic solutions is born, which finishes at the fold point F1 (p 1 ≈ 0.69494). Here, a family of stable oscillations emerges, corresponding to small-amplitude periodic orbits as can be seen at the test point P2, see Fig. 16(b). Consequently, the bifurcation points F1 and F2 defines a parameter window where both attractors coexist.
Next, we will carry out a two-parameter continuation of the fold points detected above in order to determine a region in the p 1 -p 2 plane where the small-and high-amplitude attractors of the Duffing system coexist. The result of this numerical process is presented in Fig. 16(d), where the red curve stands for a locus of fold bifurcations of limit cycles. In this picture, the yellow area enclosed by the fold curve represents the parameter regime where small-and high-amplitude oscillations coexist. The intersections of the horizontal dashed line (p 2 = 1) with the bifurcation diagram correspond to the fold bifurcations F1 and F2 found in Fig. 16(a). Furthermore, the numerical computations reveal the presence of a codimension-2 point (p 1 , p 2 ) ≈ (1.14902, 1.51194) (CP), where two branches of fold points (those corresponding to F1 and F2) join together via a cusp singularity. In this way, it is possible to determine boundaries in the considered parameter region for the application of the control mechanism proposed in this work.

Concluding remarks
This paper studied a new control method for switching stable coexisting attractors of non-autonomous smooth and non-smooth dynamical systems. Our control aim was to control an undesired coexisting attractor to a desired one by modulating a system parameter without affecting the original property of the system. To examine the proposed control concept, we implemented two control strategies with finite sampling step, namely the linear and nonlinear control strategies, where one was implemented through the external control input and the other one was applied via a system parameter. In the first part of our simulation work, two multistable scenarios (one is coexisting two attractors and another is coexisting three attractors) of the impact oscillator were studied. Our simulations show that both control strategies are effective for switching the stable coexisting attractors in the impact oscillator. The effective control region for the control parameters (excitation amplitude a and mass-spring gap e) was also found by employing the path-following methods for piecewise-smooth dynamical systems. In the second part of the simulation work, we implemented the nonlinear control strategy to the Duffing oscillator for switching a period-1 small and a period-1 large amplitude attractors. Path-following methods for limit cycles were used to identify the effective control region of the proposed control method.
Compared with the classical delay feedback control proposed by Pyragas [18], our proposed methods have many advantages on the control of coexisting attractors. First of all, the main advantage of the proposed method for control of coexisting attractors compared to delayed feedback control is its transient behaviour. Delayed feedback control locally stabilizes a target periodic orbit that is unstable without control. Its behaviour for initial conditions far away from the target is not considered in its design and, thus, can lead to undesirable arbitrarily long transients. In contrast, the proposed method contains terms explicitly driving the trajectory toward its target. Secondly, especially for periodically forced systems the delay feedback control is infeasible for switching from an initial attractor to a target with equal (or multiple) period. This is apparent in the examples in sections 3.4 (period 7 to period 7), 4.1 and 4.2 (both period 1 to period 1), where delayed feedback would be non-invasive (and, hence, ineffective) on the initial condition. Thirdly, the proposed methods do not result in a significant change on the system's states due to restricting the varying rate of the control signal. On the contrary, when the delay feedback control is just introduced into the system, the system's states have to witness a significant change due to the value of the control is very large at the beginning. Finally, the nonlinear control strategy depends only on the original properties of system parameter, and does not need to introduce any external inputs. Besides the above differences, the common point between the proposed methods and delay feedback control is needing the information of derivatives of the right-hand sides of the ODE. In details, although the proposed methods require estimates of these derivatives to generate the control signal in its implementation and the delay feedback control does not, the feasible range of control gains for which delayed feedback control is stabilizing is still limited by the same derivatives, due to the stability of controlled trajectory depends on the Jacobian [44].
Next, we consider the Taylor expansion of F (τ, Y u (τ ), u p (τ )) within the time interval as By following the same procedure in Theorem 2.1, this theorem can be proved.