A Stochastic Approach to the Synchronization of Coupled Oscillators

This paper deals with an optimal control problem associated with the Kuramoto model describing the dynamical behavior of a network of coupled oscillators. Our aim is to design a suitable control function allowing us to steer the system to a synchronized configuration in which all the oscillators are aligned on the same phase. This control is computed via the minimization of a given cost functional associated with the dynamics considered. For this minimization, we propose a novel approach based on the combination of a standard 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 employment of RBM improves the performances of the GD algorithm, reducing the computational complexity of the minimization process and allowing for a more efficient control calculation.


Introduction
Synchronization is a common phenomenon which has been observed in biological, chemical, physical and social systems for centuries and has attracted the interest of researcher in a diversified spectrum of scientific fields.
Common examples of synchronization phenomena often cited in review articles include groups of synchronously chirping crickets ( [30]), fireflies flashing in unison ( [5]), the superconducting Josephson junction ( [31]), or crowds of people walking together that will tend to synchronize their footsteps ( [27]).
Roughly speaking, synchronization means that a network of several periodic processes with different natural frequencies reaches an equilibrium configuration sharing the same common frequency as a result of their mutual interaction.
This concept is closely related to the one of consensus for multi-agent systems, widely analyzed in many different frameworks including collective behavior of flocks and swarms, opinion formation, and distributed computing (see [2,3,15,19,20]). In broad terms, consensus means to reach an agreement regarding a certain quantity 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 stability analysis of utility power grids ( [6,10,22,26]). Indeed, large networks of connected power plants need to be synchronized 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 ( [14]), describing the dynamical behavior of a (large) network of oscillators in a all-to-all coupled configuration in which every oscillator is connected with all the others. This model extends the original studies by Winfree in the context of mutual synchronization in multi-oscillator systems based on a phase description ( [32]). In particular, in Kuramoto's work, synchronization appears as an asymptotic pattern which is spontaneously reached by the system when the interactions among the oscillators are sufficiently strong.
In some more recent contribution, control theoretic methods have been employed to analyze the synchronization phenomenon. For instance, in [8] the authors design passivity-based controls for the synchronization of multi-agent systems, with application to the general problem of multi-robot coordination. In [24], feedback control laws for the stabilization of coupled oscillators are designed and analyzed. Finally, [16] deals with the problem of desynchronizing a network of synchronized and globally coupled neurons using an input to a single neuron. This is done in the spirit of dynamic 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 construction of a suitable control function via an appropriate optimization process. To this end, we propose a novel approach which combines a standard Gradient Descent (GD) methodology with the recently-developed Random Batch Method (RBM, see [13]) for the efficient numerical approximation of collective dynamics. This methodology has the main advantages of allowing to significantly reduce the computational complexity of the optimization process, especially when considering oscillators' networks of large size, yielding to an efficient control calculation.
For completeness, let us stress that stochastic approaches have been widely considered, especially by the machine learning community, for treating minimization problems depending on very large data set. In this context, 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 optimization methods in fields like empirical risk minimization ( [25]), data mining ( [28]) or artificial neural networks ( [23]).
This contribution is organized as follows: in Section 2, we present the Kuramoto model and we discuss some of its more relevant properties. We also provide there a rigorous mathematical characterization of the synchronization phenomenon. In Section 3, we introduce the controlled Kuramoto model and we describe the GD methodology for the control computation. Moreover, we briefly present the RBM approach and its inclusion into the GD algorithm. Section 4 is devoted to the numerical simulations and to the comparison of the two control techniques considered in this paper. Finally, in Section 5 we summarize and discuss our results.

The mathematical model
From a mathematical viewpoint, synchronization phenomena are most often described through the so-called Kuramoto model, consisting of a population of N ≥ 2 coupled oscillators whose dynamics are governed by the following system of non-linear first-order ordinary differential equations where θ i (t), i = 1, . . . , N , is the phase of the i-th oscillator, ω i is its natural frequency and K > 0 is the coupling strength.
The frequencies ω i are assumed to be distributed with a given probability density f (ω), unimodal and symmetric around the mean frequency In this framework, each oscillator tries to run independently at its own frequency, while the coupling tends to synchronize it to all the others. The oscillators are said to be synchronized if In other words, synchronization occurs when the phase differences given by |θ i (t) − θ j (t)| become constant asymptotically for all i, j ∈ 1, . . . , N .
In its original work [14], Kuramoto considered the continuum limit case where N → +∞ and showed that the coupling K has a key role in determining whether a network of oscillators can synchronize. In more detail, he showed that, when the coupling K is weak, the oscillators run incoherently, whereas beyond a certain threshold collective synchronization emerges spontaneously.
Later on several research works provided specific bounds for the threshold of K ensuring synchronization (see, e.g., [1,7,9,12]). In particular, in order to achieve (2.2) it is enough that where ω min < ω max are the minimum and maximum natural frequencies.
Notice, however, that (2.2) is an asymptotic characterization, meaning that is satisfied as t → +∞. In this 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 Kuramoto model (2.1).

Optimal control of the Kuramoto model
As we mentioned in Section 2, in this work we are interested in the finite-time synchronization of the Kuramoto 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 and we want to compute a control function u such that the synchronized configuration (2.2) is achieved at time T , i.e., To this end, we will follow a classical optimization approach and we will compute our control via the following optimization problem u = min subject to the dynamics (3.1).
In the cost functional (3.3), the first term is a penalization one introduced to avoid controls with a too large size, while the second term enhances the fact that all the oscillators have to synchronize at time T . Moreover, with · L 2 (0,T ;R) we denoted the following norm: Through the minimization of J(u), we will obtain a unique scalar control function u for all the oscillators included in the network. In other words, we are going to define a unique control law which is capable to act on the global dynamics in order to drive it to the desired synchronized configuration.
Let us stress that an analogous problem has already been treated for instance in [16], where the authors have considered the problem of desynchronizing a system of globally coupled neuron by constructing a single input optimal control, which is thus acting only on one component of the system, through the minimization of a suitable cost functional.
In the optimization literature, several different techniques have been proposed for minimizing the functional J(u) (see, e.g., [18]). In this work, we focus on the standard GD method, which looks for the minimum u as the limit k → +∞ of the following iterative process where η k > 0 is called the step-size or, in the machine learning context, the learning rate. The step size is typically chosen to be a constant depending on certain key parameters of the optimization problem, or following an adaptive strategy. See, e.g., [17, Section 1.2.3] for more details.
This gradient technique is typically chosen because it is easy to implement and not very memory demanding. Nevertheless, as we shall see in more detail in the next section, when applied to the optimal control of collective dynamics the GD methodology has a main drawback.
Indeed, at each iteration k the optimization scheme (3.5) requires to solve (3.1), that is, a N -dimensional non-linear dynamical system. This may rapidly become computationally very 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 combines the standard GD algorithm with the so-called Random Batch Method (RBM).
RBM is a recently developed approach which has been introduced in [13] for the numerical simulation of high-dimensional collective behavior problems.
This method uses small but random batches for particle interactions, lowering the computational cost O(N 2 ) per time step to O(N ), for systems with N particles with binary interactions. Therefore, as our numerical simulations will confirm, embedding RBM into the GD iterative scheme yields to a less expensive algorithm and, consequently, to a more efficient control computation. In what follows, we will call this approach GD-RBM algorithm, to differentiate it from the standard GD one.
3.1. The Gradient Descent approach. Let us now describe in detail the GD approach to minimize the functional (3.3), and discuss its convergence properties.
In order to fully define the iterative scheme (3.5), we need to compute the gradient ∇J(u). Since we are dealing wit a non-linear control problem, we will do this via the so-called Pontryagin maximum principle (see [21] or [29,Chapter 7]).
To this end, let us first rewrite the dynamics (3.1) in a vectorial form as follows In the control literature, (3.6) is usually called the primal system.
Using the notation just introduced, we can then see that J(u) can be rewritten in the form Let us stress that (3.8) is in the standard form to apply the Pontryagin maximum principle. Through this approach, we can obtain the following expression for the gradient of J(u) where D u F indicates the Jacobian of the vector field F , computed with respect to the variable u.
In (3.9), we denoted with p = (p 1 , . . . , p N ) the solution of the adjoint equation associated to (3.1), which is given by (3.10) where D Θ F stands again for the Jacobian of the vector field F , this time computed with respect to the variable Θ.
Taking into account the expression (3.7) of the vector field F , we can then readily check that the iterative scheme (3.5) becomes (3.12) In view of the above computations, the GD algorithm for the minimization of the cost functional J(u) can be explicitly formulated as follows: GD algorithm.
input Θ 0 : initial condition of the primal system (3.6) u 0 : initial guess for the control u k ← 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.6) Solve the the adjoint system (3.12) end for Update the control through the scheme (3.11) end while return u k+1 = 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 2N non-linear differential equations (N for the variables θ i and N for p i ).
If we introduce the time-mesh of N t points 0 = t 0 < t 1 < . . . < t Nt = T, t m = t 0 + m T N t , m = 1, . . . , N t , at each time-step t m this operation has a computational cost of O(N 2 ) and the total computational complexity for the simulation of (3.1) and (3.12) will then be O (N t N 2 ). If N is large, that is, if the amount of oscillators in the network is considerable, this will rapidly becomes very expensive.
3.2. The stochastic approach. In order to reduce the computational burden of GD for the optimization process (3.3), we propose a modification of this algorithm which includes the aforementioned Random Batch Method (RBM) for the numerical simulation of the ODE systems (3.1) and (3.12). This technique, presented in [13] for interacting particle systems, is based on the following simple idea: at each time step t m = m · dt in the mesh we employ to solve the dynamics, we divide randomly the N = nP particles into n small batches with size P < N , denoted by C q , q = 1, . . . , n and then interact only particles within the same batch. This gives the following algorithm for the numerical approximation of (3.1) and (3.12):

RBM algorithm.
for m = 1 to N t = T /dt do Divide randomly {1, . . . , N } into n batches C q , q = 1, . . . , n for q = 1 to n do Update θ i (i ∈ C q ) by solving the ODE Update p i (i ∈ C q ) by solving the ODE

end for end for
Regarding the computational complexity, according to the discussion presented in [13], at each time-step t m the RBM approach has a cost of O(P N ). The total computational complexity for the simulation of (3.1) and (3.12) is thus O(P N t N ), which is always smaller than O (N t N 2 ). Therefore, independently of the value of N , employing RGB to simulate the dynamics in each iteration of GD will yield significant improvements in terms of the computational cost.
3.3. GD convergence analysis. To complete this section, let us briefly comment about the convergence properties of the GD methodology. It is nowadays classically known that the convergence rate of the GD algorithm is determined by the regularity of the objective function. In our case, since J(u) is Lipschitz regular, the algorithm will have linear convergence, that is where, we recall, u denotes the minimum of J(u) and the norm · L 2 (0,T ;R) has been defined in (3.4).
Then, (3.13) implies that, for computing the control u up to some given tolerance ε > 0, i.e. for obtaining e k < ε, the GD algorithm requires k = O(ε −1 ) iterations.
Combining this with the per-iteration costs we gave in Sections 3.1 and 3.2, we can thus obtain the following total computational costs and we can conclude that the GD-RBM approach will be more efficient than the standard GD one to solve our optimization problem. This ins enhanced for large values of N and will be confirmed by our numerical simulations.

Numerical simulations
In this section we present our numerical results for the control of N coupled oscillators described by the Kuramoto model (3.1), following the strategy previously described. In particular, we will show how the RBM approach allows to significantly reduce the computational complexity of the GD algorithm for the calculation of the control u.
The simulations have been performed in Matlab R2018a on a laptop with Intel Core i5 − 7200U CP U @2.50GHz × 4 processor and 7.7 GiB RAM.
The oscillators are chosen such that their natural frequencies are given following the normal probability law with σ = 0.1. This means that the values of ω min and ω max are given respectively by and the coupling gain K which is necessary for synchronization has to satisfy (see In all our numerical experiments, we assumed K to satisfy the above condition. Moreover, we took the penalization parameter β = 10 −7 and we considered the time horizon T = 3s. Finally, as a stopping criterion, we chose with ε = 10 −4 , and where the the notation · 2 stands again for the L 2 (0, T R)-norm defined in (3.4). Let us stress that the stopping criterion (4.2) is consistent with (3.13) and (3.5).
In order to compare the performances of the two methodologies, we have run the simulations for different increasing values of the parameter N , namely N = 6, 20, 40, 50 and 80. When using the GD-RBM approach, we always separated the family of N oscillators into batches of size P = 2 (corresponding to n = 3, 10, 20, 25 and 40, respectively). We collected in Table 1 the computational times we obtained with the two methodologies.  Table 1. Computational times to converge to the tolerance ε = 10 −4 of the GD and GD-RBM algorithms applied with different values of N . The GD-RBM algorithm is always implemented dividing the family of N oscillators into batches of size P = 2.
In Figure 1, we show the evolution of the Kuramoto dynamics with N = 6 in the absence of control, that is, when the control function u is identically one. As we can see, the oscillators are evolving towards a synchronized configuration, which is consistent with our choice of the coupling gain K. Nevertheless, synchronization is not reached in the short time horizon we are providing. On the other hand, in Figure 2, we show the evolution of the same dynamics, this time under the action of the control function u computed through the minimization of J(u). The subplot on the left corresponds to the simulations done 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 = 3s. This means that both algorithms managed to compute an effective control.  We can see how, at the beginning of the time interval we are considering, this control is close to one and it is increasing with a small slope. On the other hand, this growth becomes more pronounced as we get closer to the final time T = 3s.
Notice that, in (3.1), u enters as a multiplicative control which modifies the strength of the coupling K. Hence, according to the profile displayed in Figure 3, the control function u we computed is initially letting the system evolving following its natural dynamics. Then, as the time evolves towards the horizon T = 3s, u enhances the coupling strength K in order to reach the desired synchronized configuration (3.2).
Finally, notice also that the control u is always positive. This is actually not surprising, if one takes into account 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 ( [11]). Moreover, according to the model (3.1), in order to maintain K positive, u has to remain positive as well.
Recall that u is computed minimizing the functional (3.3), in which the second term is a measurement of the level of synchronization in the model. Hence, since negative values of u would lead to desynchronization and to the corresponding increasing of the functional, these values remain automatically excluded during the minimization process.
In Figure 4, we show the convergence of the error when applying both the GD and GD-RBM approach. We can appreciate how, in the case of GD-RBM, this convergence is not monotonic as it is for the GD algorithm. This, however, is not surprising due to the stochastic nature of the RBM methodology. In Figure 5, we repeated the simulations for N = 80 to show that, also in this case, we are able to compute an accurate control for our system. This picture has been obtained by employing GD-RBM, since one of the main interests of this paper is to show the effectiveness of this methodology to control the Kuramoto model (3.1). Finally, in Figure 6, for increasing values of N we compared the computational times (collected in Table 1) for obtaining the control u with the GD and GD-RBM algorithms. We can see how, for low values of N , the two approaches show similar behaviors. Nevertheless, when considering a larger amount of oscillators in our system, the GD-RBM methodology allows to reduce the computational cost for the numerical resolution of (3.1) and (3.12), thus becoming substantially less expensive than the GD one. This is in accordance with our discussion in Section 3.3. Figure 6. Computational time to converge to the tolerance ε = 10 −4 of the GD and GD-RBM algorithms for the control of (3.1) with different values of N . The GD-RBM algorithm is always implemented dividing the family of N oscillators into batches of size P = 2.

Conclusions
This paper deals with the synchronization of coupled oscillators described by the Kuramoto model. In particular, we design a unique scalar control function u(t) capable os steering a N -dimensional network of oscillators to a synchronized configuration in a finite time horizon. This is done following a standard optimal control approach and obtaining the function u(t) via the minimization of a suitable cost functional.
With this approach, we computed a control which acts as a multiplicative positive force enhancing the coupling between 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, commonly employed in the optimal control community, which we have improved with the introduction of a Random Batch Method (RBM) for a more efficient numerical resolution of the Kuramoto dynamics.
Our simulations results show that the inclusion of RBM into the GD algorithm significantly reduces the computational burden for the control computation, thus allowing to deal with high-dimensional oscillators networks in an efficient way.