Delayed feedback control of unstable steady states with high-frequency modulation of the delay
Aleksandar Gjurchinovski, Thomas Jüngling, Viktor Urumov, Eckehard Schöll
aa r X i v : . [ n li n . C D ] A ug Delayed feedback control of unstable steady states with high-frequency modulation ofthe delay
Aleksandar Gjurchinovski, ∗ Thomas J¨ungling, Viktor Urumov, and Eckehard Sch¨oll † Institute of Physics, Faculty of Natural Sciences and Mathematics,Sts. Cyril and Methodius University, P. O. Box 162, 1000 Skopje, Macedonia Institute for Cross-Disciplinary Physics and Complex Systems, IFISC (UIB-CSIC),University of the Balearic Islands, 07122 Palma de Mallorca, Spain Partenij Zografski 46, 1000 Skopje, Macedonia Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, 10623 Berlin, Germany (Dated: June 10, 2018)We analyze the stabilization of unstable steady states by delayed feedback control with a periodictime-varying delay in the regime of a high-frequency modulation of the delay. The average effectof the delayed feedback term in the control force is equivalent to a distributed delay in the intervalof the modulation, and the obtained distribution depends on the type of the modulation. In ouranalysis we use a simple generic normal form of an unstable focus, and investigate the effects ofphase-dependent coupling and the influence of the control loop latency on the controllability. Inaddition, we have explored the influence of the modulation of the delays in multiple delay feedbackschemes consisting of two independent delay lines of Pyragas type. A main advantage of the variabledelay is the considerably larger domain of stabilization in parameter space.
PACS numbers: 05.45.Gg, 02.30.Ks
I. INTRODUCTION
The possibility to stabilize unstable periodic or sta-tionary states embedded in chaotic attractors has beenelaborated more than two decades ago [1]. The main ideaconsists of using small external perturbations to forcethe system to follow one of its stable manifolds. Theseperturbations are applied at specific instances when thechaotic trajectory is close to the desired periodic orbit.This seemingly straightforward theoretical concept hascaused a revolution in applied nonlinear science. It hasbeen realized that the control of chaos could have a sig-nificant outcome in real-world experiments, where onecould generate different kinds of ordered behavior froman utterly erratic one [2].A conceptually simple method to stabilize unstableequilibria and periodic orbits is the time-delayed feed-back control (TDFC) introduced by Pyragas [3, 4]. Here,the perturbation has the form of a continuous feedbackconstructed from the difference between some suitablescalar signal obtained from the system and the same sig-nal delayed by a constant time τ . The difference signal isamplified by the gain factor K and then re-injected intothe original system. For a certain choice of the feedbackgain K and the delay time τ , the control of the unstablestate can be realized, in which case the feedback forcevanishes by construction, making the method essentiallynoninvasive.Since a detailed knowledge of the target state is notrequired and the controller is very robust with respect to ∗ Electronic address: [email protected] † Electronic address: [email protected] noise, the Pyragas control has become one of the mostpopular control methods used in experiments and eventechnological applications. The application is quite di-verse, and includes stabilization of unstable states inelectronic chaotic oscillators [5–7], mechanical pendu-lums [8, 9], laser systems [10–12], electrochemical sys-tems [13, 14], drift waves in a magnetized laboratoryplasma [15], chaotic Taylor-Couette flow [16], cardiacsystems [17, 18], ferromagnetic resonance systems [19],gas discharge systems [20, 21], controlling helicopter ro-tor blades [22], controlling the walking mechanism of arobot [23–25], stabilization of cantilever oscillations in anatomic force microscope [26], and controlling librationalmotion of a tethered satellite system in an elliptic orbit[27], amongst others. For a comprehensive review of thetechnical implementation of the method, see [2, 4].In parallel to experimental realization, substantialwork has been done to understand the control mecha-nism analytically [28–31]. A notable result has been re-ported recently in the context of refuting the alleged odd-number limitation [32, 33], believed to be a severe limita-tion of the delayed feedback control technique for almosta decade [34]. A way to overcome the odd-number lim-itation was proposed by Pyragas by using an additionalunstable degree of freedom in a feedback loop [35–37].The original Pyragas method was subsequently im-proved by introducing multiple delayed signals into thefeedback control force, both with commensurate or in-commensurate delay times [38–41]. Another improve-ment of the method was achieved by introducing a time-varying delay into the delayed argument of the controlforce which could be realized experimentally by changingsome characteristics of the delay line in a deterministic(periodic) or a stochastic fashion [42–44]. For example,in coupled laser systems the delay time could be modu-lated by periodically changing the distance of the lasersor external resonator, e.g. by piezoelectric modulation.For systems acting on a slower timescale, e.g. electronicor mechanic devices, one may take advantage of digi-tal delay lines instead, where the time-functional formof the delay modulation could be controlled by an ex-ternal clock frequency modulator [45]. In this way, onecould practically realize the variable-delay feedback con-troller and stabilize unstable steady states and periodicorbits over a much larger domain of control parameters.To keep the method noninvasive also in case of unstableperiodic orbits, the delay modulation should be chosenappropriately, e.g., in a form of a (non)periodic changeof the delay between multiples of the basic period of theuncontrolled unstable orbit [46].In this paper, we aim to provide a deeper analyticalinsight into the mechanism of stabilization of unstableequilibria by the delayed feedback method with a time-dependent delay time, and extend this concept by con-sidering the influence of a phase-dependent coupling, thecontrol loop latency, and multiple delay lines on the con-trollability, which is relevant in real experiments. A gen-eral description of the variable-delay feedback controlleris given in Sec. II, where it is shown that the methodcould be explored analytically in the domain of a high-frequency modulation of the delay. In this case, we usethe formalism of distributed delays [47, 48], where theaverage contribution of the time-varying delay is repre-sented by an integral kernel describing a particular delaydistribution. In Sec. III, we consider a phase-dependentcoupling of the control force, which is a particular non-diagonal generalization of the standard diagonal cou-pling scheme. As a model subject to control we use ageneric two-dimensional linear system that has an un-stable steady state of focus type. This model mimicsthe two-dimensional center manifold of a general non-linear system, capturing the essence of the dynamics inthe vicinity of its unstable fixed point. In Sec. IV, weinvestigate the influence of control-loop latency on theefficiency of the controller. In Sec. V, we extend ourmethod to the case of multiple delay feedback by con-sidering two independent delay lines with time-varyingdelays. We conclude our findings in Sec. VI.
II. VARIABLE-DELAY FEEDBACK CONTROL
We consider a general n -dimensional dynamical system˙ x = F ( x ), where F is the vector field describing the sys-tem’s dynamics, and x = x ( t ) ∈ R n is the state columnvector. The equilibria x ∗ i of the system are the solutionsof F ( x ∗ i ) = 0, and the stability of a particular equilibrium x ∗ is determined by the eigenvalues of the Jacobian ma-trix b A = ∂ x F ( x ) calculated at x = x ∗ . In the followingwe assume that x ∗ is an unstable equilibrium, having atleast one eigenvalue with a positive real part.Under a linear variable-delay feedback control, the original system is transformed into˙ x ( t ) = F ( x ( t )) + b K [ x ( t − τ ( t )) − x ( t )] , (1)where b K is an n × n feedback gain matrix, and τ ( t ) isthe time-dependent delay time. By choosing b K and τ ( t )appropriately, we aim to stabilize the unstable equilib-rium x ∗ . Without loss of generality, we assume that theequilibrium has already been moved to the origin by atranslation of axes, such that x ∗ = 0. We thus linearizethe controlled system in the neighborhood of the originto obtain ˙ x ( t ) = b Ax ( t ) + b K [ x ( t − τ ( t )) − x ( t )] . (2)We consider a deterministic modulation of the vari-able time-delay τ ( t ) in a fixed interval around a nominal(average) delay value τ τ ( t ) = τ + εf ( νt ) , (3)where f : R → [ − ,
1] is a 2 π -periodic function with zeromean, and ε and ν are additional variable-delay controlparameters denoting the amplitude and the (angular) fre-quency of the modulation, respectively. The stability ofthe origin can be inferred by numerically integrating thelinear variable-delay system (2) for different values of thecontrol parameters b K , τ , ε and ν , thus determining thedomains in the parameter space for which the stabiliza-tion can be realized.In the limiting case of a high-frequency modulation[43, 49], when the parameter ν becomes sufficiently largecompared to the uncontrolled system’s dynamics de-scribed by an intrinsic frequency ω , the linearized system(2) with a modulated delay has the same asymptotic sta-bility properties as the averaged linear distributed-delaysystem˙ x ( t ) = b Ax ( t ) + b K (cid:18)Z ∞ ρ ( θ ) x ( t − θ ) dθ − x ( t ) (cid:19) . (4)The probability density function ρ ( θ ), i.e. thedistributed-delay kernel, is defined in a way that ρ ( θ ) dθ gives the fraction of time for which τ ( t ) lies between θ and θ + dθ , satisfying ρ ( θ ) ≥ Z ∞ o ρ ( θ ) dθ = 1 . (5)In practice, the transition to the distributed-delay limitcase does not require very large modulation frequencies,and therefore variations of the delay in the experimentare a very practical way to create different types of delaydistributions.For a continuous modulation in an ε interval around τ in form of a sawtooth-wave (triangular) modulation ofthe delay, τ ( t ) = τ + ε (cid:20) (cid:18) νt π mod 1 (cid:19) − (cid:21) , (6) ρ ( θ ) is a constant in the interval of the delay variation.Under the probability normalization condition (5) we ob-tain a uniform distribution ρ ( θ ) = ( ε , θ ∈ [ τ − ε, τ + ε ] , , elsewhere, (7)which does not depend on the skewness of the sawtoothfunction. For a sine-wave modulation, τ ( t ) = τ + ε sin( νt ) , (8)the time interval dt during which the delay τ changes by dτ is given by dt = dτνε cos( νt ) = dτν p ε − ( τ − τ ) . (9)The fraction of time dt within a half-period π/ν of thedelay function τ is equivalent to the product ρ ( τ ) dτ . Interms of our previous notation, we get ρ ( θ ) = 1 π p ε − ( θ − τ ) . (10)A periodic change of τ ( t ) between two fixed values τ − ε and τ + ε for the same time duration, τ ( t ) = τ + ε sgn[sin( νt )] , (11)results in a square wave (rectangular) modulation with atwo-peak distribution ρ ( θ ) = ( δ ( θ − τ + ε ) + δ ( θ − τ − ε )) . (12)In an analogous way, we may obtain the correspondingdistributed delay kernels for other types of delay modu-lations.Applying the exponential ansatz x ( t ) ∼ exp(Λ t ) in Eq.(4), we obtain the characteristic equation for the eigen-values Λdet (cid:20) Λ b I − b A + (cid:18) − Z ∞ ρ ( θ ) e − Λ θ dθ (cid:19) b K (cid:21) = 0 , (13)where b I is the identity matrix. Since ρ ( θ ) is nonzero onlybetween τ − ε and τ + ε , we havedet h Λ b I − b A + (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) b K i = 0 . (14)The quantity χ (Λ , ε ) = Z + ε − ε ρ ( τ + θ ) e − Λ θ dθ (15)summarizes the effect of a given modulation, and for theabove modulation types reads χ (Λ , ε ) = sinh(Λ ε )Λ ε , sawtooth-wave ,I (Λ ε ) , sine-wave , cosh(Λ ε ) , square-wave , (16) τ ε ε t x FIG. 1: (Color online) Sample of a trajectory x ( t ) (black)in the vicinity of a steady state (red/gray horizontal line) ofa system subjected to variable-delay feedback control with afast sawtooth wave. The corresponding uniform delay distri-bution (7) creates a sliding average between t − τ − ε and t − τ + ε (area highlighted in blue/gray) which at time t re-sults in the reference signal (blue/gray, oscillating). The con-trol force is constructed from the difference signal as markedby the vertical black arrow. where I is the modified Bessel function of the first kind oforder 0. In the non-modulated case, χ (Λ , ≡
1, we havethe usual TDFC. For general distributed delay, the delayterm exp( − Λ τ ) is replaced in the characteristic equation(13) by the Laplace transform of the delay kernel ρ ( θ )[47].The function χ (Λ , ε ) also allows for a compact explana-tion of the mechanism of variable-delay feedback control(VDFC). The expressions listed in Eq. (16) have regions { Λ ε ∈ C } for which | χ | <
1, mostly for Re(Λ) ≪ Im(Λ).In particular, values of χ ≈ de-structive interference of the delay signals covered by thedistribution. Positive and negative phases of the spiraloscillation in the neighborhood of the steady state canceleach other out creating a reference term x ref ( t ) ≈ x ∗ fora simplified controller˙ x ( t ) = b Ax ( t ) + b K ( x ref ( t ) − x ( t )) ≈ b Ax ( t ) + b K ( x ∗ − x ( t )) . (17)The same idea can be regarded as the motivation for theoriginal TDFC method. The advantage of VDFC liesin an improved approximation of the target state by thedelay terms. Fig. 1 illustrates the control mechanismfor an almost ideal situation, in which the steady stateis approximated by the averaged delay signal. The in-stantaneous part of the coupling term stabilizes the sys-tem while pulling it towards the reference signal x ref ( t )which due to the small χ -value has a smaller amplitudethan the remaining instantaneous x ( t ). The resultingstabilization of the steady state shows a high robustnessagainst parameter detuning, because the simplified con-trol mechanism is insensitive to the phase relation be-tween actual signal and reference signal and also worksfor a wide range in the coupling gain b K . However, thedescribed scenario relies on a small χ -value, which in gen-eral is not trivial to obtain. The delay terms can almostnever be suppressed completely, so that only a rigorousconsideration of the characteristic equation (14) revealsthe full capability of the control method. Eq. (14) is tran-scendental with respect to Λ, possessing an infinite set ofcomplex solutions { Λ i } ∈ C . The origin can be stabilizedfor those values of the control parameters ( b K , τ , ε ) forwhich all the eigenvalues { Λ i } have negative real parts.The stability domain in the parameter space ( b K , τ , ε )can be calculated numerically given the modulation dis-tribution. Note that the control parameter ν is lost inthe transition to distributed delays.As in the usual TDFC and extended TDFC controlschemes, the presence of torsion is necessary for the pro-posed control method to be able to stabilize equilibria.Torsion means that in its unstable subspace the steadystate is only of the spiral type. To show this property,we consider the characteristic quasipolynomial H (Λ) =det h Λ b I − b A + (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) b K i . It is easily shownthat this quasipolynomial is positive for Λ → ∞ , and forΛ = 0 reads H (0) = det h − b A i = Q Nn =1 ( − e n ) , where e i are the eigenvalues of b A . If b A possesses an odd numberof positive real eigenvalues, then H (0) <
0, and there ex-ists at least one positive real root of H (Λ) = 0, meaningthat the fixed point cannot be stabilized by the proposedmethod. Although the analysis is done in the distributed-delay limit, numerical simulations show that this odd-number limitation persists also for low-frequency modu-lations of the delay τ ( t ).The results of this section will be exploited in thefollowing to perform a comparative analysis betweenthe variable-delay feedback control in a distributed-delay limit and the standard delayed feedback control by us-ing a simple normal form model as a representative ofa large class of nonlinear dynamical systems. Specifi-cally, we will investigate the effects of including a phase-dependent coupling matrix and the influence of nonzerolatency times in the control scheme. III. PHASE-DEPENDENT COUPLING
We will apply the variable-delay feedback control to anunstable steady state of focus type. This system repre-sents a generic model of an unstable steady state slightlyabove a Hopf bifurcation. In center manifold coordinates,the dynamics of the system is given by Eq. (2), where x = ( x, y ) T is the two-dimensional column vector, and b A is a 2 × b A = (cid:18) λ ω − ω λ (cid:19) (18)that describes the dynamics in the absence of control.Since the stability of the free-running system is deter-mined by the eigenvalues λ ± i ω of b A , choosing λ > ω = 0 we model an unstable focus at the origin.In this section, we investigate a phase-dependent cou-pling of the control force, which is relevant, e.g., in sta-bilization of laser devices where the optical phase occursas an additional degree of freedom [50, 51]. Such cou-plings have also been used in the context of refuting theodd-number limitation of delayed feedback control in au-tonomous systems [32, 33, 52], to anticipate chaos syn-chronization [53], and more recently, to control differentsynchronous oscillatory states of delay-coupled oscillatornetworks [54]. The phase-dependent coupling is realizedvia a rotational matrix b K containing a variable phase ϕ ,and it enters as an additional multiplicative factor to thecontrol force (cid:18) ˙ x ( t )˙ y ( t ) (cid:19) = (cid:18) λ ω − ω λ (cid:19) (cid:18) x ( t ) y ( t ) (cid:19) + K (cid:18) cos ϕ − sin ϕ sin ϕ cos ϕ (cid:19) (cid:18) x ( t − τ ( t )) − x ( t ) y ( t − τ ( t )) − y ( t ) (cid:19) , (19)where K is the feedback gain. In the distributed-delay limit, the system becomes (cid:18) ˙ x ( t )˙ y ( t ) (cid:19) = (cid:18) λ ω − ω λ (cid:19) (cid:18) x ( t ) y ( t ) (cid:19) + K (cid:18) cos ϕ − sin ϕ sin ϕ cos ϕ (cid:19) (cid:18) R ∞ ρ ( θ ) x ( t − θ ) dθ − x ( t ) R ∞ ρ ( θ ) y ( t − θ ) dθ − y ( t ) (cid:19) , (20)and the stability of the origin is determined by the roots Λ of the characteristic equation (cid:2) Λ − λ + K (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) cos ϕ (cid:3) + (cid:2) ω + K (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) sin ϕ (cid:3) = 0 . (21)This equation can be further simplified toΛ + Ke ∓ iϕ (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) = λ ± iω. (22) The control method is successful if there exists anonempty set in the four-dimensional parameter space( K, τ , ε, ϕ ) for which the real parts of all the character- FIG. 2: (Color online) Control domain in the ( ϕ, K ) plane fora sawtooth-wave modulation of the delay for different mod-ulation amplitudes: (a) ε = 0, (b) ε = 0 .
3, (c) ε = 0 . ε = 0 .
9. The nominal delay is fixed at τ = 1, which isoptimal value in the non-modulated case (TDFC). The pa-rameters of the free-running system are λ = 0 . ω = π .The grayscale (color code) depicts only those values of thecontrol parameters for which the largest real part of the com-plex eigenvalues Λ is negative, indicating a successful control.The depicted solid lines that enter into a description of thestability boundary are calculated from Eqs. (26) and (27),and the dashed lines from Eqs. (26) and (28). For clarity ofthe picture, we depict only a few branches of the boundarycurves. istic roots Λ are negative.The domain of control in the plane parametrized bythe feedback phase ϕ and the feedback gain K is givenin Figs. 2 and 3. The grayscale (color code) indicatesonly those control parameters ( ϕ, K ) for which the lead-ing eigenvalues have negative real parts, and the controlis more robust as these values are more negative. Panels(a)-(d) in Fig. 2 depict control domains for a sawtooth-wave modulation corresponding to a fixed nominal delay τ = 1 and different modulation amplitudes ε : (a) ε = 0,(b) ε = 0 .
3, (c) ε = 0 . ε = 0 .
9. The parametersof the unstable focus are fixed at λ = 0 . ω = π .The chosen nominal delay τ = 1 is an optimal valuein the non-modulated control case and ϕ = 0 (diagonalcoupling). From the corresponding control domains, it isobserved that an increase of the modulation amplitude ε leads to an enlargement of the control domain in thedirection of the positive ϕ -axis, thus destroying the sym-metry of the stability island centered at ϕ = 0. A similarresult is observed for a sine-wave modulation in panels(a)-(b) in Fig. 3, corresponding to ε = 0 . ε = 0 . ε is increased from FIG. 3: (Color online) Same as Fig. 2 for a sine-wave mod-ulation of the delay: (a) ε = 0 .
5, (b) ε = 0 .
9, and for asquare-wave modulation: (c) ε = 0 .
5, (d) ε = 0 .
9. Otherparameters are the same as in Fig. 2. zero, until ε reaches a critical value after which the do-main decreases, resulting in an instability stripe at ϕ = 0when ε = 1. This non-monotonic change of the stabilitydomain is typical for a square-wave modulation, and ithas been observed in other circumstances.In order to understand the mechanisms behind thesenumerical results, we first recall the intuitive argumentabout destructive interference and the simplified con-troller from the previous section. From Eq. (20) we di-rectly see that the stabilizing diagonal elements of the in-stantaneous control term are weighted with cos ϕ . Sincethe onset of stability K min is dominated by the com-petition between λ > − K cos ϕ <
0, we ex-pect K min ∝ λ/ cos ϕ , see also the detailed reasoning inEq. (30) below. For | ϕ | > π/ τ this type of support fails and | ϕ | = π/ ω were zero, the control term wouldhave no preference for positive or negative values of ϕ because the change of ϕ → − ϕ would be equivalent tomirroring the complete system at one coordinate axis,e.g. x → − x . A non-zero value of the spiral frequency ω breaks the symmetry and leads to different resultingfrequencies Ω = Im(Λ) under control for different signsof ϕ . From the interference mechanism one can conclude,that at least for the continuous distributions (7) and (10) FIG. 4: (Color online) Characteristic exponents Λ in thecomplex plane parametrized by the feedback gain K . Left: ϕ = −
1, right: ϕ = +1. Depicted are five leading eigenvaluecurves for K ∈ [0 , K = 0 and the black dots correspond to K = 6.The eigenmodes created by the variable-delay control origi-nate from Re(Λ) = −∞ . Other parameters as in Fig. 2c. high frequencies are more favorable for successful control,because a low χ -value corresponds to many oscillationperiods covered by the distribution. Indeed, if we havea closer look at the parameter constellation underlyingFig. 2, we see that for negative ϕ the imaginary part ofthe most unstable mode tends to vanish, thus inhibit-ing the interference effect and leading to control failure.In contrast, for positive ϕ the frequency increases withincreasing K , thus improving control efficiency. A para-metric plot of different Λ i in the complex plane revealsthis property, see Fig. 4.A detailed analytical investigation of the control do-main can be done by using the characteristic equation(22) and noticing that the transition from instability tostability occurs at the boundaries of the control domain,where the leading eigenvalues Λ are purely imaginary. Setting Λ = i Ω in Eq. (22), and separating real andimaginary parts results in a system of two real-valuedequations K cos ϕ − Kχ ( i Ω , ε ) cos(Ω τ ± ϕ ) = λ, (23)Ω ∓ K sin ϕ + K sin(Ω τ ± ϕ ) χ ( i Ω , ε ) = ± ω, (24)which is an implicit parametric representation of theboundaries of the control domain in the parameter space( K, τ , ε, ϕ ), in which the eigenfrequency Ω has the roleof a parametrization variable. We have taken into ac-count that for the above mentioned distribution kernels, ρ ( θ ) is even around the nominal delay τ , meaning that χ ( i Ω , ε ) is a real-valued function χ ( i Ω , ε ) = sin(Ω ε )Ω ε , sawtooth-wave ,J (Ω ε ) , sine-wave , cos (Ω ε ) , square-wave , (25)where J is the Bessel function of the first kind of or-der 0. Considering the positivity of the nominal delay τ and the feedback gain K , as well as the positivity of theparameters ϕ, λ and ω used in the current analysis, elim-ination of the phase parameter ϕ from the last system oftranscendental equations leads to an expression for K interms of the eigenfrequency Ω K (Ω) = s λ + ( ω − Ω) sin (Ω τ ) + [ χ ( i Ω , ε ) − cos(Ω τ )] . (26)Taking into account the multivalued properties of the arc-sine function, one obtains in a similar manner the analyt-ical expressions for the phase parameter ϕ in dependenceon Ω ϕ (Ω) = arcsin (cid:18) K (Ω) χ ( i Ω , ε )[ λ sin(Ω τ ) + ( ω − Ω) cos(Ω τ )] − ( ω − Ω) λ + ( ω − Ω) (cid:19) + 2 nπ (27) ϕ (Ω) = − arcsin (cid:18) K (Ω) χ ( i Ω , ε )[ λ sin(Ω τ ) + ( ω − Ω) cos(Ω τ )] − ( ω − Ω) λ + ( ω − Ω) (cid:19) + (2 n + 1) π, (28)where n is a non-negative integer. Equations (26)–(28)describe the boundary of the control domain in Figs.2 and 3, where the solid and dashed lines representbranches of ϕ (Ω) and ϕ (Ω), respectively.In the case of zero feedback phase ( ϕ = 0), the phase-dependent feedback force is reduced to a diagonal cou-pling. In the absence of modulation (TDFC), the domainof control in the plane spanned by the feedback gain andthe nominal delay consists of several distinct stability is-lands centered at odd values of τ , separated by regions corresponding to even τ for which the stabilization failsfor any K . This observation can be inferred from Eqs.(23)-(24) setting ϕ = 0 and ε = 0. When Ω τ = (2 n +1) π ,then K min = λ/ ω . On the other hand, whenΩ τ = 2 nπ , the control fails for any finite value of K .A detailed analysis of the control boundaries in this case(TDFC) is provided in Ref. [29]. By the same reasoning,in the presence of modulation ( ε >
0) and zero phase,the corresponding minimum feedback gain values are K min ( ε ) = λ [1 + χ ( iω, ε )] , τ = (2 n + 1) π/ω,λ [1 − χ ( iω, ε )] , τ = 2 nπ/ω, (29)leading to a reconfiguration of the control domain de-pending on the type of the delay modulation [42]. Specif-ically, it is observed that the instability stripes at evenvalues of τ cease to exist, and the stability islands arestarting to interconnect as soon as ε >
0. Although for asawtooth and sine-wave modulation the enlargement ofthe control domain with increasing ε is monotonic, anoscillatory behavior is observed for a square-wave modu-lation due to the form of the χ function (i.e. cos( ωε )) inthe denominator of the expressions for K min . In the lat-ter case, the stability islands centered at τ = (2 n +1) π/ω are first enlarged and eventually interconnected with in-creasing ε , up to some ε value after which the controldomain shrinks, gradually collapsing into several stabil-ity islands centered at τ equal to even multiples of π/ω .For ϕ = 0 and a variable-delay control force ( ε > τ cannot be deduced from Eqs. (23)-(24) by followingthe previous lines of deduction. Nevertheless, we canmake approximate predictions for small values of ϕ , whenΩ τ is an integer multiple of π , and we have Ω ≈ ω fromEq. (24). Hence, in the regime of small values of ϕ , theminimum feedback gain is K min ( ε, ϕ ) = λ [1 + χ ( iω, ε )] cos ϕ , τ ≈ (2 n + 1) π/ω,λ [1 − χ ( iω, ε )] cos ϕ , τ ≈ nπ/ω. (30)The resulting expressions for K min show that in theabsence of any delay modulation and for small ϕ , theminimum feedback gain at the ”optimal” delay values τ = (2 n + 1) π/ω is K min = λ/ (2 cos ϕ ), containing thedependence on ϕ via cos ϕ in the denominator. Conse-quently, the principal stability island in ( ϕ, K ) paramet-ric plane is centered at ϕ = 0, for which K min = λ (seepanel (a) in Fig. 2). On the other hand, at τ = 2 nπ/ω ,the stabilization cannot be achieved for any K , regardlessof the value of the phase parameter ϕ . For a non-zeromodulation amplitude ε at small ϕ , the values of K min de-pend on the type and the amplitude of the delay modula-tion. Although the form of χ ( iω, ε ) for sawtooth and sine-wave modulations is such that the control in the variable-delay case is generally possible for all nominal delays τ ≥ ε , this is not the case for a square-wave modulation.Namely, in the latter case, χ ( iω, ε ) = cos( ωε ), giving anon-monotonic behavior of K min with increasing modu-lation amplitude ε . Specifically, for ε = (2 n + 1) π/ω , thecontrol fails at τ = (2 n + 1) π/ω for any K , but it isoptimal at τ = 2 nπ/ω .To further investigate the effects of the modulation am-plitude ε on the control efficiency, we have calculated the FIG. 5: (Color online) Control domain in the (
K, ε ) plane fora sawtooth-wave modulation of the delay for different valuesof the feedback phase ϕ : (a) ϕ = 0, (b) ϕ = π/
8, (c) ϕ = π/ ϕ = 3 π/
8. The nominal delay is fixed at τ = 1. Theparameters of the free-running system are as in Fig. 2. stability domains in the parametric plane ( K, ε ) for dif-ferent modulation types and different values of the phaseparameter ϕ . The results are depicted in Figs. 5, 6 and7 for sawtooth-, sine- and square-wave modulation, re-spectively. Different panels in each figure correspond todifferent phases: (a) ϕ = 0, (b) ϕ = π/
8, (c) ϕ = π/ ϕ = 3 π/
8. The parameters of the unperturbed sys-tem are λ = 0 . ω = π , and the nominal delay isfixed at τ = 1 as before. As the modulation amplitudeincreases from zero, a considerable reconfiguration of thestability domain is observed, depending on the type ofthe delay modulation. Increasing the modulation am-plitude enlarges the interval of K for successful control,although the enlargement is not necessarily monotonic,and the domain may consist of several disconnected in-tervals. The improvement of the delayed feedback con-troller from including time-varying delay is evident fromthe diagrams. For example, for a phase parameter value ϕ = 3 π/ K , the variable-delayfeedback method is able to stabilize the fixed point forcertain values of the modulation amplitude ε > K, ε ) plane can be obtained from Eqs. (23) and(24) K (Ω) = λ sin(Ω τ + ϕ ) + ( ω − Ω) cos(Ω τ + ϕ )sin(Ω τ ) , (31) χ ( i Ω , ε ) = λ sin ϕ + ( ω − Ω) cos ϕλ sin(Ω τ + ϕ ) + ( ω − Ω) cos(Ω τ + ϕ ) . (32)In the second equation, ε is implicitly contained in FIG. 6: (Color online) The corresponding control domain inthe (
K, ε ) plane for a sine-wave modulation of the delay fordifferent values of the feedback phase ϕ : (a) ϕ = 0, (b) ϕ = π/
8, (c) ϕ = π/ ϕ = 3 π/
8. The other parametersare as in Fig. 5.FIG. 7: (Color online) The control domain in the (
K, ε ) planefor a square-wave modulation of the delay for different valuesof the feedback phase ϕ : (a) ϕ = 0, (b) ϕ = π/
8, (c) ϕ = π/ ϕ = 3 π/
8. The other parameters are as in Fig. 5. χ ( i Ω , ε ), which is a real function in the delay modula-tion cases considered above (see Eq. 25). In general itis not possible to invert the function χ analytically. Oneshould treat the last system as implicit parametric rep-resentation of the control boundaries and seek for thestability curves numerically. The boundaries of the cor-responding control domains are given in Figs. 5–7 by thesolid line. IV. CONTROL-LOOP LATENCY
In an experimental realization of the control method,one has to take into account the latency of the feedbackcircuit due to the time required for the generation of thefeedback control signal and its reinjection into the sys-tem. In laser systems, the latency is associated with thetime the light takes to traverse the distance between thelaser and the Fabry-Perot controller. It has been shownthat latency time always acts destructively upon the con-trol domains, reducing the effectiveness of the controller[29, 30, 51, 55, 56].In this section, we investigate the effects of latencyon the variable-delay feedback control in the distributed-delay limit. The latency time δ enters as an additionalconstant time-delay in the control force, and the systemnow reads (cid:18) ˙ x ( t )˙ y ( t ) (cid:19) = (cid:18) λ ω − ω λ (cid:19) (cid:18) x ( t ) y ( t ) (cid:19) + K (cid:18) x ( t − δ − τ ( t )) − x ( t − δ ) y ( t − δ − τ ( t )) − y ( t − δ ) (cid:19) . (33)In the distributed-delay limit, the system becomes (cid:18) ˙ x ( t )˙ y ( t ) (cid:19) = (cid:18) λ ω − ω λ (cid:19) (cid:18) x ( t ) y ( t ) (cid:19) + K (cid:18) R ∞ ρ ( θ ) x ( t − δ − θ ) dθ − x ( t − δ ) R ∞ ρ ( θ ) y ( t − δ − θ ) dθ − y ( t − δ ) (cid:19) , (34)leading to the characteristic equation (cid:2) Λ − λ + Ke − Λ δ (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1)(cid:3) + ω = 0 , (35) that can be further simplified intoΛ + Ke − Λ δ (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) = λ ± iω. (36) FIG. 8: (Color online) Control domain in the ( τ , K ) plane fora sawtooth-wave modulation of the delay for different valuesof the latency time: (a) δ = 0 .
1, (b) δ = 0 .
2, (c) δ = 0 . δ = 0 .
4. The modulation amplitude is fixed at ε = 0 . ϕ = 0. Upon separating the real and imaginary parts, one ob-tains a system of two real-valued equations. We againconsider the control boundary given by Λ = i Ω which isgiven by the implicit parametric representation K cos(Ω δ ) − K cos(Ω( δ + τ )) χ ( i Ω , ε ) = λ, (37)Ω − K sin(Ω δ ) + K sin(Ω( δ + τ )) χ ( i Ω , ε ) = ± ω. (38) Figure 8 depicts the domain of control in the ( τ , K )parametric plane for a sawtooth-wave modulation at afixed modulation amplitude ε = 0 . δ : (a) δ = 0 .
1, (b) δ = 0 .
2, (c) δ = 0 .
3, (d) δ = 0 .
4. The control domains for a sine-wave modula-tion corresponding to δ = 0 . δ = 0 . λ = 0 . ω = π as before. For a zero latency time,the control domain fills the whole depicted range of theparametric plane down to some minimum value of K ,which for τ equal to an integer multiple of π/ω is ex-plicitly given by Eq. (29). It is observed that increasinglatency time reduces the area and the robustness of thecontrol domain. In analogy to Sec. II, one can deriveapproximate expressions for K min for small values of thelatency δK min ( ε, δ ) = λ [1 + χ ( iω, ε )] cos( ωδ ) , τ = (2 n + 1) π/ω,λ [1 − χ ( iω, ε )] cos( ωδ ) , τ = 2 nπ/ω. (39)In parallel to the conclusions of the previous section, in-clusion of a variable time-delay in the feedback controlforce generally enlarges the control domain, making thecontrol possible even for those values of τ for which thecontrol always fails at any K in the non-modulated case.To describe the boundaries of the control domain in the( τ ,K) plane analytically, we algebraically manipulate thesystem (37)-(38) and find two families of branches of so-lutions for K and τ parametrized by the eigenfrequencyΩ K (Ω) = λ cos(Ω δ ) − ( ω − Ω) sin(Ω δ ) ± p [ λ cos(Ω δ ) − ( ω − Ω) sin(Ω δ )] − (1 − χ ( i Ω , ε ))[ λ + ( ω − Ω) ]1 − χ ( i Ω , ε ) , (40) τ , (Ω) = 1Ω (cid:20) arcsin (cid:18) λ sin(Ω δ ) + ( ω − Ω) cos(Ω δ ) K (Ω) χ ( i Ω , ε ) (cid:19) + 2 nπ (cid:21) , (41) τ , (Ω) = 1Ω (cid:20) − arcsin (cid:18) λ sin(Ω δ ) + ( ω − Ω) cos(Ω δ ) K (Ω) χ ( i Ω , ε ) (cid:19) + (2 n + 1) π (cid:21) , (42)where n is a non-negative integer that takes care of thedifferent branches due to the multivalued arcsine func-tion involved in the boundary description. The controldomain boundaries in Figs. 8 and 9 are represented para-metrically by Eqs. (40)–(42), where the branches relatedto τ , are given by the solid lines, and the branches of τ , with dashed lines.The influence of the modulation amplitude ε on the stability of the fixed point can be deduced by observ-ing the control domains in the ( K, ε ) parameter plane.Figures 10, 11 and 12 depict the control domains forsawtooth-, sine- and square-wave modulation, respec-tively, for a fixed nominal delay τ = 1 and increasing val-ues of the latency parameter δ : (a) δ = 0 .
1, (b) δ = 0 . δ = 0 .
3, (d) δ = 0 .
4. The parameters of the unstablefocus were set at λ = 0 . ω = π as previously. It is0 FIG. 9: (Color online) Same as Fig. 8 for a sine-wave modu-lation of the delay for fixed ε = 0 . δ = 0 .
2, (b) δ = 0 .
4, and for a square-wavemodulation: (c) δ = 0 .
2, (d) δ = 0 .
4. Other parameters arethe same as in Fig. 8. observed that increasing the modulation amplitude ε fora constant latency value generally leads to an enlarge-ment of the K interval of successful control, dependingon the type of the modulation. For a sufficiently largelatency time when the Pyragas control fails (e.g, δ = 0 . ε >
0. Nev-ertheless, the control domain depends on the modulationtype, and it is seen that in this case it is much less pro-nounced for a sawtooth-wave modulation than for theother two delay-modulation types.In a similar fashion, from Eqs. (37) and (38) oneobtains the parametric representation of the stabilityboundaries in the (
K, ε ) plane K (Ω) = λ sin Ω( δ + τ ) + ( ω − Ω) cos Ω( δ + τ )sin(Ω τ ) , (43) χ ( i Ω , ε ) = λ sin(Ω δ ) + ( ω − Ω) cos(Ω δ ) λ sin Ω( δ + τ ) + ( ω − Ω) cos Ω( δ + τ ) , (44)in which the dependence of ε on the eigenfrequency Ωenters implicitly in the second equation via the function χ . The calculated boundaries of the control domain arethe solid lines depicted in Figs. 10–12.Overall, we see a massive destruction of the controldomains with increasing latency time. Compared to theeffect of a phase rotation in the previous section, thepresence of control-loop latency tends to spoil the con-trol mechanism radically. Besides the destructive inter-ference effect, on which variable-delay feedback controlmainly relies, an instantaneous feedback provides the FIG. 10: (Color online) Control domain in the (
K, ε ) plane fora sawtooth-wave modulation of the delay for different valuesof the latency time: (a) δ = 0 .
1, (b) δ = 0 .
2, (c) δ = 0 . δ = 0 .
4. The nominal time delay is fixed at τ = 1. Theparameters of the free-running system are as in Fig. 2.FIG. 11: (Color online) Control domain in the ( K, ε ) planefor a sine-wave modulation of the delay for different values ofthe latency time: (a) δ = 0 .
1, (b) δ = 0 .
2, (c) δ = 0 . δ = 0 .
4. The other parameters are as in Fig. 10. main source of stability which normally is reflected inan extensive coverage of the parameter space with so-lutions of successful control. If this dissipation term isreplaced by a much less effective term due to latency, thecontrol performance is consequently lost. Thus in anyexperimental application of variable-delay feedback con-trol one should split the control term if possible, so thatthe instantaneous part is implemented by a direct mod-ification of the system to be controlled, and the variabledelay part can be realized by separate devices.1
FIG. 12: (Color online) Control domain in the (
K, ε ) planefor a square-wave modulation of the delay for different valuesof the latency time: (a) δ = 0 .
1, (b) δ = 0 .
2, (c) δ = 0 . δ = 0 .
4. The other parameters are as in Fig. 10.
V. MULTIPLE DELAY FEEDBACK CONTROLWITH VARIABLE TIME-DELAYS
In order to improve the control of unstable steadystates, several extensions of the Pyragas method wereproposed by involving multiple delay feedback terms inthe control force [38–41]. A very efficient control schemeof this type was introduced by Ahlborn and Parlitz byutilizing a feedback force constructed from two or moreindependent Pyragas delayed feedback controllers withincommensurate delay times applied simultaneously inthe control circuit [40, 41]. A key result of this multiple-delay extension was a successful stabilization of chaotic intensity fluctuations of a frequency doubled ND-dopedyttrium aluminum garnet (ND:YAG) laser at higherpump rates, whose control was not achievable with a sin-gle delay controller.In this section, we show that this multiple delay feed-back control (MDFC) can be further improved by includ-ing time-varying delays in the associated feedback terms.As we have shown in the previous sections, the term im-provement refers in the first instance to an extension ofthe parameter space of successful control. Under exper-imental conditions this feature is favorable if control pa-rameters are drifting, cannot be adjusted precisely or theproperties of the unstable equilibrium are unknown sothat the optimal parameters of the constant delay controlforces cannot be determined. Concerning the robustnessof the control, we have observed for an optimal choice ofmodulation that the resulting maximum stability expo-nent in the stable domain is rather uniform. But sincethe complete response to perturbations of the controlledsystem is formed by the whole spectrum of exponents,the overall robustness can vary significantly more than re-flected only by the maximum exponent. In particular, wehave observed [45] that the excitation from the stabilizedfixed point by additive noise is decreasing with increasingfeedback strength in an optimal VDFC scheme althoughthe maximum exponent does not vary much. This is an-other typical feature of VDFC which we regard as a majorimprovement compared to non-modulated methods. Inthe following analysis of modulated MDFC we restrictourselves to the maximum stability exponent in order toallow for comparison with the previous sections. For theclarity of the presentation, we consider the simplest real-ization of MDFC consisting of two Pyragas-type delayedfeedback lines with variable time delays, but the discus-sion can be straightforwardly generalized to more thantwo delay lines. The linear normal-form system reads (cid:18) ˙ x ( t )˙ y ( t ) (cid:19) = (cid:18) λ ω − ω λ (cid:19) (cid:18) x ( t ) y ( t ) (cid:19) + K (cid:18) x ( t − τ ( t )) − x ( t ) y ( t − τ ( t )) − y ( t ) (cid:19) + K (cid:18) x ( t − τ ( t )) − x ( t ) y ( t − τ ( t )) − y ( t ) (cid:19) , (45)where τ ( t ) and τ ( t ) are two different, time-dependent, 2 π -periodic delay functions: τ ( t ) = τ + ε f ( ν t ) , (46) τ ( t ) = τ + ε f ( ν t ) . (47)The nominal delays are denoted by τ and τ , ε and ε are the corresponding modulation amplitudes, ν and ν are the modulation frequencies, and f and f are the modulation functions. For high-frequency modulations of thedelays, the system is in the distributed delay regime: (cid:18) ˙ x ( t )˙ y ( t ) (cid:19) = (cid:18) λ ω − ω λ (cid:19) (cid:18) x ( t ) y ( t ) (cid:19) + K (cid:18) R ∞ ρ ( θ ) x ( t − θ ) dθ − x ( t ) R ∞ ρ ( θ ) y ( t − θ ) dθ − y ( t ) (cid:19) + K (cid:18) R ∞ ρ ( θ ) x ( t − θ ) dθ − x ( t ) R ∞ ρ ( θ ) y ( t − θ ) dθ − y ( t ) (cid:19) . (48)Using the exponential ansatz, we obtain the characteris- tic equation (cid:2) Λ − λ + K (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) + K (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1)(cid:3) + ω = 0 , (49)2where χ (Λ , ε ) = Z ε − ε ρ ( τ + θ ) e − Λ θ dθ,χ (Λ , ε ) = Z ε − ε ρ ( τ + θ ) e − Λ θ dθ (50)are complex functions corresponding to different modu-lation types. Equation (49) can be re-cast in the simpleform Λ+ K (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) + K (cid:0) − e − Λ τ χ (Λ , ε ) (cid:1) = λ ± iω. (51)The parametric representation of the stability boundaryis obtained by looking for solutions of Eq. (51) in theform Λ = i Ω and separating real and imaginary parts, K [1 − cos(Ω τ ) χ ( i Ω , ε )]+ K [1 − cos(Ω τ ) χ ( i Ω , ε )] = λ, (52) K sin(Ω τ ) χ ( i Ω , ε )+ K sin(Ω τ ) χ ( i Ω , ε ) = ± ω − Ω . (53)The control parameter space is now six-dimensional( K , K , τ , τ , ε , ε ), but for experimental purposes itcould be reduced by matching similar types of controlparameters (e.g. K = K , or ε = ε ).To gain insight into the domain structure of the mul-tiple variable-delay feedback control, we numerically an-alyze Eq. (51) in the parametric plane spanned by thetwo nominal delays τ and τ . In Fig. 13 we depictthe resulting stability diagrams at different control pa-rameter values obtained when the time delays in bothfeedback terms are modulated with sawtooth-waves ofequal amplitudes (i.e. ε = ε at each panel). The pa-rameters of the unstable focus are taken as λ = 0 . ω = π throughout this section. The stability area isgiven in graytones (color online), corresponding to thosevalues of the control parameters for which the fixed pointcontrol could be achieved. Panels (a)–(d) correspondto increasing values of ε , for a symmetrical choice ofthe feedback gains K = K = 0 .
1. Panel (a) is theresulting stability domain without any delay variation( ε , = 0), i.e. MDFC with constant delays. The sta-bility domain is symmetrical with respect to the diag-onal τ = τ , and it is filled with oval-shaped insta-bility islands (white regions) encompassing the points( τ , τ ) = (2 nπ/ω, mπ/ω ), where n and m are non-negative integers. At these particular time delays, thefixed point is unstable for any values of the feedbackgains K and K . The result could be shown analyti-cally from Eqs. (52)–(53) by setting Ω τ = 2 nπ andΩ τ = 2 mπ . In the same manner, we derive the condi-tions for the minimum feedback gains at the nominal de-lay values between these instability points. Specifically,at ( τ , τ ) = ((2 n + 1) π/ω, mπ/ω ) the minimum feed-back gain is K min1 = λ/
2; at ( τ , τ ) = (2 nπ/ω, (2 m +1) π/ω ) the minimum feedback gain is K min2 = λ/
2; at ( τ , τ ) = ((2 n + 1) π/ω, (2 m + 1) π/ω ) the minimumfeedback gains satisfy condition K min1 + K min2 = λ/
2. Inthe latter case, by tuning the value of one of the feed-back gain parameters, a successful stabilization could beachieved at negative values of the second feedback gainparameter.The change in the stability diagrams as the amplitudes ε , increase from zero becomes evident from Fig. 13.Panels (b)–(d) show a monotonic increase of the controldomain area for increasing modulation amplitudes: (b) ε , = 0 .
25, (c) ε , = 0 .
5, (d) ε , = 1. As a conse-quence, the instability islands become smaller, and even-tually disappear at higher values of ε , (panel d). This isconfirmed analytically by deriving the minimum feedbackgains in the presence of delay variations when nominaldelay values are integer multiples of π/ω . By the samearguments as in the previous paragraph, we obtain: K min1 [1 − χ ( iω, ε )] + K min2 [1 − χ ( iω, ε )] = λ (54)at ( τ , τ ) = (2 nπ/ω, mπ/ω ), K min1 [1 + χ ( iω, ε )] + K min2 [1 − χ ( iω, ε )] = λ (55)at ( τ , τ ) = ((2 n + 1) π/ω, mπ/ω ), K min1 [1 − χ ( iω, ε )] + K min2 [1 + χ ( iω, ε )] = λ (56)at ( τ , τ ) = (2 nπ/ω, (2 m + 1) π/ω ), and K min1 [1 + χ ( iω, ε )] + K min2 [1 + χ ( iω, ε )] = λ (57)at ( τ , τ ) = ((2 n + 1) π/ω, (2 m + 1) π/ω ). It is seen thatwhile the multiple delay feedback control with constantdelays was unsuccessful at ( τ , τ ) = (2 nπ/ω, mπ/ω )for any values of the feedback gain parameters K , , in-clusion of variable delays lifts this restriction, making thecontrol possible with minimum gain parameters given byEq. (54).In panels (e)–(h) of Fig. 13 we show the correspond-ing increase of the control domain area for K = 0 .
05 and K = 0 .
2. This asymmetric choice of the feedback gainsresults in asymmetry of the stability domain with respectto the diagonal τ = τ , leading to horizontal enlarge-ment of the instability islands in the absence of delaymodulation (MDFC, panel a), which eventually connectinto horizontal instability stripes at τ = 2 mπ/ω forlower values of the gain parameter K . By including vari-able delays, the stability domain increases monotonicallywith increasing ε , , and the instability islands graduallyshrink, eventually disappearing at large ε , .A similar monotonic increase of the stability regionswith increasing modulation amplitudes in both symmet-ric and asymmetric choice of the feedback gain param-eters is also observed for sine-wave modulations of thedelays, and the results are not shown for compactnessof the presentation. For a square-wave modulation (Fig.14), we observe a non-monotonic behavior of the stabil-ity area, increasing faster with ε , , attaining its maxi-mum at ε , = 0 . FIG. 13: (Color online) Control domain in the ( τ , τ ) plane for multiple-delay feedback control (MDFC) of unstable focuswith two delay lines and time-varying delays. The delays τ ( t ) and τ ( t ) in both lines are varied with sawtooth-wave modulationswith matching modulation amplitudes. The feedback gain parameters are: (a)–(d) K = K = 0 .
1; (e)–(h) K = 0 .
05 and K = 0 .
2. The delay modulation amplitudes are: (a,e) ε , = 0 (MDFC with fixed delays); (b,f) ε , = 0 .
25; (c,g) ε , = 0 . ε , = 1. The parameters of the unstable focus are λ = 0 . ω = π . The boundary curves are given in parametricform by Eqs. (61) and (62). again to the same domain structure as in the unmodu-lated case ( ε , = 1, panels b and d), but now with theinstability islands centered at the nominal delay valuesbeing odd multiples of π/ω (compare with panels a ande in Fig. 13). This oscillatory switching of the insta-bility islands between ( τ , τ ) = (2 nπ/ω, mπ/ω ) and( τ , τ ) = ((2 n + 1) π/ω, (2 m + 1) π/ω ) continues as themodulation amplitudes ε , are further increased, and itis confirmed analytically by Eqs. (54)–(57).It is also instructive to investigate the stability dia-grams when the corresponding delays in the two feed-back control terms are modulated with different modu-lation types. A general conclusion can be extracted fromthe sample of the results provided in Fig. 15. The di-agrams are calculated for a sawtooth-wave variation of τ ( t ) and a square-wave variation of τ ( t ) at equal am-plitudes ( ε = ε ) and asymmetric gains: K = 0 . K = 0 .
05 in panels (a)–(b), and K = 0 .
05 and K = 0 . ε , = 0 . ε , = 1. The over-all behavior of the control domain with increasing ε , is mainly determined by the dominant modulation type,i.e., the modulation of the feedback with higher gain, andit is more pronounced as the difference between the re-spective feedback gain values becomes larger. By gradu-ally increasing the modulation amplitudes, the domain isenlarged monotonically in panels (a)–(b) where the feed-back term with a sawtooth-wave modulation dominatesover the square-wave term ( K > K ). The reversedcase scenario ( K > K ) in panels (c)–(d) shows a non- monotonic behavior typical for a square-wave modula-tion by the appearance of instability stripes oscillatingbetween the values of τ being even and odd multiplesof π/ω .To derive the equations for the boundary curves inFigs. (13)–(15), we follow the approach in Ref. [57] andrewrite the characteristic Eq. (51) as1 + a (Λ) e − Λ τ + b (Λ) e − Λ τ = 0 , (58)where a (Λ) and b (Λ) are given by a (Λ) = − K χ (Λ , ε )Λ − λ ∓ iω + K + K , (59) b (Λ) = − K χ (Λ , ε )Λ − λ ∓ iω + K + K . (60)At the control boundary (Λ = i Ω) the three terms in Eq.(58) can be considered as three vectors in the complexplane, with the corresponding magnitudes 1, | a ( i Ω) | and | b ( i Ω) | . According to Eq. (58), the sum of these vec-tors is a zero vector, and from the triangle formed bythe vectors it is straightforward to obtain the parametric4 FIG. 14: (Color online) Control domain in the ( τ , τ ) planefor square-wave delay modulations in both delay lines at sameamplitudes. The feedback gain parameters: (a)–(b) K = K = 0 .
1; (c)–(d) K = 0 .
05 and K = 0 .
2. The delaymodulation amplitudes: (a,c) ε , = 0 .
5; (b,d) ε , = 1. Theother parameters are as in Fig. 13.FIG. 15: (Color online) Control domain in the ( τ , τ ) plane.The delay τ ( t ) is modulated with a sawtooth-wave, and τ ( t )with a square-wave, both with the same modulation ampli-tude. The feedback gain parameters: (a)–(b) K = 0 . K = 0 .
05; (c)–(d) K = 0 .
05 and K = 0 .
2. The delay mod-ulation amplitudes: (a,c) ε , = 0 .
5; (b,d) ε , = 1. The otherparameters are as in Fig. 13. representation of τ and τ on the eigenfrequency Ω: τ (Ω) = Arg [ a ( i Ω)] + (2 u − π ± α Ω ≥ ,u = u ± , u ± + 1 , u ± + 2 . . . , (61) τ (Ω) = Arg [ b ( i Ω] + (2 v − π ∓ α Ω ≥ ,v = v ± , v ± + 1 , v ± + 2 . . . , (62)where u ± and v ± are the smallest possible integers suchthat the corresponding values of τ and τ are all non-negative, and α , α ∈ [0 , π ] are the internal angles of thetriangle formed by the vectors, calculated from the lawof cosines as: α = Arccos (cid:18) | a ( i Ω) | − | b ( i Ω) | | a ( i Ω) | (cid:19) , (63) α = Arccos (cid:18) | b ( i Ω) | − | a ( i Ω) | | b ( i Ω) | (cid:19) . (64)To gain further insight into the superiority of MDFCwith time-varying delays with respect to the fixed de-lays realization, we have numerically calculated the sta-bility diagrams in the parametric plane of the nominaldelay τ and the feedback gain K of the first feed-back line at different ε , , by fixing the nominal delay τ and the feedback gain K of the second feedback line.For the value τ we choose an even multiple of π/ω atwhich the fixed point control by MDFC with constantdelays always fails if τ is also an even multiple of π/ω .In Fig. 16 we show the corresponding control domainsfor τ = 2 π/ω = 2 and K = 0 . ε , =0 (MDFC), (b) ε , =0.25, (c) ε , =0.5,(d) ε , =1. For fixed delays (MDFC case, panel a), thecontrol domain is consisted of isolated stability islandsencompassing τ = (2 n + 1) π/ω , disconnected by in-stability stripes at τ = 2 nπ/ω at which control fails atany K , . As the modulation amplitudes increase (panelsb-d), the stability islands gradually expand into a con-nected stability region, making the control possible atany τ . This monotonic expansion of the stability do-main is sustained in the case of a sine-wave modulation,but not for a square-wave modulation, as expected fromthe earlier analysis. In the latter case, instability stripesare disappearing and re-appearing again with increasing ε , in an oscillatory manner at even multiples of π/ω along τ axis. When the delays are modulated withdifferent modulation types, the behavior of the domainstructure with increasing ε , is determined by the dom-inant feedback term.The parametric representation of the stability bound-aries in Fig. 16 parametrized by the eigenfrequency Ω5 FIG. 16: (Color online) Control domain in the ( K , τ ) planefor a sawtooth-wave modulation in both delay lines for dif-ferent values of the modulation amplitudes: (a) ε , = 0(MDFC, constant delays), (b) ε , = 0 .
25, (c) ε , = 0 . ε , = 1. The feedback gain and the nominal delay ofthe second feedback line are fixed at K = 0 . τ = 2.The solid/dashed lines that describe the stability boundaryare calculated from Eqs. (65)–(67). The parameters of thefree-running system are as in Fig. 13. can be obtained from Eqs. (52)–(53): K (Ω) = p ± p p χ ( i Ω , ε ) − [1 − χ ( i Ω , ε )] q − χ ( i Ω , ε ) , (65) τ (Ω) = 1Ω (cid:20) arcsin (cid:18) qK (Ω) χ ( i Ω , ε ) (cid:19) + 2 nπ (cid:21) , (66) τ (Ω) = 1Ω (cid:20) − arcsin (cid:18) qK (Ω) χ ( i Ω , ε ) (cid:19) + (2 n + 1) π (cid:21) , (67)where for brevity we used the notations p = λ − K [1 − cos(Ω τ ) χ ( i Ω , ε )] , (68) q = ± ω − Ω − K sin(Ω τ ) χ ( i Ω , ε ) . (69)Finally, one might ask the question, to which extentthe extended time-delayed feedback control (ETDFC) bySocolar et al. [38, 39] is also affected by a modulationon top of the discrete exponential delay distribution cre-ated by the original method. In Fig. 13 we have al-ready presented the results for an asymmetric choice ofthe coupling gains K = 0 .
05 and K = 0 .
2. With acorresponding choice of τ = 2 τ , which gives a specialsection through the shown parameter planes, the controlscheme can already be seen as a rudimentary ETDFCscheme with τ = τ , in which the longer delays nτ with n > VI. CONCLUSIONS
Extension of the delayed feedback control of unsta-ble equilibria by introducing time-varying delays in oneor more independent delayed feedback lines in the con-trol circuit significantly enlarges the stability domain inthe control parameter space. In addition to the domainenlargement, the variable-delay feedback method is suc-cessful in stabilizing unstable states even for those caseswhere the control always fails in the constant delay case.We have shown that an analytical investigation of thecontrol domains becomes possible in the range of high-frequency delay modulations, in which case the variable-delay term in the controller can be approximated by anequivalent distributed delay term. Our approach is mo-tivated by both simulations and real experiments [45],which suggest that when the frequency of the delay mod-ulation is comparable to the system frequencies, the sys-tem dynamics can be very well approximated with thedistributed-delay effect. Restricting our attention to thishigh-frequency limit, we have investigated the control ef-ficiency of the variable-delay feedback method by consid-ering a simple normal form of an unstable focus and twodifferent realizations of the controller which are exper-imentally relevant: a non-diagonal coupling of the con-trol force realized via a phase-dependent coupling matrix,and the incorporation of an additional small constant de-lay term in all the arguments of the feedback force thatrepresents the latency of the control-line realization. Inaddition, we have explored a simple realization of a mul-tiple delay feedback controller consisting of two indepen-dent delay lines of Pyragas type with time-varying de-lays. In each case, the variable-delay feedback controlwith a finite modulation of the delay is shown superiorwith respect to the constant delay case. This renders theproposed control method promising for further practicalimplementations in real experiments.The two-dimensional linear system used in the analysisis generic for all systems with unstable fixed points of afocus type, preserving the essential features of the higher-dimensional dynamics around the equilibrium. Thus,the results obtained in this paper aim to give a succinctexplanation of the delayed feedback control mechanismwith time-dependent delay in real systems (regular orchaotic) with similar bifurcation properties. As in theoriginal time-delayed feedback control, the torsion of theorbit is a necessary condition for the control method tobe able to stabilize equilibria.Since in practice the equivalence of the distributed de-lay case holds already for fairly low modulation frequen-cies, we propose that variable-delay feedback is a conve-nient experimental method for realizing distributed delayfeedback with different delay distribution.6
Acknowledgments
This work was supported by Deutsche Forschungsge-meinschaft in the framework of SFB 910: “Control of self-organizing nonlinear systems: Theoretical methodsand concepts of application.” Thomas J¨ungling acknowl-edges support by FEDER (EU) under the project FISI-COS (FIS2007-60327). [1] E. Ott, C. Grebogi and J. A. Yorke, Phys. Rev. Lett. ,1196 (1990).[2] E. Sch¨oll and H. G. Schuster (ed.), Handbook of chaoscontrol (Wiley-VCH, Weinheim, 2008), second com-pletely revised and enlarged edition.[3] K. Pyragas, Phys. Lett. A , 421 (1992).[4] K. Pyragas, Phil. Trans. R. Soc. A , 2309 (2006).[5] K. Pyragas and A. Tamaˇseviˇcius, Phys. Lett. A , 99(1993).[6] P. Celka, Int. J. Bifurc. Chaos Appl. Sci. Eng. , 1703(1994).[7] D. J. Gauthier, D. W. Sukow, H. M. Concannon, and J.E. S. Socolar, Phys. Rev. E , 2343 (1994).[8] D. J. Christini, V. In, M. L. Spano, W. L. Ditto, and J.J. Collins, Phys. Rev. E , R3749 (1997).[9] A. S. de Paula, and M. A. Savi, Chaos, Solitons andFractals , 2981 (2009).[10] S. Bielawski, D. Derozier, and P. Glorieux, Phys. Rev. E , R971 (1994).[11] T. Dahms, V. Flunkert, F. Henneberger, P. H¨ovel, S.Schikora, E. Sch¨oll, and H.-J. W¨unsche, Eur. Phys. J.ST , 71 (2010).[12] S. Schikora, H-J. W¨unsche, and F. Henneberger, Phys.Rev. E , 026203 (2011).[13] P. Parmananda, R. Madrigal, M. Rivera, L. Nyikos, I. Z.Kiss, and V. Gaspar, Phys. Rev. E , 5266 (1999).[14] Y. Zhai, I. Z. Kiss, and J. L. Hudson, Ind. Eng. Chem.Res. , 3502 (2008).[15] E. Gravier, X. Caron, G. Bonhomme, T. Pierre, and J.L. Briancon, Eur. Phys. J. D , 451 (2000).[16] O. L¨uthje, S. Wolff, and G. Pfister, Phys. Rev. Lett. ,1745 (2001).[17] K. Hall, D. J. Christini, M. Tremblay, J. J. Collins, L.Glass, and J. Billette, Phys. Rev. Lett. , 4518 (1997).[18] C. M. Berger, J. W. Cain, J. E. S. Socolar, and D. J.Gauthier, Phys. Rev. E , 041917 (2007).[19] H. Benner, and W. Just, J. Korean Phys. Soc. , 1046(2002).[20] T. Pierre, G. Bonhomme, and A. Atipo, Phys. Rev. Lett. , 2290 (1996).[21] H. Wei, Y. Chang-Xuan, Z. Jian, X. Jin-Lin, L. Wan-Dong, F. Dong-Lai, and D. Wei-Xing, Chinese Physics , 1913 (2004).[22] J. M. Krodkiewski, and J. S. Faragher, J. Sound Vib. , 591 (2000).[23] Y. Sugimoto, K. Osuka, and T. Sugie, Journal of theRobotics Society of Japan , 435 (2005).[24] S. Steingrube, M. Timme, F. W¨org¨otter, and P. Manoon-pong, Nature Physics , 224 (2010).[25] E. Sch¨oll, Nature Physics , 161 (2010).[26] K. Yamasue, K. Kobayashi, H. Yamada, K. Matsushige,and T. Hikihara, Phys. Lett. A , 3140 (2009).[27] H. Kojima, Y. Furukawa, and P. Trivailo, Journal ofGuidance, Control, and Dynamics , 998 (2012).[28] W. Just, T. Bernard, M. Ostheimer, E. Reibold, and H. Benner, Phys. Rev. Lett. , 203 (1997).[29] P. H¨ovel and E. Sch¨oll, Phys. Rev. E , 046203 (2005).[30] S. Yanchuk, M. Wolfrum, P. H¨ovel and E. Sch¨oll, Phys.Rev. E , 026201 (2006).[31] T. Dahms, P. H¨ovel and E. Sch¨oll, Phys. Rev. E ,056201 (2007).[32] B. Fiedler, V. Flunkert, M. Georgi, P. H¨ovel and E.Sch¨oll, Phys. Rev. Lett , 114101 (2007).[33] E. W. Hooton and A. Amann, Phys. Rev. Lett. ,154101 (2012).[34] H. Nakajima, Phys. Lett. A , 207 (1997).[35] K. Pyragas, Phys. Rev. Lett. , 2265 (2001).[36] K. Pyragas, V. Pyragas, and H. Benner, Phys. Rev. E , 056222 (2004).[37] K. Pyragas and V. Pyragas, Phys. Rev. E , 036215(2006).[38] J. E. S. Socolar, D. W. Sukow and D. J. Gauthier, Phys.Rev E , 3245 (1994)[39] D. W. Sukow, M. E. Bleich, D. J. Gauthier, and J. E. S.Socolar, Chaos , 560 (1997).[40] A. Ahlborn and U. Parlitz, Phys. Rev. Lett. , 264101(2004).[41] A. Ahlborn and U. Parlitz, Phys. Rev. E , 016206(2005).[42] A. Gjurchinovski and V. Urumov, Europhys. Lett. (EPL) , 40013 (2008).[43] A. Gjurchinovski and V. Urumov, Phys. Rev. E ,016209 (2010).[44] A. Gjurchinovski, T. Sandev and V. Urumov, J. Phys.A: Math. Theor. , 445102 (2010).[45] T. J¨ungling, A. Gjurchinovski and V. Urumov, Phys.Rev. E , 046213 (2012).[46] A. Gjurchinovski and V. Urumov, Rom. J. Phys. , 36(2013).[47] Y. N. Kyrychko, K. B. Blyuss, and E. Sch¨oll, Eur. Phys.J. B , 307 (2011).[48] Y. N. Kyrychko, K. B. Blyuss, and E. Sch¨oll, Phil. Trans.R. Soc. A (2013).[49] W. Michiels, V. Van Assche and S. Niculescu, IEEETrans. Autom. Control , 493 (2005).[50] S. Schikora, P. H¨ovel, H.-J. W¨unsche, E. Sch¨oll, and F.Henneberger, Phys. Rev. Lett. , 213902 (2006).[51] T. Dahms, P. H¨ovel, and E. Sch¨oll, Phys. Rev. E ,056213 (2008).[52] S. Schikora, H. J. W¨unsche and F. Henneberger, Phys.Rev. E , 026203 (2011).[53] K. Pyragas and T. Pyragiene, Phys. Rev. E , 046217(2008).[54] A. Selivanov, J. Lehnert, T. Dahms, P. H¨ovel, A. L. Frad-kov and E. Sch¨oll, Phys. Rev. E , 016201 (2012).[55] W. Just, D. Reckwerth, E. Reibold and H. Benner, Phys.Rev. E , 2826 (1999)[56] P. H¨ovel and J. E. S. Socolar, Phys. Rev. E , 036206(2003).[57] K. Gu, S.-I. Niculescu and J. Chen, J. Math. Anal. Appl.311