Time-Delayed Interactions in Networks of Self-Adapting Hopf Oscillators

Abstract

A network of coupled limit cycle oscillators with delayed interactions is considered. The parameters characterizing the oscillator’s frequency and limit cycle are allowed to self-adapt. Adaptation is due to time-delayed state variables that mutually interact via a network. The self-adaptive mechanisms ultimately drive all coupled oscillators to a consensual cyclostationary state, where the values of the parameters are identical for all local systems. They are analytically expressible. The interplay between the spectral properties of the coupling matrix and the time delays determines the conditions for which convergence towards a consensual state takes place. Once reached, this consensual state subsists even if interactions are removed. In our class of models, the consensual values of the parameters depend neither on the delays nor on the network’s topologies.

The farther back you can look, the farther forward you are likely to see
Winston Churchill

1. Introduction

The harmonic excitation of an elementary-damped harmonic oscillator 𝑥̈(𝑡)+𝖺𝑥̇(𝑡)+𝖿𝑥(𝑡)_system=𝗋sin(𝗐𝑡)_environment (1) with 𝖺, 𝖿, 𝗋 and 𝗐>0 produces the well-known asymptotic response (c.f. [1]): 𝑥(𝑡)=𝗄cos(𝗐𝑡−𝜗), (2) where 𝗄 and 𝜗 depend on the control parameters 𝖺, 𝗋, and 𝗐. By construction, the oscillating environment here materialized by the input 𝗋sin(𝗐𝑡) is totally insensitive to the oscillator 𝑥(𝑡), implying that 𝑥(𝑡) in (1) is slaved by the external forcing. The next stage of complexity is to replace the harmonic oscillator in (1) by a Lienard system: 𝑥̈(𝑡)+𝖱(𝑥(𝑡),𝑥̇(𝑡))+𝖿𝑥(𝑡)=𝗋sin(𝗐𝑡), (3) where 𝖱(𝑥(𝑡),𝑥̇(𝑡)) is a nonlinear controller. In absence of external excitation in (3) (i.e., when 𝗋=0), we assume 𝖱 to asymptotically drive the orbits towards a stable limit cycle which is independent of the initial conditions—the paradigmatic illustration being here the Van der Pol oscillator. When 𝗋≠0 and for a suitably selected range of parameters, the time asymptotic response of (3) can be qualitatively written as (c.f. [2, 3]) 𝑥(𝑡)=𝖲(𝑡) (4) with 𝖲(𝑡) being 𝖺 a synchronized signal with the same periodicity as the environment (i.e., 𝖲(𝑡+(2𝜋/𝗐))=𝖲(𝑡)). By construction, the external forcing 𝗋sin(𝗐𝑡) in (3) is, as before, insensitive to the Lienard oscillator. In the resulting synchronized regime, the oscillator 𝑥(𝑡) is caught by the external excitation—in other words, the system adjusts itself to the environment but the environment remains insensible to the system. Observe that the dynamical response given by (4) only subsists as long as 𝗋sin(𝗐𝑡) acts on the system. That is, as soon as the environment effect is removed (i.e., 𝗋=0 in (3)), the system (i.e., the limit cycle oscillator), after a transient time, recovers its original behavior—converges towards its limit cycle.

In our present paper, we will extend the previous classical system-environment relationship in order to allow more realistic situations where mutual interactions permanently affect both the system and the environment. Such adaptive mechanisms can modify individual dynamics on a permanent basis. To stylize this new situation, the basic dynamics given by (3) is modified as 𝑥̈(𝑡)(𝑥(𝑡),𝑥̇(𝑡(𝑡)𝑥(𝑡)(𝜔(𝑡)𝑡),𝑓̇(𝑡)+𝖱))+𝑓=𝗋sin=𝖠^𝖿(𝑥(𝑡),𝑥̇(𝑡)),𝜔̇(𝑡)=𝖠^𝗐(𝑥(𝑡),𝑥̇(𝑡)), (5) where 𝖿↦𝑓(𝑡) and 𝗐↦𝜔(𝑡) are no longer constant parameters but variables of the global dynamics (and hence time dependent) and the functions 𝖠^𝖿 and 𝖠^𝗐 capture the mutual adaptation of the system-environment dipole. Hence, the system of (5) has now to be considered globally—the system and its environment, are allowed to adaptively coevolve.

The particular case of (5) when 𝖠^𝖿≡0 (i.e., 𝑓(𝑡)∶=𝖿 for all 𝑡) has been studied in [4, 5]. This type of dynamical system provides a cornerstone of bioinspired robotics where legs (or arms) of robots may be modeled by damped oscillators as (1) (i.e., 𝖱(𝑥(𝑡),𝑥̇(𝑡))=𝖺𝑥̇(𝑡)). To ensure a maximum leg stride, the damped oscillators must be excited by 𝗋sin(𝜔(𝑡)𝑡) at the damped oscillator's resonant frequency (i.e., 𝜔(𝑡)⋍𝑓 for all 𝑡). However, due to structural changes on the robot (adding load, lengthening legs), the resonant frequency of the damped oscillator has to be adjusted: 𝖿↦𝖿. Hence, 𝖠^𝗐 drives 𝜔(𝑡) towards the value 𝖿 to systematically guarantee the excitation at resonant frequency and thus maximum leg stride. Frequency adaptation (as well as amplitude adaption) is also important in movement assistance (e.g., retrain the nervous system, assist people with movement disorder), where robots and human beings must work in synchronous. An example of an exoskeleton for the human elbow was studied in [6].

The case where neither 𝖠^𝖿 nor 𝖠^𝗐 is trivial has been covered in [7, 8], where the authors not only considered two adaptive coupled limit cycle oscillators, but 𝑛 mutually interacting through a network. Here, self-adapting oscillators can be applied to robot formation modeling. Each individual robot belonging to a swarm, circulating around a specific point, adapts its angular velocity in order to lower the amount of exchanged information to maintain the formation. Self-adaptation in networks is also considered in [9], where here the control signals (and not the local systems) of the variables of the CPG adapt, and this, to quickly react to new situations and produce several different behavioral patterns.

Building on what has been done in previous contributions [4, 5, 7, 8], we here consider a network of limit cycle oscillators interaction with time-delayed state variables. The general form of our dynamical system in the phase-radius coordinates (i.e., polar coordinates 𝜙_𝑘 and 𝑟_𝑘) is 𝜙̇_𝑘(𝑡)=𝖯𝜙_𝑘(𝑡),𝑟_𝑘(𝑡);Ω_𝑘−𝖼_𝑘𝜕𝖵𝜕𝜙_𝑘𝜙()()𝑟̇𝑡−𝗍,𝑟𝑡−𝗍,𝑘=1,…,𝑛,_𝑘(𝑡)=𝖱𝜙_𝑘(𝑡),𝑟_𝑘(𝑡);Ω_{𝑘localdynamics}−𝖼_𝑘𝜕𝖵𝜕𝑟_𝑘𝜙(𝑡−𝗍),𝑟(𝑡−𝗍_{couplingdynamics}),𝑘=1,…,𝑛, (6) where 𝖯 and 𝖱 govern the local dynamics, 𝜙=(𝜙₁,…,𝜙_𝑛) and 𝑟=(𝑟₁,…,𝑟_𝑛) are the state variables, Ω_𝑘 is a parameter set determining the local characteristics and 𝖼_𝑘>0 are coupling strengths. The coupling dynamics is the gradient of a potential𝖵 depending on state variables with time delay 𝗍. Adaptation of the local dynamics is accomplished by letting the constant parameters Ω_𝑘 become time-dependent (i.e., Ω_𝑘↦Ω_𝑘(𝑡)) with their own dynamics given by adaptive mechanisms 𝖠_𝑘. Delays being ubiquitous in applications, their influence on adaptive processes is worth to be investigated, and so the general formalization is 𝜙̇_𝑘(𝑡)=𝖯𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),Ω_𝑘(𝑡)−𝖼_𝑘𝜕𝖵𝜕𝜙_𝑘𝜙()()𝑟̇𝑡−𝗍,𝑟𝑡−𝗍,𝑘=1,…,𝑛_𝑘(𝑡)=𝖱𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),Ω_𝑘(𝑡)_{localdynamics}−𝖼_𝑘𝜕𝖵𝜕𝑟_𝑘𝜙(𝑡−𝗍),𝑟(𝑡−𝗍_{couplingdynamics})Ω̇,𝑘=1,…,𝑛_𝑘(𝑡)=𝖠_𝑘𝜙(𝑡−𝗍),𝑟(𝑡−𝗍_{adaptivemechanism}),𝑘=1,…,𝑛. (7) We therefore confer to Ω_𝑘 the status of variables of the whole dynamical system. Let us remark that in this present contribution, adaptation occurs in the local systems. Note that in [10], the authors introduce, with the help of the speed-gradient method, an adaptive mechanism on the coupling constant that multiplies the delayed interactions.

This paper is organized as follows: in Section 2 we define the three components that together form the global system. We then discuss the dynamics of our model in Section 3. An application is presented in Section 4 which is then followed by some numerical experiments in Section 5. Finally, we conclude in Section 6.

2. Networks of Hopf Oscillators with Adaptive Mechanisms

We now present explicitly the local and coupling dynamics as well as the adaptive mechanisms on which we will focus.

2.1. Local Dynamics

Each node of the network is equipped with a local dynamical system. In this contribution, a local system is a Hopf oscillator presented here in its polar coordinates: 𝖯𝜙_𝑘,𝑟_𝑘;Ω_𝑘=𝗐_𝑘𝖱𝜙,𝑘=1,…,𝑛,_𝑘,𝑟_𝑘;Ω_𝑘=−𝑟²_𝑘−𝗋_𝑘𝑟_𝑘,𝑘=1,…,𝑛. (8) The state variables are (𝜙_𝑘,𝑟_𝑘) and Ω_𝑘={𝜔_𝑘,𝗋_𝑘} are, for the time being, fixed and constant parameters. The parameter 𝗐_𝑘 controls the frequency of the 𝑘th oscillator given by the phase dynamics 𝖯. The radial dynamics 𝖱 produces a stable circular limit cycle with radius √𝗋_𝑘.

2.2. Coupling Dynamics

Associated to an 𝑛 vertex connected and undirected network, denote by 𝐴 the weighted adjacency matrix with positive entries 𝑎_𝑘,𝑗⩾0. Let 𝐿 be the corresponding Laplacian matrix (𝐿∶=𝐷−𝐴, where 𝐷 is the diagonal matrix with 𝑑_𝑘,𝑘∶=∑^𝑛_𝑗=1𝑎_𝑘,𝑗). The coupling dynamics is given by the gradient of the positive semidefinite function 𝖵1𝜙,𝑟∶=2𝑟∣𝐿_cos𝑟=12^𝑛_𝑘=1𝑟_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑗cos𝜙_𝑘−𝜙_𝑗⩾0 (9) with 𝜙=(𝜙₁,…,𝜙_𝑛) and 𝑟=(𝑟₁,…,𝑟_𝑛) and where the matrix 𝐿_cos has entries 𝑙_𝑘,𝑗cos(𝜙_𝑘−𝜙_𝑗). This matrix is positive semidefinite since all its eigenvalues are positive (i.e. nonnegative). This is a direct application of Gershgorin’s circle theorem [11]: for any eigenvalue 𝜁_cos of 𝐿_cos, there exists 𝑘 such that 𝜁_cos−𝑙_𝑘,𝑘⩽^𝑛_𝑗≠𝑘𝑙_𝑘,𝑗cos𝜙_𝑘−𝜙_𝑗 (10) and so |𝜁_cos−𝑙_𝑘,𝑘|⩽∑^𝑛_𝑗≠𝑘|𝑙_𝑘,𝑗cos(𝜙_𝑘−𝜙_𝑗)|⩽∑^𝑛_𝑗≠𝑘|𝑙_𝑘,𝑗|=𝑙_𝑘,𝑘. In particular, the coupling dynamics is 𝖼_𝑘𝜕𝖵𝜕𝜙_𝑘𝜙,𝑟=−𝖼_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑘𝑟_𝑗sin𝜙_𝑘−𝜙_𝑗𝖼,𝑘=1,…,𝑛,_𝑘𝜕𝖵𝜕𝑟_𝑘𝜙,𝑟=𝖼_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑗cos𝜙_𝑘−𝜙_𝑗,𝑘=1,…,𝑛, (11) where 𝖼_𝑘>0 are coupling strengths.

2.3. Adaptive Mechanisms

In this section, we now allow the fixed and constant parameters 𝗐_𝑘 and 𝗋_𝑘 to (I)become time dependent, that is, Ω_𝑘={𝗐_𝑘,𝗋_𝑘}⇝(𝜔_𝑘(𝑡),𝜌_𝑘(𝑡))=Ω_𝑘(𝑡), for 𝑘=1,…,𝑛, (II)and each of them have their own dynamics, depending solely on the state variables 𝜙 and 𝑟, that is, 𝖠_𝑘𝜙,𝑟=𝖠^𝗐_𝑘𝜙,𝑟,𝖠^𝗋_𝑘𝜙,𝑟for𝑘=1,…,𝑛. (12)

Among the numerous variants for changing the values of the parameter, we focus on those presented in [7, 8, 12], that is, 𝖠^𝗐_𝑘𝜙,𝑟=𝗌_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑘𝑟_𝑗sin𝜙_𝑘−𝜙_𝑗𝖠,𝑘=1,…,𝑛,_𝗋𝑘𝜙,𝑟=−𝗌_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟²_𝑗,𝑘=1,…,𝑛, (13) where 𝑙_𝑘,𝑗 are the entries of 𝐿 and 𝗌_𝑘 are “susceptibility constants”: the larger 𝗌_𝑘 is, the stronger the influence on 𝜔_𝑘 and 𝜌_𝑘 is. Conversely, oscillators with small 𝗌_𝑘 are reluctant to modify their frequency and their limit cycle radius.

3. Network's Dynamical System with Delay

We now discuss the resulting dynamics in the presence of a time delay 𝗍⩾0 affecting both the coupling dynamics and the adaptive mechanisms. We hence consider 𝜙̇_𝑘(𝑡)=𝜔_𝑘(𝑡)+𝖼_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑘(𝑡−𝗍)𝑟_𝑗()𝜙𝑡−𝗍×sin_𝑘(𝑡−𝗍)−𝜙_𝑗(),𝑟̇𝑡−𝗍_𝑘(𝑡)=−𝑟_𝑘(𝑡)²−𝜌_𝑘(𝑡)𝑟_𝑘(𝑡)−𝖼_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑗(𝑡−𝗍)cos𝜙_𝑘(𝑡−𝗍)−𝜙_𝑗(),𝜔̇𝑡−𝗍_𝑘(𝑡)=𝗌_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑘(𝑡−𝗍)𝑟_𝑗(𝑡−𝗍)sin𝜙_𝑘(𝑡−𝗍)−𝜙_𝑗(),𝜌̇𝑡−𝗍_𝑘(𝑡)=−𝗌_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑗(𝑡−𝗍)², (14) for 𝑘=1,…,𝑛. For (14), we have the following.

Two Constants of Motion. The functions 𝖩(𝜔)=^𝑛_𝑗=1𝜔_𝑗𝗌_𝑗,𝖪𝜌=^𝑛_𝑗=1𝜌_𝑗𝗌_𝑗 (15)

are constants of motion: in other words, if 𝜔(𝑡)=(𝜔₁(𝑡),…,𝜔_𝑛(𝑡)) and 𝜌(𝑡)=(𝜌₁(𝑡),…,𝜌_𝑛(𝑡)) are orbits of (14), then d[𝖩(𝜔(𝑡))]=⟨(𝜔(𝑡∣𝜔̇(𝑡=d𝖪𝜌(𝑡d𝑡∇𝖩)))⟩)=𝜌(𝑡)∣𝜌̇(𝑡d𝑡∇𝖪)=0. (16)

Existence of a Consensual Oscillatory State. We can explicitly exhibit a consensual oscillatory state. Indeed, for given 𝜔_𝑐 and 𝜌_𝑐>0, 𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),𝜔_𝑘(𝑡),𝜌_𝑘(𝑡)∶=𝜔_𝑐𝑡,√𝜌_𝑐,𝜔_𝑐,𝜌_𝑐 (17)

for all 𝑡∈[−𝗍,0]∪ℝ_⩾0 is a consensual orbit of (14).

Observe that in absence of the radial component 𝑟_𝑘 and without adaptation, (14) yields the famous Kuramoto model with delays [13, 14]. In these contributions, the authors considered the coupling dynamics with delay of the form: ∑^𝑛_𝑗=1𝑙_𝑘,𝑗sin(𝜙_𝑘(𝑡)−𝜙_𝑗(𝑡−𝗍)), that is, delays concern the “exterior” variables 𝜙_𝑗 and not the “local” variable 𝜙_𝑘. This can be done here for the coupling dynamics but not for the adaptive mechanism on 𝜔_𝑘, if we require 𝖩 to be a constant of motion 𝖩.

In the absence of time delay (i.e., 𝗍=0) and under appropriate conditions (c.f. [12] for details), the adaptive mechanisms tune the value of frequencies 𝜔_𝑘 and the radii 𝜌_𝑘 of the attractors so that the global dynamical system is driven into a consensual oscillatory state. In other words, we have the following limit: lim_𝑡→∞𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),𝜔_𝑘(𝑡),𝜌_𝑘(𝑡)−𝜔_𝑐𝑡,√𝜌_𝑐,𝜔_𝑐,𝜌_𝑐=0∀𝑘, (18) with constants 𝜔_𝑐 and 𝜌_𝑐 and where ‖⋅‖ is the Euclidean norm. The consensual state is permanent (i.e., even if interactions are switched off, all local dynamics still oscillate with the same frequency and same amplitude). Let us now discuss the conditions for which Limit (18) holds when the global dynamics is affected by a time delay (i.e., t>0).

Convergence towards a Consensual Oscillatory State. The Limit (18) raises two issues: (1) the existence itself and (2) the limit values 𝜔_𝑐 and 𝜌_𝑐. For expository reasons, we first discuss the limit values and then the convergence conditions.

Limit Values. Thanks to the constant of motions in (15), we have 𝖩(𝜔(0))=𝖩(𝜔(𝑡))∀𝑡,𝖪𝜌(0)=𝖪𝜌(𝑡)∀𝑡, (19) with 𝜔(𝑡)=(𝜔₁(𝑡),…,𝜔_𝑛(𝑡)) and 𝜌(𝑡)=(𝜌₁(𝑡),…,𝜌_𝑛(𝑡)) orbits of (14) with given initial conditions. Supposing that Limit (18) holds, we hence have 𝖩(𝜔(0))=lim_𝑡→∞𝖩(𝜔(𝑡))=𝜔_𝑐^𝑛_𝑗=11𝗌_𝑗,𝖪𝜌(0)=lim_𝑡→∞𝖪𝜌(𝑡)=𝜌_𝑐^𝑛_𝑗=11𝗌_𝑗, (20) and so the asymptotic values are analytically expressed as 𝜔_𝑐=∑^𝑛_𝑗=1𝜔_𝑗(0)/𝗌_𝑗∑^𝑛_𝑗=11/𝗌_𝑗,𝜌_𝑐=∑^𝑛_𝑗=1𝜌_𝑗(0)/𝗌_𝑗∑^𝑛_𝑗=11/𝗌_𝑗. (21) It is important to emphasize that the consensual values 𝜔_𝑐 and 𝜌_𝑐 do not depend on the network topology (i.e., not on 𝐿), nor on the initial conditions of the state variables (i.e., not on (𝜙_𝑘(0),𝑟_𝑘(0)) nor on the time delay (i.e., not on 𝗍).

Convergence Conditions. To this aim we study the first-order approximation of (14) in the vicinity of Solution (17) and assume that linear stability analysis is sufficient to infer convergence conditions for the nonlinear system. Accordingly, we study the asymptotic behavior of the small perturbations 𝜖_𝜙𝑘(𝑡),𝜖_𝑟𝑘(𝑡),𝜖_𝜔𝑘(𝑡), and 𝜖_𝜌𝑘(𝑡) and write 𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),𝜔_𝑘(𝑡),𝜌_𝑘(𝑡)=𝜔_𝑐𝑡+𝜖_𝜙𝑘(𝑡),√𝜌_𝑐+𝜖_𝑟𝑘(𝑡),𝜔_𝑐+𝜖_𝜔𝑘(𝑡),𝜌_𝑐+𝜖_𝜌𝑘(𝑡). (22) Taking into account the constant of motions, we impose that ^𝑛_𝑗=1𝜖_𝜔𝑗(0)𝗌_𝑗=0,^𝑛_𝑗=1𝜖_𝜌𝑗(0)𝗌_𝑗=0. (23)

First Order Approximation. Rearranging the variables (i.e., the first 𝑛 are the 𝜙_𝑘, the second 𝑛 are the 𝑟_𝑘, the third 𝑛 are the 𝜔_𝑘, and finally the last 𝑛 are the 𝜌_𝑘), the first-order approximation of (14) is 𝜖̇_𝜙𝜖̇_𝑟𝜖̇_𝜔𝜖̇_𝜌=𝟎𝟎𝐼𝑑𝟎𝟎−2𝜌_𝑐𝐼𝑑𝟎√𝜌_𝑐𝜖𝐼𝑑𝟎𝟎𝟎𝟎𝟎𝟎𝟎𝟎^𝜙𝜖_𝑟𝜖_𝜔𝜖_𝜌+−𝜌_𝑐[𝖼][𝖼]𝐿𝟎𝟎𝟎𝟎−𝐿𝟎𝟎−𝜌_𝑐[𝗌]√𝐿𝟎𝟎𝟎𝟎−2𝜌_𝑐[𝗌]𝜖̌𝐿𝟎𝟎^𝜙𝜖̌_𝑟𝜖̌_𝜔𝜖̌_𝜌, (24) with the 𝑛×𝑛 identity matrix 𝐼𝑑, diagonal matrices [𝖼] and [𝗌] with, respectively, the coupling strengths and the susceptibility constants as entries and 𝜖_𝜙∶=(𝜖_𝜙1,…,𝜖_𝜙𝑛), 𝜖_𝑟∶=(𝜖_𝑟1,…,𝜖_𝑟𝑛), 𝜖_𝜔∶=(𝜖_𝜔1,…,𝜖_𝜔𝑛) and 𝜖_𝜌∶=(𝜖_𝜌1,…,𝜖_𝜌𝑛), and where the delayed perturbations are 𝜖̌_𝜙=𝜖_𝜙(𝑡−𝗍), 𝜖̌_𝑟=𝜖_𝑟(𝑡−𝗍), 𝜖̌_𝜔=𝜖_𝜔(𝑡−𝗍)𝜖̌_𝜌=𝜖_𝜌(𝑡−𝗍).

Diagonalization. Suppose now that [𝗌]=𝗊[𝖼] for some positive constant 𝗊 and let 𝑂 denote an orthogonal matrix (i.e., 𝑂^⊤𝑂=𝑂𝑂^⊤=𝐼𝑑) with real entries such that 𝑂^⊤[𝖼]^1/2𝐿[𝖼]^1/2𝑂=[𝜁], with [𝜁] being a diagonal matrix with the eigenvalues of the symmetric matrix [𝖼]^1/2𝐿[𝖼]^1/2 on its diagonal. The signs of these coincide with those of the eigenvalues of 𝐿: they are all strictly positive except for one, that is, zero. Hence, without loss of generality, one takes 𝜁₁=0 and 𝜁_𝑘>0 for 𝑘=2,…,𝑛.

Changing the basis of System (24) with a 4×4 bloc matrix (each bloc of size 𝑛×𝑛) with 𝑂^⊤[𝖼]^−1/2 on its diagonal, we can decompose the original system into 2𝑛  2-dimensional systems of the form𝜀̇_𝜙𝑘𝜀̇_𝜔𝑘=01𝜀00𝜀^𝜙𝑘_𝜔𝑘+−𝜌_𝑐𝜁_𝑘0−𝜌_𝑐𝗊𝜁_𝑘0𝜀̌_𝜙𝑘𝜀̌_𝜔𝑘𝜀̇_𝑟𝑘,𝜀̇_𝜌𝑘=−2𝜌_𝑐√𝜌_𝑐𝜀00𝜀^𝑟𝑘_𝜌𝑘+−𝜁_𝑘0−2√𝜌_𝑐𝗊𝜁_𝑘0𝜀̌_𝑟𝑘𝜀̌_𝜌𝑘, (25a)(25b)

with 𝜀_𝜙∶=𝑂^⊤[𝖼]^−1/2𝜖_𝜙 (resp., for 𝜖_𝑟,𝜖_𝜔,𝜖_𝜌) and 𝜀_𝜙̌∶=𝑂^𝖳[𝖼]^−1/2𝜖̌_𝜙 (resp. for 𝜖̌_𝑟,𝜖̌_𝜔,𝜖̌_𝜌) for delayed perturbations obtained after the change of basis.

The case 𝑘=1 is worked out in the appendix. For 𝑘≠1, let us focus on the 2-dimensional systems and we rewrite (25a)-(25b) as linear second-order time delayed differential equations 𝜀̈_𝜙𝑘(𝑡)+𝜌_𝑐𝜁_𝑘𝜀̇_𝜙𝑘(𝑡−𝗍)+𝜌_𝑐𝗊𝜁_𝑘𝜀_𝜙𝑘()𝜀̈𝑡−𝗍=0,_𝑟𝑘(𝑡)+2𝜌_𝑐𝜀̇_𝑟𝑘(𝑡)+𝜁_𝑘𝜀̇_𝑟𝑘(𝑡−𝗍)+2𝜌_𝑐𝗊𝜁_𝑘𝜀_𝑟𝑘(𝑡−𝗍)=0. (26) The convergence towards a consensual state is hence determined by the asymptotic stability of the zero solution of (26). Stability follows if, and only if, all roots of the corresponding characteristic equations have strictly negative real parts (c.f. [15] for details). For (25a), one can apply Theorem 3.3 in [16] which states that in this case the zero solution is asymptotically stable if and only if(I)𝗍𝜌_𝑐𝜁_𝑘<𝜋2, (27a)(II)𝗍²𝜌_𝑐𝗊𝜁_𝑘𝑧²(𝑧)𝜋<cos,where𝑧istheuniquesolutionin0,2ofsin(𝑧)=𝗍𝜌_𝑐𝜁_𝑘𝑧, (27b)for 𝑘=2,…,𝑛. We emphasize that the consensual value 𝜌_𝑐 influences the condition for convergence whereas it is not the case for 𝜔_𝑐. This leads to the idea that shaping the attractor is more delicate than tuning the angular velocity. This has been observed in [17].

Adaptation only on 𝜔_𝑘. If there is no adaptation on the radii (i.e., 𝜌_𝑘(𝑡)∶=𝗋 for all 𝑘 and 𝑡), then (25b) reduced to 𝜀̇_𝑟𝑘(𝑡)=−2𝗋𝜀_𝑟𝑘(𝑡)−𝜁_𝑘𝜀_𝑟𝑘(𝑡−𝗍) (28) for which the zero solution is asymptotically stable provided (c.f. [18] for details) 4𝗋²<𝜁²_𝑘,𝗍<𝗍̌with𝗍̌=cot⁻¹−2𝗋/𝜁²_𝑘−4𝗋²𝜁²_𝑘−4𝗋². (29) The zero solution is unstable when 𝗍>𝗍̌. Note that stability is guaranteed for any 𝗍 when 𝜁²_𝑘<4𝗋² or 𝜁²_𝑘=4𝗋²≠0.

No Time Delay in the Coupling Dynamics. If there is no time delay in the coupling dynamics (i.e., the coupling dynamics is defined as in (11) with no delay), then (26) becomes, for 𝑘=2,…,𝑛, 𝜀̈_𝜙𝑘(𝑡)+𝜌_𝑐𝜁_𝑘𝜀̇_𝜙𝑘(𝑡)+𝜌_𝑐𝗊𝜁_𝑘𝜀_𝜙𝑘()𝜀̈𝑡−𝗍=0,_𝑟𝑘(𝑡)+2𝜌_𝑐+𝜁_𝑘𝜀̇_𝑟𝑘(𝑡)+2𝜌_𝑐𝗊𝜁_𝑘𝜀_𝑟𝑘(𝑡−𝗍)=0. (30a)(30b)Invoking Theorem 3.5 in [16], the zero solution for (30a) and (30b) is asymptotically stable if and only if.

For (30a)𝗍²𝜌_𝑐𝗊𝜁_𝑘+𝑧²+𝗍²𝜌_𝑐𝜁_𝑘²(𝑧)𝜋cos>0,where𝑧istheuniquesolutionin0,2of𝑧sin(𝑧)=t𝜌_𝑐𝜁_𝑘cos(𝑧), (31)

for (30b) −2𝗍²𝜌_𝑐𝗊𝜁_𝑘+𝑧²+𝗍²2𝜌_𝑐+𝜁_𝑘²(𝑧)𝜋cos>0,where𝑧istheuniquesolutionin0,2of𝑧sin(𝑧)=𝗍2𝜌_𝑐+𝜁_𝑘cos(𝑧), (32) for 𝑘=2,…,𝑛.

Summary. For a network (with arbitrary topology but with symmetric, positive entries adjacency matrix) of Hopf oscillators (as defined in Section 2.1 by (8)) interacting through time delayed Kuramoto type coupling (as defined in Section 2.2 by (11)) and with time-delay adaptive mechanisms (as defined in Section 2.3 by (13)) on the frequencies and amplitudes of the local systems—in other words, for (14), we have the following. (I)two constants of motions (c.f. (15)), (II)the existence of a consensual oscillatory state (c.f. (17)),(III)the consensual oscillatory state is linearly (i.e., locally) stable if all the roots of characteristic equations corresponding to (26) have strictly negative real parts, (IV)if there is no adaptation on the radii, the consensual oscillatory state is linearly (i.e., locally) stable if (27a), (27b), and (29) hold, (V)if there is no delay in the coupling dynamics, the consensual oscillatory state is linearly (i.e., locally) stable if (31), and (32) hold.

3.1. Miscellaneous Remark: Delayed Stabilization Mechanism

In this section we discuss the particular case arising when the delay is introduced only in the stabilization mechanism (i.e., dissipative part) of the local dynamics. The dynamical system is 𝜙̇_𝑘(𝑡)=𝜔_𝑘(𝑡)+𝖼_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑘(𝑡)𝑟_𝑗(𝑡)sin𝜙_𝑘(𝑡)−𝜙_𝑗(𝑡),𝑟̇𝑘=1,…,𝑛,_𝑘(𝑡)=−(𝑟_𝑘(𝑡−𝗍)²−𝜌_𝑘(𝑡)𝑟_𝑘(𝑡)−𝖼_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑗(𝑡)cos𝜙_𝑘(𝑡)−𝜙_𝑗(𝑡)𝜔̇,𝑘=1,…,𝑛,_𝑘(𝑡)=𝗌_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑘(𝑡)𝑟_𝑗(𝑡)sin𝜙_𝑘(𝑡)−𝜙_𝑗(𝑡)𝜌̇,𝑘=1,…,𝑛,_𝑘(𝑡)=−𝗌_𝑘^𝑛_𝑗=1𝑙_𝑘,𝑗𝑟_𝑗(𝑡)²,𝑘=1,…,𝑛. (33)

Equations (33) still admit the existence of a consensual oscillatory state and two constants of motion as in (15) and in (17), respectively. Linear stability analysis of the consensual state reduces to the study of 𝜀̈_𝜙𝑘(𝑡)+𝜌_𝑐𝜁_𝑘𝜀̇_𝜙𝑘(𝑡)+𝜌_𝑐𝗊𝜁_𝑘𝜀_𝜙𝑘(𝑡)𝜀̈=0,_𝑟𝑘(𝑡)+2𝜌_𝑐𝜀̇_𝑟𝑘(𝑡−𝗍)+2𝜌_𝑐𝗊𝜁_𝑘𝜀_𝑟𝑘(𝑡)=0, (34a)(34b)for 𝑘=2,…,𝑛. The zero solution for (34a) is asymptotically stable. For (34b), we apply Theorem 3.4 in [16] which states that in this case the zero solution is asymptotically stable if and only if(I)∃𝑧∈ℤ_≥0{}()𝜋∶=0,1,2,…suchthat2𝑧−1𝜋+2<𝗍2𝜌_𝑐𝗊𝜁_𝑘𝜋<2𝑧𝜋+2 (35)(II)0>−2𝗍𝜌_𝑐)𝜋>max(2𝑧−1𝜋+2−2𝗍²𝜌_𝑐𝗊𝜁_𝑘()(),−𝜋2𝑧−1𝜋+𝜋/22𝑧𝜋+2+2t²𝜌_𝑐𝗊𝜁_𝑘((,2𝑧𝜋+𝜋/2)) (36) for 𝑘=2,…,𝑛.

4. Applications

Conceptually, the problem of reliably distributing time and frequency among several spatial remote locations is a “leitmotiv” in applications ranging from basic metrology, navigation and position determination, signal processing, computer communications, energy distribution networks, swarms robotics, bioengineering, multiagents systems, life sciences, acoustics, and musical art to give but only a highly nonexhaustive list. Presently, a strong research impetus is devoted to complex interacting oscillating systems able to exhibit self-adaptive capabilities leading to a resilient consensual dynamic. Whatever the configurations under study, communication delays between the collection of interacting subparts of the global systems are physically unavoidable. Depending on the underlying time scales, delays do strongly affect the resulting dynamics. Our class of models explicitly study the influence of delays and in particular their destabilizing effects, that modify the instantaneous behavior. By an appropriate tuning of control parameters (e.g., susceptibility constants), our class of models offer, via a unique formalism, the possibility to continuously explore interacting configurations ranging from slave-master (i.e., system-environment relationship) to fully decentralized regimes.

Alternatively, we may view this problematic in the context of soft-controlled systems which presently receive a sustain attention [19]. Here, a swarm of agents is infiltrated by a lure agent (sometimes called a shill in economy). While the lure exhibits all the features of any ordinary agent, it can be externally controlled by an operator. As the interactions between the lure and any agent of the swarm remain unaffected (i.e., the lure remains incognito to ordinary agents), the external control of the lure can ultimately drive the whole population to a specific configuration. In our class of dynamics, a suitable choice of the susceptibility constant of a given local system (i.e., oscillator) may convert it into a shill. Indeed, in view of (21), the ultimate consensual values 𝜔_𝑐 and 𝜌_𝑐 are weighted averages. Such weighted averages can be made to be strongly dependent on a very insensitive shill-stubborn to any external influence (i.e., with a very low susceptibility constant).

In absence of time delays, convergence towards a consensual state is observed even for large heterogeneities (widely dispersed initial frequencies and radii and large discrepancies among the susceptibility constants). Hence, a shill agent can be easily introduced. However, our present study shows, that time delays restrict the conditions for convergence towards a consensual state. As a consequence, the implementation of a lure is more delicate matter (i.e., here, the dynamics is far more sensitive to the value taken by the susceptibility constants).

5. Numerical Simulations

Adaptation on 𝜔_𝑘 and 𝜌_𝑘. We consider four Hopf oscillators interacting on a network with topology as in Figure 1. We choose the coupling strengths and susceptibility constants as {𝖼₁,…,𝖼₄}={0.1,2,5,3} and 𝗌_𝑘=0.8c_𝑘 for 𝑘=1,…,4 (i.e., 𝗊=0.8). The time delay is 𝗍=0.1. The initial function (i.e., history) is (𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),𝜔_𝑘(𝑡),𝜌_𝑘(𝑡))=(𝑡,1,1,1) for 𝑡∈[−𝗍,0[, having a jump at 𝑡=0 with values (𝜙_𝑘(0), 𝑟_𝑘(0), 𝜔_𝑘(0), 𝜌_𝑘(0)) that are randomly uniformly drawn from ]−0.1,0.1[×]0.9,1.1[×]0.9,201.1[×]0.9,1.1[ with the exception for 𝜔₁(0)=0.9. The 𝜌_𝑘(0) are rescaled such that the consensual value 𝜌_𝑐 is one.

Figure 1 

Network topology.

The resulting dynamics is shown in Figure 2. With the same initial conditions, we carry out another numerical simulation with here 𝗍=0.12. This violates Condition (27b) for 𝑘=4 and hence the network does not converge towards a consensual state. This is shown in Figure 3.

(a)

(b)

Figure 2 

Time evolution of 𝜔𝑘 (Figure 2(a)) and 𝜌𝑘 (Figure 2(b)) for four Hopf oscillators, interacting through the network in Figure 1. In both figures, all variables converge towards a constant and common consensual value.

(a)

(b)

Figure 3 

Time evolution of 𝜔𝑘 (Figure 3(a)) and 𝜌𝑘 (Figure 3(b)) for four Hopf oscillators, interacting through the network in Figure 1. In both figures, no common consensual value is reached.

Note that in Figure 2(a), the 𝜔_𝑘(0) converge close to 0.9, that is, close to the initial value 𝜔₁(0). This because the first oscillator’s susceptibility is “small” and hence it is this local system that acts as a shill. It interacts with its neighbor in the same way as they act with it. It can control the behavior of the network and this without being connected to all other local systems.

In the extreme case, when 𝗌_𝑘=0 for all 𝑘 (i.e., no adaptation), (14) describe the dynamics of coupled oscillators with different limit cycles and frequencies and delayed interactions. For small frequency and attractor's shape heterogeneities, the network is able to synchronize. In this case (i.e., without adaptive mechanisms for 𝜔_𝑘 and 𝜌_𝑘), several numerical simulations show that for delayed time 𝗍=0.1175 synchronization is not systematically attained (i.e., it depends on the initial conditions)—whereas it is attained in the absence of the time delay (i.e., 𝗍=0). On the other hand, when the adaptive mechanisms are switched on (i.e., 𝗌_𝑘=0.8c_𝑘>0 for all 𝑘), a consensual state is reached—the linear stability criteria are still satisfied with 𝗍=0.1175. Summarizing, here adaptation enhances the syntonization capability of the network.

Adaptation only on 𝜔_𝑘. Two Hopf oscillators, both having the same radius for the attractor (i.e., 𝜌_𝑘=𝗋=0.1 for 𝑘=1,2), are coupled with coupling strengths and susceptibility constants as {𝖼₁,𝖼₂}={1,14} and 𝗌_𝑘=𝖼_𝑘 for 𝑘=1,2 (i.e., 𝗊=1). The time delay is 𝗍=0.1057. The initial function (i.e., history) is (𝜙_𝑘(𝑡),𝑟_𝑘(𝑡),𝜔_𝑘(𝑡))=(𝑡,√0.1,1) for 𝑡∈[−𝗍,0[, having a jump at 𝑡=0 with values (𝜙_𝑘(0),𝑟_𝑘(0),𝜔_𝑘(0)) that are randomly uniformly drawn from ]−0.1,0.1[×]√0.1−0.05, √0.1+0.05[×]0.9,1.1[. Under this configuration, Conditions (27a)-(27b) are satisfied (hence the zero solution of (25a) is asymptotically stable) but not Condition (29)—the time delay is too large. Figure 4 displays the resulting dynamics. Observe that the radii do not converge towards their attractor. However, the oscillators still adapt their frequencies.

(a)

(b)

Figure 4 

Time evolution of 𝑟𝑘 (Figure 4(a)) and 𝜔𝑘 (Figure 4(b)) for two Hopf oscillators. The two oscillators manage to adapt their frequencies while their radii do not converge towards their respective attractor.

6. Conclusion

It is more the rule than the exception that parameter adaptation in dynamical systems can be achieved via delayed mechanisms. This therefore converts ordinary differential equations arising in absence of delay to functional differential equations. Stability issues become more difficult to discuss since, dealing with functional differential equations, an infinite number of degrees of freedom is introduced into the dynamics. For the class of oscillatory networks with parametric adaptation we here considered, we are able to observe how the time delay affects the adaptation mechanisms. While it is intuitively expected, that large delays are likely to destabilize the dynamics, we are here able to analytically quantify the underlying critical delays. The analytical linear stability discussion is made possible since, for our class of dynamics, stability issues can be reduced to the study of two linear second order functional differential equations for which suitable theorems can be found. Finally, we emphasize that numerical simulations show that adaptation may enhance the emergence of common dynamical pattern whereas in classical synchronization, time delays may be too large for synchronous motion to be attained.

Appendix

Assuming that all perturbations 𝜀_𝜙𝑘,𝜀_𝑟𝑘,𝜀_𝜔𝑘,𝜀_𝜌𝑘,𝜀̌_𝜙𝑘,𝜀̌_𝑟𝑘,𝜀̌_𝜔𝑘 and 𝜀̌_𝜌𝑘 for 𝑘=2,…,𝑛 converge to zero, let us now study the case for 𝑘=1. Here, 𝜁₁=0 and so 𝜀̇_𝜙1=𝜀_𝜔1,𝜀̇_𝜔1𝜀̇=0,_𝑟1=−2𝜌_𝑐𝜀_𝑟1+√𝜌_𝑐𝜀_𝜌1,𝜀̇_𝜌1=0, (A.1) and therefore 𝜀_𝜔1(𝑡)=𝜀_𝜔1(0) and 𝜀_𝜌1(𝑡)=𝜀_𝜌1(0) for all 𝑡. Both of these constants 𝜀_𝜔1(0) and 𝜀_𝜌1(0) are zero. This is because the first orthonormal base vector (i.e., the normalized eigenvector for the eigenvalue 𝜁₁=0) is ℭ(1/√𝖼₁,…,1/√𝖼_𝑛) (with ℭ∶=(∑^𝑛_𝑗=1(1/𝖼_𝑗))^−1/2) and the first coordinate of 𝑂^⊤[𝖼]^−1/2𝜖_𝜔 and 𝑂^⊤[𝖼]^−1/2𝜖_𝜌 is, respectively, 𝜀_𝜔1(0)=ℭ^𝑛_𝑗=1𝜖_𝜔𝑗(0)𝖼_𝑗=ℭ𝗊^𝑛_𝑗=1𝜖_𝜔𝑗(0)𝗌_𝑗,𝜀_𝜌1(0)=ℭ^𝑛_𝑗=1𝜖_𝜌𝑗(0)𝖼_𝑗=ℭ𝗊^𝑛_𝑗=1𝜖_𝜌𝑗(0)𝗌_𝑗, (A.2) since we supposed that 𝗌_𝑗=𝗊𝖼_𝑗 for all 𝑗. These two sums are zero according to (23). Therefore, 𝜀̇_𝑟1=−2𝜌_𝑐𝜀_𝑟1 (i.e., lim_𝑡→∞𝜀_𝑟1(𝑡)=0) and 𝜀_𝜙1(𝑡)=𝜀_𝜙1(0) for all 𝑡. This allows one to conclude that all perturbations 𝜀_𝑟𝑘,𝜀_𝜌𝑘, and 𝜀_𝜔𝑘 decay for all 𝑘. We now need to study how the perturbations on the phases evolve. Since 𝜖_𝜙=[𝖼]^1/2𝑂𝜖_𝜙, then lim_𝑡→∞𝜖_𝜙𝑘(𝑡)=lim_𝑡→∞√𝖼_𝑘^𝑛_𝑗=1𝑜_𝑘,𝑗𝜀_𝜙𝑗(𝑡)=√𝖼_𝑘𝑜_𝑘,1lim_𝑡→∞𝜀_𝜙1(𝑡)=ℭ𝜀_𝜙1(0)=ℭ^2𝑛_𝑗=1𝜖_𝜙𝑗(0)𝖼_𝑗=∑^𝑛_𝑗=1𝜖_𝜙𝑗(0)/𝖼_𝑗∑^𝑛_𝑗=11/𝖼_𝑗, (A.3) since lim_𝑡→∞𝜀_𝜙𝑘(𝑡)=0 for 𝑘=2,…,𝑛, 𝑜_𝑘,1=ℭ1/√𝖼_𝑘 (𝑘=1,…,𝑛) and the first coordinate of the product 𝑂^⊤[𝖼]^−1/2𝜖_𝜙 is 𝜀_𝜙1(0)=ℭ∑^𝑛_𝑗=1(𝜖_𝜙𝑗(0)/𝖼_𝑗). Hence, all perturbations converge towards zero except those on the phase that all converge towards a constant (i.e., average phase perturbation). This corresponds to a phase shift.

Acknowledgments

M.-O. Hongler acknowledges partial support from the ESF (European Science Foundation) under the project entitled Exploring the Physics of Small Devices. J. Rodriguez acknowledges the support from the DFG-IRTG 1132 (Deutsche Forschungsgemeinschaft—International Research Training Group) under the project entitled Internationales Graduiertenkolleg—Stochastics and Real World Models.