A stochastic approach to the synchronization of coupled oscillators
AA STOCHASTIC APPROACH TO THE SYNCHRONIZATION OF COUPLEDOSCILLATORS
UMBERTO BICCARI AND ENRIQUE ZUAZUA
Abstract.
This paper deals with an optimal control problem associated with the Kuramoto model describingthe dynamical behavior of a network of coupled oscillators. Our aim is to design a suitable control functionallowing us to steer the system to a synchronized configuration in which all the oscillators are aligned on thesame phase. This control is computed via the minimization of a given cost functional associated with thedynamics considered. For this minimization, we propose a novel approach based on the combination of astandard Gradient Descent (GD) methodology with the recently-developed Random Batch Method (RBM)for the efficient numerical approximation of collective dynamics. Our simulations show that the employmentof RBM improves the performances of the GD algorithm, reducing the computational complexity of theminimization process and allowing for a more efficient control calculation. Introduction
Synchronization is a common phenomenon which has been observed in biological, chemical, physical andsocial systems for centuries and has attracted the interest of researcher in a diversified spectrum of scientificfields.Common examples of synchronization phenomena often cited in review articles include groups of syn-chronously chirping crickets ([44]), fireflies flashing in unison ([5]), superconducting Josephson junction ([45]),or crowds of people walking together that will tend to synchronize their footsteps ([36]).Roughly speaking, synchronization means that a network of several periodic processes with differentnatural frequencies reaches an equilibrium configuration sharing the same common frequency as a result oftheir mutual interaction.This concept is closely related to the one of consensus for multi-agent systems, widely analyzed in manydifferent frameworks including collective behavior of flocks and swarms, opinion formation, and distributedcomputing (see [2, 3, 23, 27, 28]). In broad terms, consensus means to reach an agreement regarding a certainquantity of interest that depends on the state of all agents.Synchronization is also a key issue in power electrical engineering, for instance in the model and stabilityanalysis of utility power grids ([6, 13, 30, 35]). Indeed, large networks of connected power plants need to besynchronized to the same frequency in order to work properly and prevent the occurrence of blackouts.Synchronization phenomena are most often characterized by the so-called Kuramoto model ([21]), describingthe dynamical behavior of a (large) network of oscillators in a all-to-all coupled configuration in which everyoscillator is connected with all the others. This model extends the original studies by Winfree in the contextof mutual synchronization in multi-oscillator systems based on a phase description ([46]). In particular, inKuramoto’s work, synchronization appears as an asymptotic pattern which is spontaneously reached by thesystem when the interactions among the oscillators are sufficiently strong.In some more recent contribution, control theoretic methods have been employed to analyze the synchro-nization phenomenon. For instance, in [8] the authors design passivity-based controls for the synchronization
Key words and phrases. coupled oscillators, Kuramoto model, optimal control, synchronization, Gradient Descent, RandomBatch Method.This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020research and innovation programme (grant agreement NO. 694126-DyCon). The work of both authors was partially supported bythe Grant MTM2017-92996-C2-1-R COSNET of MINECO (Spain) and by the Air Force Office of Scientific Research (AFOSR)under Award NO. FA9550-18-1-0242. The work of E.Z. was partially funded by the Alexander von Humboldt-Professorshipprogram, the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grantagreement No.765579-ConFlex, the Grant ICON-ANR-16-ACHN-0014 of the French ANR and the Transregio 154 Project‘‘Mathematical Modelling, Simulation and Optimization Using the Example of Gas Networks’’ of the German DFG. a r X i v : . [ n li n . AO ] M a y f multi-agent systems, with application to the general problem of multi-robot coordination. In [33], feedbackcontrol laws for the stabilization of coupled oscillators are designed and analyzed. In [29, 42], the authorspropose methods for the suppression of synchrony in a globally coupled oscillator network, based on (possiblytime-delayed) feedback schemes. Finally, [24] deals with the problem of desynchronizing a network ofsynchronized and globally coupled neurons using an input to a single neuron. This is done in the spirit ofdynamic programming, by minimizing a certain cost function over the whole state space.In this work, we address the synchronization problem for coupled oscillators through the constructionof a suitable control function via an appropriate optimization process. To this end, we propose a novelapproach which combines a standard Gradient Descent (GD) methodology with the recently-developedRandom Batch Method (RBM, see [20]) for the efficient numerical approximation of collective dynamics. Thismethodology has the main advantages of allowing to significantly reduce the computational complexity ofthe optimization process, especially when considering oscillator networks of large size, yielding to an efficientcontrol calculation.At this regard, we shall mention that GD methodologies have already been applied in the contextof the Kuramoto model. For instance, in [38], the author develop GD algorithms to efficiently solveoptimization problems that aim to maximize phase synchronization via network modifications. Moreover, in[22], optimization and control theory techniques are applied to investigate the synchronization properties of ageneralized Kuramoto model in which each oscillator lives on a compact, real Stiefel manifold. Nevertheless,to the best of our knowledge, the employment of stochastic techniques such as RBM to improve the efficiencyof the GD strategy has never been proposed in the context of the Kuramoto model.For completeness, let us stress that stochastic approaches have been widely considered, especially by themachine learning community, for treating minimization problems depending on very large data set. In thiscontext, they have shown amazing performances in terms of the computational efficiency (see, for instance,[4] and the references therein). Nowadays, stochastic techniques are among the preeminent optimizationmethods in fields like empirical risk minimization ([34]), data mining ([39]) or artificial neural networks ([31]).This contribution is organized as follows: in Section 2, we present the Kuramoto model and we discusssome of its more relevant properties. We also provide there a rigorous mathematical characterization of thesynchronization phenomenon. In Section 3, we introduce the controlled Kuramoto model and we describethe GD methodology for the control computation. Moreover, we briefly present the RBM approach and itsinclusion into the GD algorithm. Section 4 is devoted to the numerical simulations and to the comparison ofthe two optimization techniques considered in this paper. Finally, in Section 5 we summarize and discuss ourresults. 2. The mathematical model
From a mathematical viewpoint, synchronization phenomena are most often described through the so-calledKuramoto model, consisting of a population of N ≥ ˙ θ i ( t ) = ω i + KN N (cid:88) j =1 sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) , i = 1 , . . . , N, t > θ i (0) = θ i , (2.1)where θ i ( t ), i = 1 , . . . , N , is the phase of the i -th oscillator, ω i is its natural frequency and K is the couplingstrength.The frequencies ω i are assumed to be distributed with a given probability density f ( ω ), unimodal andsymmetric around the mean frequency Ω = 1 N N (cid:88) i =1 ω i , that is, f (Ω + ω ) = f (Ω − ω ).In this framework, each oscillator tries to run independently at its own frequency, while the coupling tendsto synchronize it to all the others. n the literature, many notions of synchronization have been considered. For identical oscillators (i.e., thosein which ω i = (cid:98) ω for every i = 1 , . . . , N ), one often studies whether the network can reach a configuration inwhich all the phases converge to the same value, that islim t → + ∞ | θ i ( t ) − θ j ( t ) | = 0 , for all i, j = 1 , . . . , N. (2.2)For systems with heterogeneous dynamics, such as when the natural frequencies ω i are not all identical(which is typical in real-world scenarios), this definition of synchronization is too restrictive (see [37]). Inthese cases, (2.2) is replaced by the alignment conditionlim t → + ∞ | ˙ θ i ( t ) − ˙ θ j ( t ) | = 0 , for all i, j = 1 , . . . , N, (2.3)according to which synchronization occurs when the phase differences given by | θ i ( t ) − θ j ( t ) | becomeconstant asymptotically for all i, j ∈ , . . . , N . This notion (2.3), which in some references is called completesynchronization (see for instance [15]), is the one that we will consider in this work.In its original work [21], Kuramoto considered the continuum limit case where N → + ∞ and showed thatthe coupling K has a key role in determining whether a network of oscillators can synchronize. In moredetail, he showed that, when the coupling K is weak, the oscillators run incoherently, whereas beyond acertain threshold collective synchronization emerges.Later on several research works provided specific bounds for the threshold of K ensuring synchronization(see, e.g., [1, 7, 9, 11, 12, 19]). In particular, in order to achieve (2.3) it is enough that K > K ∗ = | ω max − ω min | , (2.4)where ω min < ω max are the minimum and maximum natural frequencies.Notice, however, that (2.3) is an asymptotic characterization, meaning that is satisfied as t → + ∞ . Inthis work we are rather interested in the possibility of synchronizing the oscillators in a finite time horizon T .As we will discuss in the next section, this may be achieved by introducing a control into the Kuramotomodel (2.1). 3. Optimal control of the Kuramoto model
As we mentioned in Section 2, in this work we are interested in the finite-time synchronization of theKuramoto model. In particular, we aim at designing a control capable to steer the Kuramoto dynamics (2.1)to synchronization in a final time horizon T . In other words, we are going to consider the controlled system ˙ θ i ( t ) = ω i + Ku ( t ) N N (cid:88) j =1 sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) , i = 1 , . . . , N, t > θ i (0) = θ i , (3.1)and we want to compute a control function u such that the synchronized configuration (2.3) is achieved attime T , i.e., | ˙ θ i ( T ) − ˙ θ j ( T ) | = 0 , for all i, j = 1 , . . . , N. (3.2)From the practical applications viewpoint, this problem may be assimilated for instance to the necessityof synchronizing all the components of an electric grid after a black-out. In this interpretation, the differentelements in the grid are represented by the oscillators in (3.1), and T is the time horizon we provide for theblack-start, being therefore an external input to our problem. The objective is then to complete restoring thenetwork in a finite (possibly small) time T , which can be done by introducing a control in the system.To compute this optimal control allowing us to reach the synchronized configuration (3.2) we will adopt aclassical optimization approach based on the resolution of the following optimization problem (cid:98) u = min u ∈ L (0 ,T ; R ) J ( u ) J ( u ) = 12 N (cid:88) i,j =1 sin (cid:0) θ j ( T ) − θ i ( T ) (cid:1) + β (cid:107) u (cid:107) L (0 ,T ; R ) , (3.3) ubject to the dynamics (3.1). Here, with L (0 , T ; R ) we denoted the space of all functions u : (0 , T ) → R forwhich the following norm is finite: (cid:107) u (cid:107) L (0 ,T ; R ) := (cid:32)(cid:90) T | u ( t ) | dt (cid:33) . (3.4)In what follows, we will use the abridged notation (cid:107) u (cid:107) := (cid:107) u (cid:107) L (0 ,T ; R ) .In the cost functional (3.3), the first term enhances the fact that all the oscillators have to synchronize attime T . In particular, the optimal control (cid:98) u will yield to a dynamics in whichsin( θ j ( T ) − θ i ( T )) = 0 ⇒ θ j ( T ) − θ i ( T ) = kπ, k ∈ Z . (3.5)This is consistent with (3.2). For completeness, let us also stress that, in the case of identical oscillators,it has been shown for instance in [14] that, at least asymptotically, the two notions (3.2) and (3.5) coincide.The second term in (3.3) is introduced to avoid controls with a too large size. In it, β > (cid:98) u . Roughly speaking, thesmaller is β the larger will be (cid:98) u . A more detailed discussion on this point will be presented in Section 4.Through the minimization of J ( u ), we will obtain a unique scalar control function (cid:98) u : (0 , T ) → R , (cid:98) u > K violates the condition (2.4) and the uncontrolled dynamics runs incoherently towards adesynchronized configuration.Nevertheless, our proposed methodology may have the disadvantage of being less flexible than the otherswe mentioned above. In particular, it does not allow to control only a specific component of the network andthis may be a limitation in certain practical applications.Let us stress that the above considerations are merely heuristic and should be corroborated by a deeperanalysis based, for instance, on sharp numerical experiments. Notwithstanding that, in the present work wewill not address this specific issue, since our main interest is not to compare the performances of differentcontrol strategies but rather to present an efficient way to tackle the control problem (3.1).In the optimization literature, several different techniques have been proposed for minimizing the functional J ( u ) (see, e.g., [26]). In this work, we focus on the standard GD method, which looks for the minimum u asthe limit k → + ∞ of the following iterative process u k +1 = u k − η k ∇ J ( u k ) , (3.6)where η k > k the optimization scheme (3.6) requires tosolve (3.1), that is, a N -dimensional non-linear dynamical system. This may rapidly become computationallyvery expensive, especially when the number N of oscillators in our system is large.In order to reduce this computational burden, in this work we propose a novel methodology which combinesthe standard GD algorithm with the so-called Random Batch Method (RBM).RBM is a recently developed approach which has been introduced in [20] for the numerical simulationof high-dimensional collective behavior problems. This method uses small but random batches for particle nteractions, lowering the computational cost O ( N ) per time step to O ( N ), for systems with N particleswith binary interactions. Therefore, as our numerical simulations will confirm, embedding RBM into theGD iterative scheme yields to a less expensive algorithm and, consequently, to a more efficient controlcomputation. In what follows, we will call this approach GD-RBM algorithm, to differentiate it from thestandard GD one.3.1. The Gradient Descent approach.
Let us now describe in detail the GD approach to minimize thefunctional (3.3), and discuss its convergence properties.In order to fully define the iterative scheme (3.6), we need to compute the gradient ∇ J ( u ). Since we aredealing with a non-linear control problem, we will do this via the so-called Pontryagin maximum principle (see [41, Chapter 4, Section 4.8] or [40, Chapter 7]).To this end, let us first rewrite the dynamics (3.1) in a vectorial form as follows (cid:40) ˙Θ( t ) = Ω + F (cid:0) Θ( t ) , u ( t ) (cid:1) , t > , (3.7)with Θ := ( θ , . . . , θ N ) (cid:62) , Θ := ( θ , . . . , θ N ) (cid:62) and Ω := ( ω , . . . , ω N ) (cid:62) , and where F is the vector field givenby F = ( F , . . . , F N ) , F i := Ku ( t ) N N (cid:88) j =1 sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) , i = 1 , . . . , N. (3.8)In the control literature, (3.7) is usually called the primal system .Using the notation just introduced, we can then see that J ( u ) can be rewritten in the form J ( u ) = (cid:90) T L ( u ( t )) dt + φ (Θ( T )) , (3.9)with L ( u ( t )) = β | u ( t ) | and φ (Θ( T )) = 12 N (cid:88) j =1 sin (cid:0) θ j ( T ) − θ i ( T ) (cid:1) . Let us stress that (3.9) is in the standard form to apply the Pontryagin maximum principle. Through thisapproach, we can obtain the following expression for the gradient of J ( u ) ∇ J ( u ) = βu + ( D u F ) (cid:62) p, (3.10)where D u F indicates the Jacobian of the vector field F , computed with respect to the variable u .In (3.10), we denoted with p = ( p , . . . , p N ) the solution of the adjoint equation associated with (3.1),which is given by (cid:40) − ˙ p = ( D Θ F ) (cid:62) pp ( T ) = ∇ Θ( T ) φ (Θ( T )) , (3.11)where D Θ F stands again for the Jacobian of the vector field F , this time computed with respect to thevariable Θ.Taking into account the expression (3.8) of the vector field F , we can then readily check that the iterativescheme (3.6) becomes u k +1 = u k − η k βu k + KN N (cid:88) i =1 p i N (cid:88) j =1 sin( θ j − θ i ) , (3.12) ith − ˙ p i = − Kup i N N (cid:88) i (cid:54) = j =1 cos (cid:0) θ j − θ i (cid:1) + KuN N (cid:88) i (cid:54) = j =1 p j cos (cid:0) θ j − θ i (cid:1) , i = 1 , . . . , N, t > p i ( T ) = 12 N (cid:88) i (cid:54) = j =1 sin (cid:0) θ i ( T ) − θ j ( T ) (cid:1) . (3.13)In view of the above computations, the GD algorithm for the minimization of the cost functional J ( u ) canbe explicitly formulated as follows:GD algorithm. input Θ : initial condition of the primal system (3.7) u : initial guess for the control uk ←
0: iteration counter k max : maximum number of iterations allowed tol : tolerance while STOP-CRIT and k < k max do k ← k + 1 for j = 1 to N do Solve the the primal system (3.7)Solve the the adjoint system (3.13) end for
Update the control through the scheme (3.12) end whilereturn u k +1 = (cid:98) u : minimum of the functional J ( u ).In particular, we see that the control computation through the above algorithm requires, at each iteration k , to solve 2 N non-linear differential equations ( N for the variables θ i and N for p i ). If we introduce thetime-mesh of N t points0 = t < t < . . . < t N t = T, t m = t + m TN t , m = 1 , . . . , N t , at each time-step t m this operation has a computational cost of O ( N ) and the total computational complexityfor the simulation of (3.1) and (3.13) will then be O ( N t N ). If N is large, that is, if the number of oscillatorsin the network is considerable, this will rapidly become very expensive.3.2. The Random Batch Method.
In order to reduce the computational burden of GD for the optimizationprocess (3.3), we propose a modification of this algorithm which includes the aforementioned
Random BatchMethod (RBM) for the numerical simulation of the ODE systems (3.1) and (3.13).This technique, presented in [20] for interacting particle systems, is based on the following simple idea: ateach time step t m = m · dt in the mesh we employ to solve the dynamics, we divide randomly the N particlesinto n small batches with size 2 ≤ P < N , denoted by C q , q = 1 , . . . , n , that is C q = { i q , . . . , i q P } ⊂ { , . . . , N } , for all q = 1 , . . . , nC q ∩ C r = ∅ , for all q, r = 1 , . . . , n n (cid:91) q =1 C q = { , . . . , N } . Notice that the last batch C n may have size smaller than P if nP (cid:54) = N .Once this partition of { , . . . , N } has been performed, we solve the dynamics by interacting only particleswithin the same batch. This gives the following algorithm for the numerical approximation of (3.1) and(3.13): BM algorithm. for m = 1 to N t = T /dt do Divide randomly { , . . . , N } into n batches C q , q = 1 , . . . , n for q = 1 to n do Update θ i ( i ∈ C q ) by solving the ODE ˙ θ i = ω i + KuP (cid:88) j ∈ C q sin (cid:0) θ j − θ i (cid:1) θ i (0) = θ i . Update p i ( i ∈ C q ) by solving the ODE − ˙ p i = − Kup i P (cid:88) j ∈ C q cos (cid:0) θ j − θ i (cid:1) + KuP (cid:88) i (cid:54) = j ∈ C q p j cos (cid:0) θ j − θ i (cid:1) p i ( T ) = 12 (cid:88) i (cid:54) = j ∈ C q sin (cid:0) θ i ( T ) − θ j ( T ) (cid:1) . end forend for Regarding the complexity, note that random division into n batches of can be implemented using randompermutation. In Matlab, this can be done by using the function randperm(N) . Then, the first P elements areconsidered to be in the first batch, the second P elements are in the second batch, and so on. According to thediscussion presented in [20], at each time step t m this procedure yields to a cost of O ( P N ) for approximatingthe dynamics with RBM.If one is to simulate up to time T , the total number of time steps is N t as in the algorithm above. Then,the total computational complexity for the simulation of (3.1) and (3.13) is O ( P N t N ). Notice that, since P < N , this is always smaller than O ( N t N ).Summarizing, with the GD and GD-RBM methodologies we obtain the following per-iteration costs: • GD −→ cost GD = C GD N t N . • GD-RBM −→ cost GD − RBM = C GD − RBM
P N t N .Therefore, independently of the value of N , employing RBM to simulate the dynamics in each iteration ofGD yields improvements in terms of the computational cost.For completeness, we shall mention that the above considerations are simply heuristic and would require adeeper analysis. As a matter of fact, to have a rigorous validation of the reduction in the computationalcomplexity when using RBM one should have more precise information on the two constants C GD and C GD − RBM , and be sure that the difference among them does not overwhelm the help that the batchingprocedure is providing.At this regard, let us stress that the RBM method has been developed only recently in [20] and, at presenttime, there is not a well-established qualitative analysis on its computational cost, going in more detail thanwhat we mentioned above. The evidence that in our case of the Kuramoto model (3.1) the GD-RBM methodallows for a more efficient control computation, in particular for large oscillator networks, will then be giventhrough the numerical simulations in Section 4.3.3.
Convergence analysis.
To complete this section, let us briefly comment about the convergenceproperties of the GD methodology. It is nowadays classically known that the convergence rate of the GDalgorithm is determined by the regularity of the objective function. In our case, since J ( u ) is L-smooth, thatis (cid:107)∇ J ( u ) − ∇ J ( v ) (cid:107) ≤ L (cid:107) u − v (cid:107) , it can be proven that (cid:107)∇ J ( u k ) (cid:107) → k → + ∞ (3.14) nd (cid:107) J ( u k ) − J ( (cid:98) u ) (cid:107) = O (cid:18) k (cid:19) , (3.15)where, we recall, (cid:98) u denotes the minimum of J ( u ) and the norm (cid:107) · (cid:107) L (0 ,T ; R ) has been defined in (3.4).In particular, (3.15) implies that for achieving ε -optimality, i.e. for obtaining (cid:107) J ( u k ) − J ( (cid:98) u ) (cid:107) < ε , theGD algorithm requires k = O ( ε − ) iterations.Combining this with the per-iteration costs we gave at the end of Section 3.2, we can thus obtain thefollowing total computational costs • GD: O (cid:16) N t N ε (cid:17) • GD-RBM: O (cid:0) P N t Nε (cid:1) ,and we can conclude that the GD-RBM approach will be more efficient than the standard GD one to solveour optimization problem. This is enhanced for large values of N and will be confirmed by our numericalsimulations. 4. Numerical simulations
We present here our numerical results for the control of N coupled oscillators described by the Kuramotomodel (3.1), following the strategy previously described. This section is divided into two parts:1. In a first moment, we will show that the optimization problem (3.3) indeed allows to compute aneffective control function which is capable to steer the Kuramoto model (3.1) to a synchronizedconfiguration. This will be done both for a strong coupling K > K ∗ (see (2.4)) and for a weakcoupling K < K ∗ . Besides, we will also briefly analyze the role of the parameter β in the optimizationprocess. Finally, we will show the efficacy of our control strategy in the more realistic cases of asparse interaction network and for a second-order Kuramoto model with damping.2. Once the effectiveness of the control strategy we propose has been corroborated, we will comparethe GD and GD-RBM algorithms for the minimization of J ( u ). In particular, we will show how theRBM approach allows to significantly reduce the computational complexity of the GD algorithm forthe calculation of the control u , especially when considering oscillator networks of large dimension.The oscillators are chosen such that their natural frequencies are given following the normal probabilitylaw f ( ω ) = 15 σ √ π e − ω σ , (4.1)with σ = 0 .
1. This means that the values of ω min and ω max are given respectively by ω min = min ω ∈ R f ( ω ) = 0 ω max = max ω ∈ R f ( ω ) = f (0) = 2 √ π and the coupling gain K which is necessary for synchronization in the absence of a control has to satisfy (see(2.4)) K > | ω max − ω min | = 2 √ π . The initial datum θ is chosen following a normal distribution as well, in such a way that | θ i − θ j | < π for all i, j = 1 , . . . , N . Let us stress that this choice of θ allows the synchronization of the uncontrolledmodel, as it has been shown for instance in [10].Without loss of generality, we considered the time horizon T = 3 s for completing the synchronization.That is, we want all the oscillators in our model to reach the configuration (3.2) in three seconds.Finally, we used an explicit Euler scheme for solving the direct and adjoint dynamics (3.1) and (3.13)during the minimization of J ( u ), and we chose as a stopping criterion e k := (cid:107)∇ J ( u k ) (cid:107) (cid:107) u k (cid:107) < ε, (4.2) ith ε = 10 − , and where the notation (cid:107) · (cid:107) stands again for the L (0 , T R )-norm defined in (3.4). Let usstress that the stopping criterion (4.2) is consistent with (3.14) and (3.15).4.1. Computation of the optimal control.
In this section, we show that through the optimizationproblem (3.3) we are able to compute an effective control function which is capable to steer the Kuramotomodel (3.1) to a synchronized configuration in a given time horizon T .We performed the simulations in Matlab R2018a on a laptop with Intel Core i − U CP U @2 . GHz × . N = 10 oscillators in an all-to-all coupled configuration andwith a coupling gain K > K ∗ . Moreover, we set the penalization parameter β in (3.3) to take the value β = 10 − When using the GD-RBM approach, the family of N = 10 oscillators has been separated into n = 5batches of size P = 2.In Figure 1-top, we show the evolution of the uncontrolled dynamics, which corresponds to taking u ≡ K . Nevertheless, synchronization is not reached in the short time horizonwe are providing. At this regard, let us remark that, for the uncontrolled Kuramoto model with a sufficientlystrong coupling gain, synchronization is expected to be reached only asymptotically, i.e. when t → + ∞ .In Figure 1-bottom, we show the evolution of the same dynamics, this time under the action of the controlfunction u computed through the minimization of J ( u ). The subplot on the left corresponds to the simulationsdone with the GD approach, while the one on the right is done employing GD-RBM. We can clearly see how,in both cases, the oscillators are all synchronized at the final time T = 3 s . This means that both algorithmsmanaged to compute an effective control. Figure 1.
Top: evolution of the free dynamics of the Kuramoto model (3.1) with N = 10oscillators. Bottom: evolution of the controlled dynamics of the Kuramoto model (3.1) with N = 10 oscillators. The control function (cid:98) u is obtained with the GD (left) and the GD-RBM(right) approach.In Figure 2, we show the convergence of the error in logarithmic scale when applying both the GD andGD-RBM approach. We can appreciate how, in the case of GD-RBM, this convergence is not monotonicas it is for the GD algorithm. This, however, is not surprising due to the stochastic nature of the RBMmethodology. igure 2. Convergence of the error e k (see (4.2)) in logarithmic scale with the GD (left)and GD-RBM (right) algorithm.In Figure 3, we display the behavior of the control function (cid:98) u computed via the GD-RBM algorithm. Wecan see how, at the beginning of the time interval we are considering, this control is close to one and it isincreasing with a small slope. On the other hand, this growth becomes more pronounced as we get closer tothe final time T = 3 s . Figure 3.
Control function (cid:98) u obtained through the GD-RBM algorithm applied to theKuramoto model (3.1) with N = 10 oscillators.Notice that, in (3.1), (cid:98) u enters as a multiplicative control which modifies the strength of the coupling K .Hence, according to the profile displayed in Figure 2, the control function (cid:98) u we computed is initially lettingthe system evolving following its natural dynamics. Then, as the time evolves towards the horizon T = 3 s , (cid:98) u enhances the coupling strength K in order to reach the desired synchronized configuration (3.2).Finally, notice also that the control (cid:98) u is always positive. This is actually not surprising, if one takes intoaccount the following observation.In the Kuramoto model, in order to reach synchronization the coupling strength K needs to be positive.Otherwise, the system would converge to a desynchronized configuration (see [18]). Moreover, according tothe model (3.1), if we start from K >
0, in order to maintain this coupling positive (cid:98) u has to remain positiveas well.Recall that (cid:98) u is computed minimizing the functional (3.3), in which the second term is a measurement ofthe level of synchronization in the model. Hence, since negative values of (cid:98) u would lead to desynchronizationand to the corresponding increasing of the functional, these values remain automatically excluded during theminimization process.Let us now discuss briefly the role of the penalization parameter β in the computation of the optimalcontrol. To this end, we have run simulations with different values of β = 10 − , − , − and 10 − .As we already mentioned in Section 3, in the cost functional (3.3) β is a (usually small) penalizationparameter which allows to tune the norm of the optimal control (cid:98) u , that is, the amount of energy that thecontrol introduces into the system. Roughly speaking, the smaller is β the larger will be (cid:98) u .This is clearly seen in Figure 4. In particular, we can appreciate how, for β = 10 − , the computed controlremains smaller than in the other cases. igure 4. Control function (cid:98) u obtained through the GD-RBM algorithm applied to theKuramoto model (3.1) with N = 10 oscillators and different values of β .We already mentioned above that the effect of the control in (3.1) is to enhance synchronization bymodifying the strength of the coupling K . Hence, we can expect that, if (cid:98) u is small (in particular, if it remainsclose to one), the synchronization properties of the Kuramoto model (3.1) will be worst than when applyinga larger control. At this regard, let us recall that the level of synchronization in (3.1) can be analyzed interms of the quantity r ( t ) := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N N (cid:88) j =1 e iθ j ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , measuring the coherence of the oscillator population (see [1]). In particular we always have 0 ≤ r ( t ) ≤ r reaches the value one.In Figure 5, we show the behavior of r ( t ) with respect to the parameter β . On the one hand, in all thecases displayed we can clearly see that r ( T ) = 1. This means that all the computed controls will be effectiveto steer the system (3.1) to its synchronized configuration at time T . On the other hand, we can also noticehow, when decreasing β , the function r ( t ) reaches the value one faster, meaning that the correspondingcontrol is expected to yield to better synchronization properties. Figure 5.
Behavior of the synchronization function r ( t ) corresponding to the controlleddynamics (3.1) for different values of the parameter β .Let us now conclude this section by showing that the control strategy that we propose in this paperis effective also in situations in which the coupling gain among the oscillators is too weak to ensuresynchronization for the uncontrolled dynamics of the Kuramoto model. This corresponds to taking K < K ∗ (see (2.4)). In particular, we will consider the case K < t increases. On the other hand, when applying the control (cid:98) u , the system is once again steered to asynchronized configuration. igure 6. Evolution of the uncontrolled (left) and controlled (right) dynamics of theKuramoto model (3.1) with N = 10 oscillators and K < K so that the oscillators are all synchronized at time T and that, for the uncontrolleddynamics, synchronization requires K > Figure 7.
Control function (cid:98) u obtained through the GD-RBM algorithm applied to theKuramoto model (3.1) with N = 10 oscillators and K <
The case of a sparse interaction network.
We start by considering the case of a sparse interactionnetwork in our Kuramoto model (3.1). In other words, we are considering here the following system ˙ θ i ( t ) = ω i + Ku ( t ) N N (cid:88) j =1 a i,j sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) , i = 1 , . . . , N, t > θ i (0) = θ i , (4.3)with a i,j = (cid:40) , if θ i is connected with θ j , if θ i is not connected with θ j A schematic representation of the network considered in our simulations is given in Figure 8, in which theblue dots correspond to a i,j = 1.The simulations have been performed with the same initial datum and time horizon we considered in ourprevious experiments. Moreover, we addressed here only the case of a strong coupling gain K > K ∗ . Theminimization of the functional J ( u ) has been performed with the GD algorithm.In Figure 9, we show the evolution of the uncontrolled and controlled dynamics. As we can see, while inthe absence of a control the oscillators are evolving towards a desynchronized configuration, when applyingthe control function we computed the system still reaches synchronization at time T . igure 8. Sparse interaction scheme for the Kuramoto model (4.3).
Figure 9.
Evolution of the uncontrolled (left) and controlled (right) dynamics of theKuramoto model (4.3) with N = 10 oscillators, K > K ∗ and interactions as in Figure 8.The control function obtained for these numerical experiments is plotted in Figure 10. We can observe how,differently from what is shown in Figures 3 and 4, this time (cid:98) u reaches larger values. This is not surprising, ifwe consider that now our model has a lower level of interactions and if we recall our previous discussion onhow our control affects the Kuramoto dynamics. Figure 10.
Control function (cid:98) u for the Kuramoto model (4.3) with N = 10 oscillators, K > K ∗ and interactions as in Figure 8. .1.2. A second-order model with damping.
We consider here the second-order Kuramoto model with damping ¨ θ i ( t ) + ˙ θ i ( t ) = ω i + Ku ( t ) N N (cid:88) j =1 sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) , i = 1 , . . . , N, t > θ i (0) = θ i , ˙ θ i (0) = θ i , which can be rewritten as the following first-order system ˙ θ i ( t ) = ξ i ( t ) , i = 1 , . . . , N, t > ξ i ( t ) = − ξ i ( t ) + ω i + Ku ( t ) N N (cid:88) j =1 sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) , i = 1 , . . . , N, t > θ i (0) = θ i , ξ i (0) = θ i . (4.4)In the context of power grids, this model has been firstly introduced in the work [13]. Later on, in [32], ithas been extended to consider more complex scenarios in which the dynamics of the voltage amplitude istaken into account.Also in this case, we are interested in computing a control capable to steer the system to the synchronizedconfiguration (3.2). This can be done once again by solving the optimal control problem (3.3), this timeunder the dynamics (4.4).The simulations have been performed with the same initial datum Θ = ( θ i ) Ni =1 we considered in ourprevious experiments and with Θ = ( θ i ) Ni =1 = (0 , , . . . , (cid:62) . The time horizon is once again T = 3 s .Moreover, we addressed here both the cases of a strong coupling gain K > K ∗ and of a negative one K < J ( u ) has been performed with the GD algorithm.At this regard, we shall mention that, in [43], the GD methodology has been applied to obtain synchro-nization in a sparse network of Kuramoto oscillators with damping under the action of an additive control,i.e. the following model ˙ θ i ( t ) = ξ i ( t ) , i = 1 , . . . , N, t > ξ i ( t ) = − ξ i ( t ) + ω i + KN N (cid:88) j =1 a i,j sin (cid:0) θ j ( t ) − θ i ( t ) (cid:1) + u i ( t ) , i = 1 , . . . , N, t > θ i (0) = θ i , ξ i (0) = θ i . The advantages and disadvantages of this additive control action with respect to the multiplicative onewe propose have been discussed in Section 3. In particular, our control strategy allows us to deal also withnegative coupling gain
K <
0, while in [43] only
K >
K > K ∗ and K <
0, respectively. Also in this case, the proposed control strategy allows us to compute aneffective control function (cid:98) u which steers the system to a synchronized configuration in time T .Finally, the control functions obtained for these numerical experiments are plotted in Figure 12. Also inthese cases we can observe different behaviors than what is shown in Figures 3 and 4. In particular, this time (cid:98) u changes sign in the time horizon (0 , T ). At this regard, let us mention that for the second-order Kuramotomodel (4.4), due to the presence of the damping term, our previous considerations on the sign of the controldo not apply anymore. Hence, it is not surprising to obtain a behavior as the one displayed.4.2. Comparison of GD and GD-RBM.
We conclude this section on the numerical experiments bycomparing the performances of the GD and GD-RBM algorithms for the computation of the optimal control (cid:98) u . To this end, we run simulations for increasing values of N , namely N = 10 , , , , T = 3 s and a penalization parameter β = 10 − . Moreover, we considered the case ofa large coupling gain K > K ∗ .In what follows, we focus only on the simple case of the first-order Kuramoto model (3.1) with an all-to-allinteraction network. We will briefly comment about possible extensions to the more realistic scenariosdescribed in Sections 4.1.1 and 4.1.2 in the the last part of this paper, devoted to conclusions and openproblems. igure 11. Evolution of the uncontrolled (left) and controlled (right) dynamics of theKuramoto model (4.4) with N = 10 oscillators and K > K ∗ (top) and K <
Figure 12.
Control function (cid:98) u for the Kuramoto model (4.4) with N = 10 oscillators and K > K ∗ (left) and K < N , the two approaches show similar behaviors. Nevertheless,when increasing the number of oscillators in our system, the advantages of the GD-RBM methodologywith respect to GD become evident. In particular, the growth of the computational time for GD-RBM issignificantly less pronounced than for GD. As a matter of fact, in the case of N = 1000 oscillators, we decidednot to perform the simulations with the GD algorithm since its behavior with smaller values of N alreadysuggested that this experiment would be computationally too expensive.On the other hand, even with N = 1000 oscillators in the system, the GD-RBM approach turns out to beable to compute an effective control for the Kuramoto model (3.1) (see Figure 14) in about 29 seconds. D GD-RBM N Time (sec.) Time (sec.)10 5 . .
450 11 . . . . . . . Table 1.
Computational times required by the GD and GD-RBM algorithm to computethe optimal control (cid:98) u with increasing values of N . Figure 13.
Computational times required by the GD and GD-RBM algorithm to computethe optimal control (cid:98) u with increasing values of N . Figure 14.
Evolution of the controlled dynamics of the Kuramoto model (3.1) with N =1000 oscillators. The control has been computed with the GD-RBM algorithm.5. Conclusions
This paper deals with the synchronization of coupled oscillators described by the Kuramoto model. Inparticular, we design a unique scalar control function u ( t ) capable of steering a N -dimensional network ofoscillators to a synchronized configuration in a finite time horizon. This is done following a standard optimalcontrol approach and obtaining the function u ( t ) via the minimization of a suitable cost functional. ith this approach, we computed a control which acts as a multiplicative force enhancing the couplingbetween the oscillators in the network, thus favoring the synchronization in the time horizon provided.To carry out this minimization process, we used a Gradient Descent (GD) methodology, commonlyemployed in the optimal control community, which we have coupled with the novel Random Batch Method(RBM) for a more efficient numerical resolution of the Kuramoto dynamics. The main purpose of this workhas been to show how the introduction of RBM into GD may yield considerable improvement in term of thecomputational complexity, in particular for large oscillator networks.Our simulations results have shown the following main facts: • The proposed control strategy is indeed effective to reach a synchronized configuration in a finite timehorizon. Moreover, it allows to deal efficiently with large networks of oscillators (namely, N = 1000)and with the case of a low coupling gain in the network, when the uncontrolled dynamics is notexpected to reach a synchronized configuration. • For large values of N , the inclusion of RBM into the GD algorithm significantly reduces thecomputational burden for the control computation, thus allowing to deal with high-dimensionaloscillators networks in a more efficient way.In conclusion, the study conducted in this paper suggests that the proposed methodology, based on thecombination of the standard GD optimization algorithm with the novel RBM method for the numericalresolution of multi-agent dynamics, may help in significantly reduce the computational complexity for thecontrol computation of the Kuramoto model, in particular in the case of a high-dimensional system.At this regard, we shall stress that the analysis in this paper has been developed mostly in a simplifiedframework of a network with all-to-all coupling, in which the oscillators are all interacting among themselves.The following interesting questions remain unaddressed: • To study whether our methodology remains valid in more complex scenarios of networks with sparseinteraction topologies, perturbations due to disconnections, or rewiring. In particular, it would berelevant to determine if the reduction in computational complexity that we obtained through theGD-RBM algorithm is related with the density of the network or, instead, this approach allows todeal efficiently also with the case of a low number of interactions. A starting point for this analysiswould be to determine whether the GD-RBM methodology can be successfully applied in the scenarioswe addressed in Sections 4.1.1 and 4.1.2. • To analyze what happens if, instead of selecting the oscillators uniformly at random during thebatching process in RBM, we organize them in groups with similar frequencies. At this regard, itwould be important to understand if the methodology we propose is still effective or, instead, somemodifications need to be introduced. • To analyze if our methodology may be applied also in different framework than the ones consideredin this paper. For instance, a similar formalism than the one we have proposed has been developedfor computing rare-events in oscillator networks due to noise. See for instance [17, 16]. In thesementioned contributions, the objective functional is the probability for a rare event to happen. Theactions are more complicated than the simple L -norm, but the batch techniques we employed maybe useful in this context as well.All these are key open problem which will be considered in future works. Acknowledgment
The authors wish to acknowledge Jes´us Oroya and Dongnam Ko (Chair of Computational Mathematics,Fundaci´on Deusto, Bilbao, Spain) for interesting discussion on the topics of the present paper.
References [1]
Acebr´on, J. A., Bonilla, L. L., Vicente, C. J. P., Ritort, F., and Spigler, R.
The Kuramoto model: A simpleparadigm for synchronization phenomena.
Rev. Modern Phys. 77 , 1 (2005), 137.[2]
Ben-Naim, E.
Opinion dynamics: rise and fall of political parties.
Europhys. Lett. 69 , 5 (2005), 671.[3]
Biccari, U., Ko, D., and Zuazua, E.
Dynamics and control for multi-agent networked systems: A finite-differenceapproach.
Math. Models Methods Appl. Sci. 29 , 04 (2019), 755--790.[4]
Bottou, L., Curtis, F., and Nocedal, J.
Optimization methods for large-scale machine learning.
SIAM Rev. 60 , 2 (2018).[5]
Buck, J.
Synchronous rhythmic flashing of fireflies II.
Quart. Rev. Bio. 63 , 3 (1988), 265--289. Chassin, D. P., and Posse, C.
Evaluating north american electric grid reliability using the Barab´asi-Albert network model.
Phys. A 355 , 2-4 (2005), 667--677.[7]
Chopra, N., and Spong, M. W.
On synchronization of Kuramoto oscillators. In
Proceedings of the 44th IEEE Conferenceon Decision and Control (2005), IEEE, pp. 3916--3922.[8]
Chopra, N., and Spong, M. W.
Passivity-based control of multi-agent systems. In
Adv. Robot Control . Springer, 2006,pp. 107--134.[9]
Chopra, N., and Spong, M. W.
On exponential synchronization of Kuramoto oscillators.
IEEE Trans. Autom. Control54 , 2 (2009), 353--357.[10]
Dong, J.-G., and Xue, X.
Synchronization analysis of Kuramoto oscillators.
Commun. Math. Sci. 11 , 2 (2013), 465--480.[11]
D¨orfler, F., and Bullo, F.
Synchronization of power networks: Network reduction and effective resistance.
IFACProceedings Volumes 43 , 19 (2010), 197--202.[12]
D¨orfler, F., Chertkov, M., and Bullo, F.
Synchronization in complex oscillator networks and smart grids.
Proc. Nat.Acad. Sci. 110 , 6 (2013), 2005--2010.[13]
Filatrella, G., Nielsen, A. H., and Pedersen, N. F.
Analysis of a power grid using a Kuramoto-like model.
Europ.Phys. J. B 61 , 4 (2008), 485--491.[14]
Ha, S.-Y., Kim, H. K., and Park, J.
Remarks on the complete synchronization of Kuramoto oscillators.
Nonlinearity 28 , 5(2015), 1441.[15]
Ha, S.-Y., Kim, H. K., and Ryoo, S. W.
Emergence of phase-locked states for the Kuramoto model in a large couplingregime.
Commun. Math. Sci. 14 , 4 (2016), 1073--1091.[16]
Hindes, J., Jacquod, P., and Schwartz, I. B.
Network desynchronization by non-gaussian fluctuations.
Phys. Rev. E100 , 5 (2019), 052314.[17]
Hindes, J., and Schwartz, I. B.
Rare slips in fluctuating synchronized oscillator networks.
Chaos 28 , 7 (2018), 071106.[18]
Hong, H., and Strogatz, S. H.
Kuramoto model of coupled oscillators with positive and negative coupling parameters:an example of conformist and contrarian oscillators.
Phys. Rev. Lett. 106 , 5 (2011), 054102.[19]
Jadbabaie, A., Motee, N., and Barahona, M.
On the stability of the Kuramoto model of coupled nonlinear oscillators.In
Proc. 2004 American Control Conference (2004), vol. 5, IEEE, pp. 4296--4301.[20]
Jin, S., Li, L., and Liu, J.-G.
Random Batch Methods (RBM) for interacting particle systems.
J. Comput. Phys. 400 (2020), 108877.[21]
Kuramoto, Y.
Self-entrainment of a population of coupled non-linear oscillators. In
International symposium on mathe-matical problems in theoretical physics (1975), Springer, pp. 420--422.[22]
Markdahl, J., Thunberg, J., and Goncalves, J.
High-dimensional Kuramoto models on Stiefel manifolds synchronizecomplex networks almost globally.
Automatica 113 (2020), 108736.[23]
Mehyar, M., Spanos, D., Pongsajapan, J., Low, S. H., and Murray, R. M.
Distributed averaging on asynchronouscommunication networks. In
Proceedings of the 44th IEEE Conference on Decision and Control (2005), IEEE, pp. 7446--7451.[24]
Nabi, A., and Moehlis, J.
Single input optimal control for globally coupled neuron networks.
J. Neural Eng. 8 , 6 (2011),065008.[25]
Nesterov, Y.
Applied optimization . Introductory Lectures on Convex Optimization: A Basic Course. Kluwer AcademicPublishers Boston, Dordrecht, London, 2004.[26]
Nocedal, J., and Wright, S.
Numerical optimization . Springer Science & Business Media, 2006.[27]
Olfati-Saber, R.
Flocking for multi-agent dynamic systems: Algorithms and theory.
IEEE Trans. Autom. Control 51 , 3(2006), 401--420.[28]
Olfati-Saber, R., Fax, J. A., and Murray, R. M.
Consensus and cooperation in networked multi-agent systems.
Proc.IEEE 95 , 1 (2007), 215--233.[29]
Rosenblum, M., and Pikovsky, A.
Delayed feedback control of collective synchrony: An approach to suppression ofpathological brain rhythms.
Phys. Rev. E 70 , 4 (2004), 041904.[30]
Sachtjen, M., Carreras, B., and Lynch, V.
Disturbances in a power transmission system.
Phys. Rev. E 61 , 5 (2000),4877.[31]
Schmidhuber, J.
Deep learning in neural networks: an overview.
Neural networks 61 (2015), 85--117.[32]
Schmietendorf, K., Peinke, J., Friedrich, R., and Kamps, O.
Self-organized synchronization and voltage stability innetworks of synchronous machines.
Eur. Phys. J. Spec. Top. 223 , 12 (2014), 2577--2592.[33]
Sepulchre, R., Paley, D., and Leonard, N.
Collective motion and oscillator synchronization. In
Cooperative control .Springer, 2005, pp. 189--205.[34]
Shalev-Shwartz, S., and Zhang, T.
Accelerated proximal stochastic dual coordinate ascent for regularized loss minimiza-tion. In
International Conference on Machine Learning (2014), pp. 64--72.[35]
Strogatz, S. H.
Exploring complex networks.
Nature 410 , 6825 (2001), 268--276.[36]
Strogatz, S. H., Abrams, D. M., McRobie, A., Eckhardt, B., and Ott, E.
Crowd synchrony on the millennium bridge.
Nature 438 , 7064 (2005), 43--44.[37]
Sun, J., Bollt, E. M., and Nishikawa, T.
Master stability functions for coupled nearly identical dynamical systems.
Europhys. Lett. 85 , 6 (2009), 60011.[38]
Taylor, D., Skardal, P. S., and Sun, J.
Synchronization of heterogeneous oscillators under network modifications:Perturbation and optimization of the synchrony alignment function.
SIAM J. Appl. Math. 76 , 5 (2016), 1984--2008.[39]
Toscher, A., and Jahrer, M.
Collaborative filtering applied to educational data mining.
KDD cup (2010). Tr´elat, E.
Contrˆole optimal: th´eorie & applications . Vuibert, Collection Math´ematiques Concr`etes, 2005.[41]
Tr¨oltzsch, F.
Optimal control of partial differential equations: theory, methods, and applications . American MathematicalSoc., 2010.[42]
Tukhlina, N., Rosenblum, M., Pikovsky, A., and Kurths, J.
Feedback suppression of neural synchrony by vanishingstimulation.
Phys. Rev. E 75 , 1 (2007), 011918.[43]
Tumash, L., Olmi, S., and Sch¨oll, E.
Stability and control of power grids with diluted network topology.
Chaos: AnInterdisciplinary Journal of Nonlinear Science 29 , 12 (2019), 123105.[44]
Walker, T. J.
Acoustic synchrony: two mechanisms in the snowy tree cricket.
Science 166 , 3907 (1969), 891--894.[45]
Wiesenfeld, K., Colet, P., and Strogatz, S. H.
Frequency locking in Josephson arrays: connection with the Kuramotomodel.
Phys. Rev. E 57 , 2 (1998), 1563.[46]
Winfree, A. T.
Biological rhythms and the behavior of populations of coupled oscillators.
J. Theor. Bio. 16 , 1 (1967),15--42.
Umberto Biccari, Chair of Computational Mathematics, Fundaci´on Deusto, Avda. de las Universidades 24,48007 Bilbao, Basque Country, Spain.Umberto Biccari, Universidad de Deusto, Avda Universidades 24, 48007 Bilbao, Basque Country, Spain.
E-mail address : [email protected],[email protected] Enrique Zuazua, Chair in Applied Analysis, Alexander von Humboldt-Professorship, Department of MathematicsFriedrich-Alexander-Universit¨at Erlangen-N¨urnberg, 91058 Erlangen, Germany.Enrique Zuazua, Chair of Computational Mathematics, Fundaci´on Deusto, Avda. de las Universidades 24, 48007Bilbao, Basque Country, Spain.Enrique Zuazua, Departamento de Matem´aticas, Universidad Aut´onoma de Madrid, 28049 Madrid, Spain.
E-mail address : [email protected]@fau.de