Periodic Coupling inhibits Second-order Consensus on Networks
PPeriodic Coupling inhibits Second-order Consensus on Networks
Fabian Baumann, ∗ Igor M. Sokolov,
1, 2 and Melvyn Tyloo Institut f¨ur Physik, Humboldt-Universit¨at zu Berlin, Newtonstraße 15, 12489 Berlin, Germany IRIS Adlershof, Humboldt-Universit¨at zu Berlin,Zum Großen Windkanal 6, 12489, Berlin, Germany School of Engineering, University of Applied Sciences of Western Switzerland HES-SO, CH-1951 Sion, Switzerland.
Consensus algorithms on networks have received increasing attention in recent years for variousapplications ranging from animal flocking to multi-vehicle co-ordination. Building on the establishedmodel for second-order consensus of multi-agent networks, we uncover a mechanism inhibiting theformation of collective consensus states via rather small time-periodic coupling modulations. Wetreat the model in its spectral decomposition and find analytically that for certain intermediatecoupling frequencies parametric resonance is induced on a network level – at odds with the expectedemergence of consensus for very short and long coupling time-scales. Our formalism precisely pre-dicts those resonance frequencies and links them to the Laplacian spectrum of the static backbonenetwork. The excitation of the system is furthermore quantified within the theory of parametricresonance, which we extend to the domain of networks with time-periodic couplings.
I. INTRODUCTION
De-centralized collective behaviors arise in a variety ofdifferent contexts, where single units interact without acentralized authority controlling the overall process. Of-ten such interactions correspond to exchanges of informa-tion that are modeled using networks. Depending on thecontext, the information flow between two network nodesmay be interpreted differently – ranging from single neu-rons coupled via action potentials [1, 2], animal groupsexchanging cues on prey or predators [3, 4], to commu-nicating robots within vehicular platoons [5–8]. Notably,it has been shown that for many coupled dynamical sys-tems, the arising collective states strongly depend on thecharacteristics of the network topology [9].Traditionally, special attention has been paid to theemergence of consensus [10–12]. The outstanding roleof consensus states for many systems is often justified bycollective decision making, which requires agreement on apopulation level to reach a collective goal [11]. For exam-ple, fish schools have been shown to reach increased lev-els of collective sensitivity if individuals are more stronglyaligned constituting a consensus on the collective headingdirection [13]. In social systems language-based commu-nication relies on sets of quasi-static vocabularies, wherea consensus on the meaning of a word has been achieved[14]. And in vehicular platoons, consensus is often at-tributed to situations, in which all involved robotic unitsneed to be located at the same position, point in the samedirection or move at the same speed [5–7, 15, 16].Taking into account the specific physical constraintsof single units (or agents) gives rise to different typesof consensus algorithms. Considerable efforts have beendevoted to second-order consensus protocols – referringto the Newtonian dynamics of single agents – both indiscrete and continuous time [15, 17–26]. Illustrative ex- ∗ Corresponding author: [email protected] amples include applications to leader-following control[26], flocking [16] and mulit-vehicle co-orperative con-trol [18]. More specifically, second-order consensus al-gorithms were investigated for limited interaction ranges[25], time delays [19] and for the case of different networksfor position and velocity coupling [24]. Most importantlyit was shown, that networks with a (directed) spanningtree do not yield a sufficient condition for achieving con-sensus in second-order systems [27] – at odds with first-order counterparts.While the vast majority of research on consensus fo-cused on static networks [15, 17–19, 22–26], less atten-tion has been paid to the effects arising on temporal, i.e.time-varying networks. Exceptions include investigationsof multi-agent systems on switching topologies [20, 21],where the interaction network changes in discrete steps.Similar to the case of static networks – and again in con-trast to first-order systems – it was found that consensusis not necessarily achieved if infinitely many unions ofconsecutive network snapshots have a directed spanningtree [21]. Related to consensus formation is the conceptof synchronization on networks of nonlinear oscillators[28, 29], which was investigated for a further type oftemporal network. Here – in contrast to switching net-works – the coupling strengths are time-dependent, whilethe topology of the network does not change over time.In Refs. [30, 31] such network coupling was assumed tobe periodically modulated and it was found, that suchmodulation may strongly impact the synchronization be-havior of the considered systems. In Ref. [31], usingthe classical R¨ossler oscillator, network synchronizabilitycould be significantly improved by tuning the couplingfrequency. In contrast, the same periodic coupling wasfound to suppress synchronization in ensembles of glob-ally coupled first-order Kuramoto oscillators [30].In this paper we uncover yet another mechanism in-hibiting the formation of collective second-order consen-sus states on networks. It is induced, not by a switchingnetwork topology [20, 21], but the time-periodic modula-tion of coupling strengths as proposed in [30, 31]. To this a r X i v : . [ n li n . AO ] A ug FIG. 1.
Dynamics of the DO and SOC model for resonant and off-resonant coupling modulation. For ease of illustrationswe consider a small system of three agents connected via a complete graph A . The coupling strength between any two agents ismodulated time-periodically according to W ( t ) = f ( t ) A , cf. Eq. (1). In the case of a static network A for ω = 0, or sufficientlyoff-resonant modulations ( ω (cid:54) = ω ∗ ), the agents will reach a consensus in the limit of t → ∞ . In the case of the DO modelvelocities vanish and all agents come to rest in the same position. Instead, for the SOC model consensus is defined as a statein which all agents move together (being in the same place) at equal velocities. In contrast, for resonant coupling modulations( ω = ω ∗ α ) parametric excitation, on a collective network level, inhibits the establishment of global consensus states and theagents’ amplitudes diverge, as depicted in the right panel. end we investigate a general model for coupled oscillatorswith Newtonian dynamics [32] in two different settings:i) linearly coupled damped oscillators (DO) [33], and ii)the classical second-order consensus model (SOC) withvelocity alignment [17–19, 22–25]. For both models, theinteraction strengths are periodically modulated aroundwell defined mean values, which are encoded in a staticand symmetric backbone network. In the limits of veryfast and very slow modulations, we find that global con-sensus states emerge – as expected for the time-averagedstatic connectivity [18]. Interestingly, however, specificintermediate modulation frequencies, lead to an excita-tion of the system, and thereby inhibit the formation ofa global consensus.As we demonstrate, this collective dynamical featurecan be understood in terms of parametric resonance [34]on the level of network modes, revealed by the spectraldecomposition of the system. In our analysis we link thecorresponding resonance frequencies (of such parametricoscillators) to the Laplacian spectrum of the static back-bone network. We furthermore estimate the growth ratesof oscillation amplitudes in the case of resonance. As wewill briefly discuss in the Appendix using the Kuramotomodel [35] with inertia [36–38], our framework may alsobe applied to non-linearly coupled systems, close to asynchronization fixed point.In the left panel of Fig. 1 the dynamics of both consid-ered models towards consensus is illustrated. While forthe damped oscillator model (DO) a motionless consen-sus via dissipation (orange lines) emerges, the consensusstate of the SOC model preserves the total kinetic energy,such that all agents travel as one cluster with a specificfinite velocity (purple lines). In the case of resonant mod-ulation, consensus is inhibited for both models, and the agents’ amplitudes diverge, cf. right panel of Fig. 1.The rest of the paper is structured as follows. In Sec. IIwe introduce the considered general model of coupledagents. Using a spectral decomposition we derive de-coupled mode equations, which can be treated withinthe theory of parametric oscillators. This allows to linkthe resonance frequencies of both, the DO and the SOCmodel, to the Laplacian spectrum of the backbone ma-trix. In Sec. III we validate and complement our theoret-ical considerations using numerical simulations on differ-ent networks, and briefly discuss the theory’s applicabil-ity to non-linearly coupled oscillators, cf. App. A. Thework is summarized and concluded in Sec. IV. II. MODEL AND THEORY
Let us consider a network of n coupled units (oragents). The coupling matrix, W ( t ) = f ( t ) A , which en-codes the time-dependent connection strengths betweentwo agents is composed of two parts: i) a static andsymmetric backbone network, described by the adjacencymatrix A , and ii) the periodic modulation [30, 31] f ( t ) = 1 + h sin( ωt ) , (1)with frequency ω and amplitude h . If two agents i and j are connected the matrix elements A ij (= A ji ) = 1 and A ij (= A ji ) = 0, otherwise. To exclude repulsive couplingwe assume h ∈ [0 ,
1] in the following. Note, that due tothe periodicity of f ( t ), the network backbone and thecoupling matrix are connected via the temporal average A = lim t →∞ T − (cid:82) T W ( t )d t . Each agent i is associatedwith a distinct node of the network and is characterizedby a time-dependent state variable, x i ( t ), where we willomit the time dependence for brevity in the following.To model the collective dynamics we consider theestablished model for second-order consensus on net-works [17–19, 22–25], where we additionally assume time-periodic couplings. In its most general form includingdissipation [39], it is described by the following set ofsecond-order differential equations¨ x i + d ( t ) ˙ x i = − γ n (cid:88) j =1 W ij ( t ) ( x i − x j ) − µ n (cid:88) j =1 W ij ( t ) ( ˙ x i − ˙ x j ) , (2)where W ij ( t ) is the ij –th element of the symmetric cou-pling matrix W at time t and d ( t ) denotes the dampingcoefficient. The coefficients γ, µ > d > µ = 0). Thiscorresponds to an ensemble of n one-dimensional har-monic oscillators, coupled via linear springs with a time-dependent spring constant γW ij ( t ). Similar harmonicoscillators with time-invariant couplings have been inves-tigated on networks using the classical mass spring model[33]. More importantly, however, the DO model describesthe linearized dynamics of non-linearly coupled oscilla-tors, such as the second-order Kuramoto model, whichwas previously used to investigate synchronization phe-nomena on power-grids [37, 38]. In Appendix A, we willdemonstrate that the developed formalism may also beapplied to such non-linearly coupled systems. Note, thatthe agents in the damped oscillator (DO) model achieveconsensus due to the finite dissipation of d > c and consensus is characterized as x i = c , ∀ i and ˙ x i ( t → ∞ ) = 0 , ∀ i . The dynamics is illustrated inthe left panel of Fig. 1.Subsequently, we investigate the effects of periodic cou-pling for the classical model for second-order consensus(SOC) on networks [15, 17–26], where agents do notdissipate energy ( d = 0) but may establish consensusdue to a velocity alignment mechanism, for µ >
0, cf.Eq. (2). Here, consensus is usually defined as situationswith | x i − x j | = 0 , ∀ i, j and | ˙ x i − ˙ x j | = 0 , ∀ i, j – a statewith equal but potentially finite velocities [18]. In thefollowing, we will show that both models are susceptibleto time-periodic coupling, which inhibits the formation ofconsensus states. The presented formalism furthermorereveals an intimate relation between dissipation and ve-locity alignment and their roles in consensus formation.To proceed further, we express the general model, Eqs. (2), in vectorial form as¨ x + d ( t ) ˙ x = − L ( t )( γ x + µ ˙ x ) , (3)where L ( t ) denotes the time-dependent network Lapla-cian of the system. Analogous to the temporal couplingmatrix, W ( t ), L ( t ) factorizes into the time-independentLaplacian L (0) and the coupling modulation function f ( t ), i.e. L ( t ) = L (0) f ( t ). The elements of L (0) are ob-tained from the static backbone connectivity A as L (0) ij = (cid:26) − A ij , i (cid:54) = j , (cid:80) k A ik , i = j . (4)Note, that the modulation f ( t ) acts upon every connec-tion of the static backbone network A simultaneously[30, 31]. Therefore, also the time-dependent eigenvaluesof L ( t ), λ α ( t ), can be expressed as λ α ( t ) = λ (0) α f ( t ) where { λ (0) α } is the set of eigenvalues of L (0) .A spectral decomposition of Eq. (3) using x ( t ) = (cid:80) α c α ( t ) u α , involving the eigenvectors of L (0) , u α , givesrise to the following second-order differential equation¨ c α + k ( t ) ˙ c α + γλ α ( t ) c α = 0 , (5)for the α –th expansion coefficient, c α , of x ( t ) on set ofeigenvectors { u α } , where we defined k ( t ) = d ( t )+ µλ α ( t ).In this form, the effect of a time-varying coupling onthe dynamics of the system becomes particularly evi-dent: equation (5) describes the collective dynamics onthe level of network modes and manifests itself as para-metric oscillator [40]. Such oscillators are not excited byan explicit external force, but are susceptible to periodicvariations of the system’s parameters – a phenomenonwhich is called parametric resonance. Therefore, para-metric oscillators cannot be driven out of a stable equi-librium, but existing deviations from such a stable pointmay be amplified exponentially by carefully tuned peri-odic modulations, via k ( t ) or λ α ( t ). The composition of k ( t ) reveals the analogy between dissipation d > µ >
0. Both mecha-nisms may contribute to the damping of the α –th mode,and will therefore influence the dynamics in similar ways.In what follows we treat Eq. (5) analytically and iden-tify, for both model versions (DO and SOC), certain res-onant modulation frequencies ω ∗ [cf. Eq. (1)], for whichthe corresponding networked systems cannot establish aglobal consensus. At a first glance this may come as a sur-prise: diffusively coupled agents with W ( t ) >
0, gener-ally aim to assimilate, hence, one might expect the emer-gence of consensus. Equation (5), however, suggests thatperiodically modulated coupling strengths translate intocollective parametric resonance, which is, as we will seein the following, characterized by an exponential growthof the agents’ amplitudes. Here, we establish an analyt-ical link between the spectrum of the static Laplacian L (0) and the expected resonance frequencies ω ∗ α . Finally,we quantify the exponential excitation in the case of res-onant modulation.To obtain the set of resonance frequencies for bothmodels we proceed by expressing Eq. (5) in a more suit-able form using the following variable transformations, c α ( t ) = e − K ( t ) q α ( t ) , (6) K ( t ) = 12 (cid:90) t k ( t (cid:48) )d t (cid:48) . (7)Substituting those into Eq. (5) yields¨ q α + Ω α ( t ) q α = 0 , (8)the Hill equation [41], with the time-dependent frequencyΩ α ( t ), defined asΩ α ( t ) = γλ α ( t ) − k ( t ) / − ˙ k ( t ) / . (9)Note, that the time dependence in Eq. (9) either resultsfrom a periodic modulation of the damping d ( t ), the α –th time-dependent eigenvalue λ α ( t ), or both. Up tothis point, the framework includes both considered cases,which differ regarding Eq. (9). Below, we will thereforediscuss both models separately. A. Damped oscillator model (DO)
In the case of the DO model we assume for simplicitya time-independent damping coefficient d . Then Eq. (9)becomes Ω α ( t ) = γλ (0) α [1 + h sin( ωt )] − d / , (10)where we have used λ α ( t ) = λ (0) α [1 + h sin( ωt )]. Defining ω α = γλ (0) α − d / q α + ω α (cid:34) γλ (0) α hω α sin( ωt ) (cid:35) q α = 0 , (11)such that parametric excitation for the α –th mode is ex-pected for ω ∗ α ω ∗ α = 2 ω α = 2 (cid:115) γλ (0) α − (cid:18) d (cid:19) , (12)corresponding to twice the natural frequency of the sys-tem ω α – as in the case of a one-dimensional parametri-cally driven oscillator [34]. Note, that not all eigenvalues λ (0) α are expected to give rise to a resonance – only thoseeigenvalues, which fulfill the following requirements:4 γλ (0) α > d , (13) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γλ (0) α h ω α (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) > d . (14)Indeed, Eq. (13) ensures that ω ∗ α in Eq. (12) is purelyreal while Eq. (14) yields an exponential growth of the oscillations (see below). Interestingly, as both conditionsdepend on the mode index α , only a subset of the modesmay be activated by the periodic coupling depending onthe parameter values of d and h . More precisely thoseconditions set a lower bound under which natural fre-quencies cannot be activated. Apart from locating theresonance frequencies ω ∗ α in the system, it is instructiveto quantify the exponential growth of the α –th mode dueto resonant modulation.To this end, we first assume small oscillation ampli-tudes in Eq. (11), i.e. γλ (0) α h/ω α (cid:28) q α ( t ) = a α ( t ) cos[( ω α + ε ) t ]+ b α ( t ) sin[( ω α + ε ) t ] . (15)where ε defines the detuning from ω ∗ α . Note, that thisansatz will not yield an exact solution to Eq. (11), asterms with multiples of the frequency (2 ω α + (cid:15) ) have beenneglected, it is, however, sufficient in the case of smallmodulation amplitudes [34]. Inserting the expression for q α ( t ) into Eq. (11) for the particular case of ω ∗ α = 2 ω α + (cid:15) ,and neglecting terms of order O ( h ), yields − a α + a α γλ (0) α h ω α − εb α = 0 , b α + b α γλ (0) α h ω α − εa α = 0 . (16)We search the solutions to Eqs. (16) in the form of a α ( t ) , b α ( t ) ∼ e s α t , respectively, which yields the fol-lowing expression for the exponent s α s α = (cid:34) γλ (0) α h ω α (cid:35) − ε = γλ (0) α h (cid:113) γλ (0) α − d / − ε . (17)Hence, resonance for the α –th mode is expected to appearon the following interval of (cid:15) , − γλ (0) α h (cid:113) γλ (0) α − d / < ε < γλ (0) α h (cid:113) γλ (0) α − d / , (18)centered around the resonance frequency ω ∗ α . Note, thatthe width of the interval depends linearly on h , such thatthe excitation region is increased for larger modulationamplitudes. In the case of perfectly resonant modulation( ε = 0) the amplitudes a α ( t ) and b α ( t ), respectively, areexpected to grow exponentially as e s α t . Re-substitutingthis into the initial expression for the expansion coef-ficients yields c α ( t ) ∝ e ( s α − d ) t , from which conditionEq. (14) is derived. B. Second-order consensus model (SOC)
In a similar manner we may approach the model ofsecond-order consensus, which will be discussed in thefollowing. Here, we consider Eq. (3) without damping( d = 0) and with diffusive velocity coupling ( µ > k ( t ) = µλ (0) α [1 + h sin( ωt )], we getΩ α ( t ) = γλ (0) α − (cid:32) µλ (0) α (cid:33) (19)+ γλ (0) α h − h (cid:32) µλ (0) α (cid:33) sin( ωt ) − hωµλ (0) α ωt ) − (cid:32) hµλ (0) α (cid:33) sin ( ωt ) . Small modulation amplitudes, h (cid:28)
1, allow to neglectterms of the order O ( h ) in the expression for Ω α ( t ).Combining the linear cosine and sine terms, then yieldsΩ α ( t ) (cid:39) ω α (cid:32) Bω α sin( ωt + ϕ ) (cid:33) (20)with ω α = γλ (0) α − (cid:32) µλ (0) α (cid:33) , (21) B = h (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116)(cid:32) ωµλ (0) α (cid:33) + γλ (0) α − (cid:16) µλ (0) α (cid:17) (22)and ϕ = tan − ωµλ (0) α γλ (0) α − (cid:16) µλ (0) α (cid:17) . (23)As in the previous case of the DO model, the resonancefrequency corresponds to twice the natural frequency ofthe system ω α , i.e. ω ∗ α = 2 ω α = 2 (cid:118)(cid:117)(cid:117)(cid:116) γλ (0) α − (cid:32) µλ (0) α (cid:33) (24)and we get s α = (cid:20) B ω α (cid:21) − ε , (25) where c α ( t ) is expected to grow exponentially if s α t − K ( t ) >
0. The detuning range for the SOC is limitedto the interval − B/ (4 ω α ) < (cid:15) < B/ (4 ω α ) around ω ∗ α .Evaluating the integral in Eq. (7) yields K ( t ) = µλ (0) α (cid:16) t + hω (1 − cos( ωt )) (cid:17) , (26)which reduces for h (cid:28) K ( t ) (cid:39) µλ (0) α t . (27)Resonances are therefore expected if the following twoinequalities are fulfilled: Bω α − µλ (0) α > γ − µ λ (0) α > . (29)Note, that as in the case of the DO model, the velocityalignment, controlled by the parameter µ , is hamperingthe emergence of parametric resonance and therefore fos-ters consensus.In Fig. 1 the two different situations, consensus andresonance are illustrated for both models. In each casethe basic model setup consists of the same static back-bone network, an all-to-all network of three nodes. Thepairwise coupling strengths, W ij ( t ), are either modulatedoff-resonantly ( ω (cid:28) ω ∗ ) (left panel) or in a perfectly res-onant manner as ω = ω ∗ (right panel). While the sys-tem approaches a global consensus in the former case,consensus is inhibited by parametric excitation for reso-nant modulation. Instead, the nodes’ amplitudes growexponentially and do not reach a steady state where | x i − x j | = 0 ∀ i, j .The presented results are rather general. They allow todetermine the resonance frequencies ω ∗ α for systems on ar-bitrary static backbone structures A , where the couplingbetween agents is modulated by a periodic function f ( t )of small amplitudes. In the following we investigate thoseresonance effects for both models on different backbonenetworks. While the eigenvalues λ (0) α of L (0) , required forthe determination of { ω ∗ α } , can be obtained analyticallyfor simple regular graphs, one needs to determine themnumerically for more complex or random networks. III. NUMERICAL RESULTS
In order to validate the theoretical analysis we per-form numerical simulations. As suggested by the theory,the agents in both models will or will not reach a globalconsensus, depending on the modulation frequency ω . Todetect the corresponding ranges of resonance numericallywe shall proceed as follows.We integrate Eqs. (3) using a 4th–order Runge-Kuttascheme [42], by which we obtain the time evolution of FIG. 2.
Amplitude growth for resonant modulation with ω = ω ∗ and d = µ = 0. The exponential growth of theagents’ amplitudes is quantitatively captured by the time-evolution of the α –th expansion coefficient c α ( t ) ∼ e st (blackline). The green trajectory corresponds to the time evolutionof a single agent i , and its local maximum values are depictedas green dots. The inset shows the exponential growth as asubset of local maximum values of the same agent (green dots)and c α ( t ) (black line) over many orders of magnitude. Theshown results are for a complete graph with n = 100 nodes. the system, x ( t ), for a specific modulation frequency ω .Averaging the integral (cid:82) x ( t ) d t , obtained for each mod-ulation frequency ω , over multiple simulation runs withrandom initial conditions (IC), yields A ( ω ) = (cid:42) T (cid:90) x ( t ) d t (cid:43) { IC } , (30)which quantifies the excitation of the system in termsof the nodes’ squared amplitudes and (implicitely) as afunction of ω . In such a spectrum A ( ω ), resonances ap-pear as pronounced peaks, while regions of consensus areflat. Note that, the value of A (0) corresponds to the caseof no coupling modulation, i.e. the reference model of astatic network A . Without loss of generality we initial-ize both models as resting non-consensus states, i.e. weset the state variables at time t = 0 randomly on theinterval x i (0) ∈ I init and assume vanishing initial veloci-ties ˙ x i (0) = 0 , ∀ i . Unless indicated differently we take I init = [ − , h = 0 . d = 0 . γ = 1) and( γ = 1, µ = 0 .
01) for the DO and SOC model, respec-tively. For reasons of comparability between the mod-els, we chose equal values for the damping and velocityalignment coefficients. As a consequence of this choice,the predicted resonance frequencies are similar for bothmodels and it is ensured that the conditions (13)-(14)and (28)-(29) are satisfied. In the following we discussour theoretical results on complete graphs, ring networksand finally on random Erd˝os-R´enyi networks [43].
FIG. 3.
Resonance on a complete graph.
The vertical or-ange and dashed magenta lines locate the predicted resonancefrequency ω ∗ for each model on a complete graph of n = 5nodes. The orange and magenta spectra correspond to theDO and SOC model, respectively, on the interval ω ∈ [0 , ε around ω ∗ , cf.Eq. (18) and Eq. (25). The inset shows, for the DO model,the theoretically predicted values of ω ∗ (black line), and asextracted from simulations (orange crosses) as functions ofthe system size. Amplitude growth
In a first step we check the theoretical predictions withrespect to the exponential growth of the system’s ampli-tudes in the case of resonant modulation. Here, we con-sider for simplicity the case without both damping ( d =0) and velocity alignment ( µ = 0), where the two consid-ered models coincide: the expressions for the resonancefrequencies, cf. Eqs. (12),(24), yield ω ∗ α = 2 (cid:113) γλ (0) α , andthe exponent s α becomes s α = h (cid:113) γλ (0) α . (31)The case of such parametric excitation is shown in Fig. 2for a complete graph of n = 100 nodes, where the singleresonance is obtained as ω ∗ = 20. For clarity the tra-jectory and the local maximum values of a single node( i = 1) are shown in green. The trajectories x j ( t ) ofthe remaining agents are colored in grey. Clearly, forresonant modulation the nodes’ amplitudes increase ex-ponentially over time. The growth is well described by c α ( t ) ∼ e st (black line), where s is given in Eq. (31).As the inset of Fig. 2 demonstrates, the behavior is notonly captured for short times, but the approximation alsoholds over many orders of magnitude for larger times t (cid:29) FIG. 4.
Resonances on ring networks.
The orange and dashed magenta vertical lines locate the analytically predictedresonance frequencies for both models on a ring network of n = 5 (a) and n = 25 (b) nodes, respectively. The resonance spectraare computed as averages over 10 random initial conditions and a simulation time of T = 100. The inset in panel (a) depictsthe ratio A ( ω ∗ ) / A ( ω ∗ ) for the DO model, as a function of the simulation time T (orange rhombuses). The time-dependentratio is well described by e s − s ) t (black line), where s and s are computed using ω and ω , cf. Eq. (17). Complete graphs
In order to verify the theoretical predictions for theresonance frequencies, we first consider complete graphs,or all-to-all connected networks – a general topology, of-ten considered in studies of coupled oscillators [38]. Theeigenvalues of L (0) for a complete graph are given as λ = 0 and λ = n , where λ has the algebraic multi-plicity n − ω ∗ = 2 (cid:115) γn − (cid:18) d (cid:19) , (32)for the DO model and ω ∗ = 2 (cid:114) γn − (cid:16) µn (cid:17) . (33)in the case of the SOC model, respectively.In Fig. 3 we depict the corresponding resonance spectraon a complete graph of n = 5 nodes. Here, and in all sub-sequent figures the orange and magenta lines depict thespectra of the DO and SOC model, respectively. Verticallines locate the predicted values for ω ∗ α , colored accord-ing to the associated model. For both models, the singlepeaks of integrated square amplitudes A ( ω ) are predictedremarkably well by the theoretical values of ω ∗ . Further-more, the finite widths of the peaks are approximatelycaptured by the theoretical bounds of detuning ε .The thermodynamic limit of n → ∞ yields differentbehaviors for both models. For fixed values of µ and γ , no resonance is expected for the SOC model, as theresonance conditions are violated for large n . Instead, inthe case of the DO model, the damping coefficient doesnot influence the expression for ω ∗ significantly. Hence, the resonance frequency is well described by ω ∗∞ (cid:39) √ n ,which is supported by numerical simulations for networksizes up to n = 1000, see inset of Fig. 3(a). For growingsystem sizes the establishment of a global consensus inthe DO model is therefore increasingly robust against lowfrequency modulations in f ( t ). Solely high frequency-components may parametrically excite the system andinhibit consensus. Rings
In contrast to the previous case of complete graphs,the Laplacian spectrum of a ring network may containmore than two distinct eigenvalues, which may generallyleads to multiple resonances. The set { λ (0) α } distributesover the interval λ (0) α ∈ [0 ,
4] and the eigenvalues are an-alytically given as, λ (0) α = 2 − k α ) , with k α = 2 π ( α − n . (34)This leads to the following expressions for the resonancefrequencies: ω ∗ α = 2 (cid:115) − k α ) − (cid:18) d (cid:19) (35)in the case of the DO model, and like-wise ω ∗ α = 2 (cid:112) − µ + (2 µ −
2) cos( k α ) − µ cos ( k α ) (36)for the SOC model. Note, that in the descriptive case ofa ring network of n = 5 nodes, both distinct eigenvalues λ (0) α > ω ∗ α ,i.e. all resonance conditions are satisfied for d = µ =0 . FIG. 5.
Resonances on Erd˝os-R´enyi networks.
The orange and dashed magenta vertical lines locate the analyticallypredicted resonance frequencies ω ∗ α for both models on Erd˝os-R´enyi networks of n = 25 (a) and n = 100 (b) nodes, respectively.The spectra are averaged over 10 random initial conditions, and the shaded blue area shows the corresponding standarddeviation. In panel (b) the modulation amplitude is h = 0 .
4. The inset in panel (b) shows the spectrum obtained for a largermodulation amplitude h = 0 .
6, where half resonant frequencies excite the system also for the SOC model.
The corresponding spectra are depicted in Fig. 4(a),where both resonances show up as pronounced peaks,matched by their theoretical predictions. For the chosensimulation time of T = 100, the heights of the resonancepeaks, A ( ω ∗ ) and A ( ω ∗ ), are separated by several ordersof magnitude, where the separation is larger in the case ofthe DO model. In the inset of Fig. 4(a) we demonstrate,for the DO model, that the ratio of observed peak heights,may also be captured by the theory. It is explained bydifferent exponential growth rates s α , and can, for large t , be approximated as (cid:82) c d t/ (cid:82) c d t ∼ ( s /s ) e s − s ) t (black line), obtained from Eq. (17). It captures well thenumerical results over many orders of magnitude.In Fig. 4(b) we show the resonance spectra for largerring networks of n = 25 nodes, where an important con-sequence of the bounded Laplacian spectrum becomesapparent. At odds with the behavior of the DO modelon complete graphs, increasing the number of nodes forring networks does not increase the highest resonancefrequency beyond all limits. Instead, for n → ∞ , thebounded intervals of resonance frequencies are graduallyfilled with resonance peaks, as suggested by Eqs. (34).Once two resonant frequencies are sufficiently close toeach other, i.e. within their respective ε –ranges, alsointermediate frequencies may parametrically excite thesystem, such that each resonance spectrum becomes aquasi-continuous one. The effect becomes pronouncedfor higher modulation frequencies ω →
4, where the in-creased values of A ( ω ), between two predicted resonancefrequencies ω ∗ α , indicate a partially off-resonant excita-tion of the system. Random Networks
Finally we demonstrate that the theory holds on ran-dom networks, i.e. that resonance frequencies (based on L (0) ) can be predicted reliably on networks with disor-der. In the previous cases of regular ring networks andcomplete graphs, the Laplacian spectrum was given an-alytically. This is not possible for random networks andthe corresponding eigenvalues need to be determined nu-merically.In panels (a) and (b) of Fig. 5 we depict the resonancespectra for two Erd˝os-R´enyi (ER) networks of n = 10and n = 100 nodes, respectively. Especially, in the for-mer case, the random network structure is revealed bythe rather irregular positions of the resonance peaks –an observation, which is directly related to the randomspectrum of the corresponding Laplacian L (0) . Never-theless, the resonances for both models, are successfullypredicted by Eqs. (12) and (21).In the case of the larger ER network, cf. panel (b) ofFig. 5, the resonance frequencies are very close to eachother, such that the single vertical lines are barely dis-tinguishable. The overlapping detuning ranges (cid:15) , hence,lead to a quasi-continuous excitation of the system onthe whole interval of roughly ω ∈ [8 ,
13] for both models.Interestingly, the resonance spectrum of the DO model,shows a further feature, known for general parametricoscillators: the excitation of the system at half the reso-nance frequency [34]. As we show in the inset of Fig. 5for h = 0 .
6, this effect becomes increasingly pronouncedfor larger modulation amplitudes h , and is not limited tothe DO model, but also appears for the SOC model. IV. SUMMARY AND CONCLUSIONS
In this work we studied parametric resonances emerg-ing in a general model of networked agents, which we dis-cussed in two different settings: i) a model for dampedharmonic oscillators, and ii) a well established model forsecond-order consensus formation, without dissipationand with velocity alignment. We considered those mod-els for time-periodic coupling modulations, and foundthat, for certain intermediate time scales of such mod-ulation, the collective dynamics is not captured by thetime-averaged coupling A = (cid:104) W ( t ) (cid:105) t , which would leadto a global consensus state. Instead, periodic couplingmodulations may induce parametric resonance on a col-lective level, where single network modes are excited andthe emergence of consensus is inhibited. The developedformalism correctly predicts those resonance frequencies,within the theory of parametric resonance, which we ex-tended to the domain of networks with time-varying cou-pling. As we demonstrated, resonant modulation givesrise to exponentially growing mode amplitudes, whichcan be quantified in terms of the associated eigenvaluesof the static Laplacian L (0) . In this context, both dissi-pation and velocity alignment counteract the excitationof the system and therefore foster the establishment ofglobal consensus states.In comparison both considered models, DO and SOC,are similarly susceptible to time-periodic coupling mod-ulations. The excitation for a fixed resonance frequency,however, is increased for the DO model, as suggested bythe numerically obtained resonance spectra. For smallsystem sizes and rather weak dissipation, as well as veloc-ity alignment, the expected resonance frequencies are inclose vicinity. Using the example of complete graphs, weillustrated that this generally changes for large networks.Increasingly large eigenvalues would not give rise to res-onances in the SOC model, which emphasized the en-hanced positive effect on consensus formation for velocityalignment, compared to dissipation, cf. condition (29).Interestingly, we found that the formalism also holdson networks of non-linearly coupled second-order Ku-ramoto oscillators close to a synchronization fixed point.It correctly predicts modulation frequencies, for whichthe synchronized fixed point becomes unstable. This didnot yield exponentially growing but finite oscillation am-plitudes.Our work might have important implications for therobustness of second-order consensus algorithms appliedin real systems. For example, in order to disrupt con-sensus states within networked vehicular platoons or co-operating autonomous robots, signals which periodicallyweaken the units’ sensors, may be used. Note that para-metric resonance might also be used to infer characteris-tic networks features, in cases where the probing of thesystem is not localized on certain nodes, but performedon the overall coupling. ACKNOWLEDGMENTS
This work was developed within the scope of theIRTG 1740/TRP 2015/50122-0 and funded by theDFG/FAPESP and the Swiss National Science Founda-tion under grant No. 200020 182050. We thank AdrianPacheco, J¨org N¨otel and Thomas Peron for fruitful dis- cussions and helpful comments on the manuscript.
Appendix A: Non-linear oscillators
In the main text we considered (linear) diffusively cou-pled second-order consensus systems. However, Eq. (3)may also be regarded as the linearization of non-linearlycoupled oscillators.Traditionally, one of the most studied models of non-linearly coupled oscillators is the Kuramoto model [35].Second-order Kuramoto models were previously used toinvestigate synchronization phenomena, such as partic-ular cases of consensus formation [45] or frequency syn-chronization in power grids [46]. By exchanging the diffu-sive coupling by a sinusoidal coupling, setting µ = 0 andswitching to the common notation for phase oscillators,i.e. x i → φ i , Eq. (3) yields the second-order Kuramotomodel as¨ φ i + d ˙ φ i = ω i − N (cid:88) j =1 W ij ( t ) sin( φ i − φ j ) . (A1)At odds with the linear model of Eq. (3) the boundednessof the sinusoidal coupling inhibits an arbitrarily large ex-citation of the system in the case of resonant modulation.Instead, the nodes’ trajectories x i are characterized bylarge but bounded oscillation amplitudes.Without loss of generality, we take (cid:80) i ω i = 0 whichyields a synchronous state with ˙ φ i = 0 ∀ i . Furthermorewe assume that Eq. (A1) has a stable fixed point φ i, fora specific backbone network A , which yields ω i = N (cid:88) j =1 W ij ( t ) sin( φ i, − φ j, ) . (A2)For small modulation amplitudes in f ( t ), i.e. h (cid:28) δφ i ( t ) = φ i ( t ) − φ i, into Eq. (A1) yields δ ¨ φ i + d δ ˙ φ i = − N (cid:88) j =1 W ij ( t ) cos( φ i, − φ j, )( δφ i, − δφ j, ) . (A3)In vectorial form we have δ ¨ φ + dδ ˙ φ = − L ( t ; φ ) δ φ , (A4)where L ( t ; φ ) denotes a Laplacian matrix, weighted bythe angle differences at the fixed point φ . Therefore, therelations derived in Sec. II still hold – where { λ (0) α } nowdenotes the eigenvalue spectrum of the weighted Lapla-cian L ( t ; φ ).Note, that we focus on deviations from a global con-sensus for which holds φ ,i = φ ,j , ∀ i (cid:54) = j . In this casethe weighted Laplacian reduces to the one used in thelinear model [Eq. (3)], i.e. we have L ( t, φ ) = L ( t ).0 FIG. 6.
Resonance for the second-order Kuramotomodel.
The red line locates the predicted resonance fre-quency ω ∗ for in the case of a complete graph of n = 5 nodesclose to the consensus fixed point. The blue line depicts theresonance spectrum for the second-order Kuramoto model onthe interval ω ∈ [0 , I init = [ − . , . φ i ( t ) for resonant modulationwith ω ∗ (cid:39) . The inset of Fig. 6 shows the time evolution of theagents’ phases φ i for resonant modulation ( ω = ω ∗ ) ona complete graph of n = 5 nodes. As expected, inthe case of resonance, consensus is inhibited. In con-trast to the linear models (DO and SOC), however, os-cillation amplitudes do not grow exponentially withoutbounds. Instead, phase oscillations reach finite ampli-tudes, which would vanish for off-resonant modulation,where agents will reach a global consensus, as in thelinear case, depicted in Fig. 1. This behavior is con-sistent for the whole range of relevant modulation fre-quencies as shown in Fig. 6, where a pronounced peakappears for the same frequency as predicted for theDO model, cf. Fig. 3. Note, that the linearizationof the sinusoidal coupling only holds for the dynamicsclose to the fixed point of consensus. To observe reso-nance in the numerical simulations, the Kuramoto modelneeds therefore to be probed close to this fixed point( | φ i (0) − φ j (0) | (cid:28) | ˙ φ i (0) − ˙ φ j (0) | (cid:28) ∀ i (cid:54) = j ), forwhich we decreased the width of the initial interval of φ i (0) to I init = [ − . , . [1] D. S. Bassett and O. Sporns, Nature neuroscience ,353 (2017).[2] E. M. Izhikevich, Dynamical systems in neuroscience (MIT press, 2007).[3] S. B. Rosenthal, C. R. Twomey, A. T. Hartnett, H. S.Wu, and I. D. Couzin, Proceedings of the NationalAcademy of Sciences , 4690 (2015).[4] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, andI. D. Couzin, Proceedings of the National Academy ofSciences , 18720 (2011).[5] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, in , Vol. 2 (IEEE, 2003) pp.2010–2015.[6] V. Gupta, B. Hassibi, and R. M. Murray, in , Vol. 1 (IEEE, 2003) pp. 504–509.[7] A. Rodriguez-Angeles and H. Nijmeijer, in , Vol. 2 (IEEE, 2003) pp. 1514–1519.[8] J. Werfel, K. Petersen, and R. Nagpal, Science , 754(2014).[9] M. A. Porter and J. P. Gleeson, Frontiers in Applied Dy-namical Systems: Reviews and Tutorials (2016).[10] M. H. DeGroot, Journal of the American Statistical As-sociation , 118 (1974).[11] A. Baronchelli, Royal Society open science , 172189(2018).[12] C. Castellano, S. Fortunato, and V. Loreto, Reviews ofmodern physics , 591 (2009).[13] I. Couzin, Nature , 715 (2007).[14] M. Hechter and K.-D. Opp, Social norms (Russell SageFoundation, 2001).[15] W. Ren, R. W. Beard, and E. M. Atkins, IEEE Control systems magazine , 71 (2007).[16] H.-T. Zhang, Z. Cheng, G. Chen, and C. Li, IEEE Trans-actions on Circuits and Systems I: Regular Papers ,1599 (2015).[17] W. Ren and E. Atkins, in AIAA Guidance, Navigation,and Control Conference and Exhibit (2005) p. 6238.[18] W. Ren and E. Atkins, International Journal of Robustand Nonlinear Control: IFAC-Affiliated Journal , 1002(2007).[19] W. Yang, A. L. Bertozzi, and X. Wang, in (IEEE, 2008)pp. 2926–2931.[20] J. Qin, C. Yu, and S. Hirche, IEEE Transactions onIndustrial Informatics , 986 (2012).[21] W. Ren, in (IEEE,2007) pp. 1431–1436.[22] W. Yu, G. Chen, M. Cao, and J. Kurths, IEEE Trans-actions on Systems, Man, and Cybernetics, Part B (Cy-bernetics) , 881 (2009).[23] W. Yu, G. Chen, and M. Cao, Automatica , 1089(2010).[24] D. Goldin and J. Raisch, Asian Journal of Control ,30 (2014).[25] X. Ai, S. Song, and K. You, Automatica , 329 (2016).[26] J. Hu, G. Chen, and H.-X. Li, in Proceedings of the 30thChinese control conference (IEEE, 2011) pp. 4819–4824.[27] W. Ren and R. W. Beard, IEEE Transactions on auto-matic control , 655 (2005).[28] A. Pikovsky, J. Kurths, M. Rosenblum, and J. Kurths, Synchronization: a universal concept in nonlinear sci-ences , Vol. 12 (Cambridge university press, 2003).[29] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, andC. Zhou, Physics reports , 93 (2008).[30] S. Li, X. Wang, and S. Guan, New Journal of Physics , 113013 (2018). [31] S. Li, N. Sun, L. Chen, and X. Wang, Phys. Rev. E ,012304 (2018).[32] H. G. Oral, E. Mallada, and D. F. Gayme, in (IEEE, 2017) pp. 1688–1694.[33] M. Zhan, S. Liu, and Z. He, PloS one , e82161 (2013).[34] L. Landau and E. Lifschitz, “Lehrbuch der theoretischenphysik. band 1: Mechanik,” (1984).[35] Y. Kuramoto, in International symposium on mathemat-ical problems in theoretical physics (Springer, 1975) pp.420–422.[36] H.-A. Tanaka, A. J. Lichtenberg, and S. Oishi, PhysicaD: Nonlinear Phenomena , 279 (1997).[37] F. D¨orfler, M. Chertkov, and F. Bullo, Proceedings ofthe National Academy of Sciences , 2005 (2013).[38] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, Physics Reports , 1 (2016).[39] T. W. Grunberg and D. F. Gayme, IEEE Transactionson Control of Network Systems , 456 (2018).[40] L. Turyn, Quarterly of applied mathematics , 389(1993).[41] G. W. Hill, Acta mathematica , 1 (1886).[42] W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T.Vetterling, et al. , “Numerical recipes,” (1989).[43] P. Erd˝os and A. R´enyi, Publ. Math. Inst. Hung. Acad.Sci , 17 (1960).[44] F. R. Chung and F. C. Graham, Spectral graph theory ,92 (American Mathematical Soc., 1997).[45] S. Patterson and B. Bamieh, Proc. of the 49 th IEEE CDC(2010).[46] A. R. Bergen and D. J. Hill, IEEE Transaction on PowerAppartus and Systems