Optimal synchronization in pulse-coupled oscillator networks using reinforcement learning

Abstract

Spontaneous synchronization is ubiquitous in natural and man-made systems. It underlies emergent behaviors such as neuronal response modulation and is fundamental to the coordination of robot swarms and autonomous vehicle fleets. Due to its simplicity and physical interpretability, pulse-coupled oscillators has emerged as one of the standard models for synchronization. However, existing analytical results for this model assume ideal conditions, including homogeneous oscillator frequencies and negligible coupling delays, as well as strict requirements on the initial phase distribution and the network topology. Using reinforcement learning, we obtain an optimal pulse-interaction mechanism (encoded in phase response function) that optimizes the probability of synchronization even in the presence of nonideal conditions. For small oscillator heterogeneities and propagation delays, we propose a heuristic formula for highly effective phase response functions that can be applied to general networks and unrestricted initial phase distributions. This allows us to bypass the need to relearn the phase response function for every new network.

Introduction

Recent research has turned to using the pulse-coupled oscillator (PCO) model (1–10) that was initially proposed to describe biological neuronal networks (11–22) and cardiac pacemakers (23–26), but has subsequently found applications in many other systems, including artificial neural networks (27), social self-organization (28), and clock coordination of wireless sensor networks (29, 30). The PCO model explicitly incorporates the hybrid nature of network dynamics (in contrast to other models such as the Kuramoto oscillators (31)) and promises great potential for synchronizing engineered networks (29). For example, the PCO-based strategies have been found to be highly successful in achieving motion coordination in robot swarms (32–34). Since only simple and content-free pulses are sent between agents, PCO-based synchronization is naturally resilient to message corruption in communications and incurs little communication overhead, leading to improved network robustness and reduced communication latency (35, 36). These advantages make PCO-based synchronization particularly appealing for the coordination in engineering networks such as robot swarms and vehicle fleets that are subject to stringent reliability and real-time constraints.

The synchronizability of PCO networks depends crucially on the phase response function (PRF), which characterizes how an oscillator adjusts its phase when a pulse is received from a neighboring oscillator (37). The amount of adjustment depends on the current state of the receiving oscillator. Many analytical results on PRF have been reported under ideal conditions. For example, Klinglmayr et al. (38) and Lyu (39) investigated the synchronization of PCOs with stochastic PRFs. Lyu (40) further analyzed the synchronization time of PCOs on tree networks. Wang and Doyle (41) proposed a PRF that can maximize the speed of synchronization of PCOs when the initial oscillator phases are distributed within a half-cycle. However, these results are often obtained under ideal conditions, including zero time delays and identical oscillator frequencies. In fact, nonideal factors such as propagation delays and heterogeneous oscillator frequencies render the analytical design of PRF extremely difficult, if at all possible (9, 42–49, 50–52).

In this article, we propose a reinforcement learning (RL) approach to determine a highly effective PRF under both ideal and nonideal conditions. The interaction strategy found by our RL approach improves synchronization probability compared to previously proposed PRFs in (41, 42, 43, 47). Moreover, our results provide insights on a general design principle for PRF that can be adapted to general network topologies. Finally, the flexibility of our RL framework allows oscillators to adapt to changing network structures and environmental noise, ensuring robust synchronization under real-world conditions.

Results

Pulse-coupled oscillators

Let us consider a network of N PCOs, where ${\varphi }_i\in {\mathbb{S}}^1=[0,2\pi )$ is the phase of oscillator $i\in \mathcal{V}=\{1,2,\dots ,N\}$⁠. Each oscillator evolves its phase at a frequency ${\omega }_i$⁠. When ${\varphi }_i$ reaches the threshold value $2\pi$⁠, oscillator i fires a pulse and resets its phase to zero. Neighboring oscillators receive this pulse after some (random) time delay ${\tau }_{ij}$⁠. This delay in receiving a pulse is primarily due to the finite propagation speed of the pulse, but it may also include the time required for a node to process the incoming pulse.

An oscillator responds to a received pulse by changing its phase ${\varphi }_i$ by

where ${\varphi }_i^{-}$ and ${\varphi }_i^{+}$ represent the phase of oscillator i immediately before and after receiving a pulse, respectively. The function $F(\varphi )$⁠, which determines the amount that an oscillator will adjust its phase as a function of its phase value upon which the pulse is received, is called the phase response function (PRF). The jump in the value of the phase ${\psi }_i$ is determined not only by the PRF, but also by the coupling strength $l\in (0,1)$⁠, which is introduced to facilitate the analysis and design of PCO-based synchronization in engineered systems (41, 42, 49, 53). It is worth noting that PRF is related to the phase transition curve (PTC) in (43) as $\mathrm{PTC}={\varphi }^{-}+l\mathrm{PRF}$⁠.

To quantify the degree of synchronization of an oscillator network, we define the containing arc $\mathrm{\Lambda }(\varphi )$ as the smallest interval in ${\mathbb{S}}^1$ that contains all oscillator phases in the network:

Following (54), we define an arc as a connected subset of the one dimension torus ${\mathbb{S}}^1$⁠. Thus, ${\upsilon }_i(\varphi )=\underset{j\ne i}{min}\{({\varphi }_j-{\varphi }_i)\text{mod}2\pi \}$ in the preceding Eq. (2) denotes the length of the arc along ${\mathbb{S}}^1$ to the first oscillator ahead in the phase of oscillator i. Note that $\sum _{i\in \mathcal{V}}{\upsilon }_i(\varphi )=2\pi$ always holds. Hence, $\mathrm{\Lambda }(\varphi )=2\pi -\underset{i\in \mathcal{V}}{max}{\upsilon }_i(\varphi )$ is the smallest arc containing all oscillators and can be used to quantify the degree of synchronization.

We use RL to determine an optimal PRF $F(\varphi )$ that maximizes the probability of synchronization for a given number of oscillators:

where G is an underlying graph with edges drawn randomly according to some distribution, $\varphi (0)$ is an initial phase configuration drawn uniformly at random, and $\varphi (t)$ is the phase trajectories of all oscillators under a given PRF and initial phase distribution.

Reinforcement learning

A number of works have recently been reported that use learning methods to investigate the dynamics of coupled oscillators (see, e.g. 55–58). However, all of these results consider smooth-interaction oscillators like Kuramoto oscillators. Recently the authors of (59) propose to use learning methods to predict if an oscillator network can synchronize or not, and they consider both smooth and pulse interactions. The approach in (59) considers interaction mechanisms that are given and predetermined. In this work, we leverage the exploration and adaptation properties of RL to optimize the interaction mechanism of PCOs. More specifically, we use RL to determine an optimal PRF $F(\varphi )$ that maximizes the probability to synchronize under both ideal and nonideal conditions. It is worth noting that formal analysis of PCO networks under general initial phase distributions and practical nonideal conditions, such as coupling delay and frequency heterogeneity, remains out of reach with current analytic techniques (47, 60). With RL strategies, we can model the nonideal factors, and let the oscillators evolve naturally in the network and gradually optimize their response to maximize synchronization probability.

A schematic of the approach is given in Fig. 1. Specifically, for each oscillator $i\in \mathcal{V}$⁠, the state is its phase value $s={\varphi }_i$⁠. When a pulse is received from a neighboring oscillator, oscillator i changes its phase value by an action $a=F({\varphi }_i)$⁠. The value of each possible state-action pair is described as a matrix $Q(s,a)$⁠. Based on the current state-action values, oscillator i chooses a policy ${\pi }_i$ under its environment, which consists of its neighboring oscillators along with the network dynamics. Each oscillator receives a set of reward values, which are used to update the state-action values. The RL process repeats until an optimal policy ${\pi }_i^{*}$ is obtained. By taking the average of ${\pi }_i^{*}$⁠, we find a best piece-wise linear fit as the optimal PRF. The detailed RL protocol, such as the designs of reward values and Q-value update are described in the “Materials and Methods” section.

Optimal phase response functions under ideal conditions

We first consider the ideal case where oscillators have identical frequency ${\omega }_i=2\pi$ and pulses are received instantaneously with no delay. We begin with a PCO network of $N=2$ oscillators, where the containing arc is always within a half of a cycle. Fig. 2 shows the average learned policy over 10 experiments, which closely approximates the proposed PRF in (41) obtained analytically under the assumption that the containing arc is less than half of a cycle:

We next consider $N=6$ oscillators with all-to-all coupling (more results for $N=3$⁠, 4, 5 and ER random graphs see Fig. S1–S6 in the supplementary materials). Fig. 3 shows the optimal average learned policy over 18 experiments, which is different from the analytical result in (41), primarily due to the unrestricted phase distributions. During the training phase, the randomness in the policy selection does not guarantee that the oscillators stay within a half-cycle,even if the phase values were initialized within a half-cycle. Thus, the assumptions for the analytical result in (41) do not hold in this case.

Design principle of phase response functions

The average learned policies from both Figs. 2 and 3 can be modeled using a simple form

where $c_1<2\pi$ and $c_2<2\pi$ are positive constants that will be determined later. The best-fits to the learned policies using Eq. (5) are shown in Figs. 2 and 3.

This PRF model offers important insight into the design principles of the highly effective phase response policy. When its phase value is close to the start or end of a cycle, an oscillator learns to take the maximum phase adjustment toward the threshold value, which has been shown analytically to decrease the synchronization time in (41). But when the phase value is near the middle of the cycle, the oscillator adjusts its phase proportionally to the distance to the end of the cycle, similar to the strategy used in (43), which has been shown to almost always lead to synchronization. This combination of two different strategies gives a phase response function that achieves synchronization efficiently.

We now use Erdős–Rényi–Gilbert (ERG) networks with parameters chosen at random to find the best parameters for Eq. (5). Fig. 4 shows the best-fit values for the parameters $c_1$ and $c_2$⁠, which fit well to exponential functions of the oscillator indegree ${\delta }^{-}$⁠:

where $b_{i,1}$ and $b_{i,2}$ are constants. This formula reduces the optimization of a function to the much simpler problem of optimizing two parameters. Equation (5) combined with Eq. (6) provides a powerful heuristic formula that can be used to predict the optimal PRF for oscillators in general networks without repeating the RL process.

Optimal phase response functions under nonideal conditions

When we consider experiments using various nonidentical oscillator frequencies (with differences within $10\mathrm{\%}$ of the nominal frequency) and nonzero propagation delays (with delays within $10\mathrm{\%}$ of an oscillation cycle), we find that learned PRFs are very similar to the ones for ideal environments (see Fig. S7–S14 in the supplementary materials for more details). Thus, we conclude that the same PRFs are optimal in achieving synchronization when the network is subject to moderate nonideal conditions.

Next, we compare our learned PRF based on Eq. (5) to previously proposed PRFs in (41, 42, 43, 47) under nonideal conditions. Since the synchronization strategies by Mirollo and Strogatz (43), Nishimura and Friedman (47), and Klinglmayr (42) do not include an explicit coupling strength parameter, we modify these algorithms to scale the change in phase adjustment by the coupling strength l, as in Eq. (1). These algorithms are recovered under $l=1.0$⁠.

For our test, we consider an ER random graph with $N=10$ oscillators and edge probability $p=0.3$⁠, as illustrated in Fig. 5. We set the oscillator frequency ${\omega }_i$ uniformly distributed in the range $[1.9\pi ,2.1\pi ]$ and propagation delays ${\tau }_{ij}$ uniformly distributed in the range $[0.01,0.08]$ cycles. We randomly assign initial oscillator phase values in ${\mathbb{S}}^1$⁠.

Fig. 6 shows the average value of the containing arc at steady-state after 2,000 runs over a broad range of coupling strength values. Our learned PRF, Wang and Doyle’s delay-advance PRF in (41), and Nishimura and Friedman’s “strong type II” PRF in (47) are able to synchronize the network more closely than the other synchronization strategies due to the similarity among these PRFs. The addition of the coupling strength parameter allows for additional tuning to achieve a greater level of synchronization, especially in nonideal environments.

Moreover, if we look at how often the network is able to synchronize, we see that our learned PRF is able to synchronize more often at low coupling strengths than any of the other synchronization strategies, as shown in Fig. 7. Similar results are obtained when we vary the network size and topology, the amount of oscillator frequency heterogeneity, and the amount of coupling delay (see Figs. S15–S25 in the supplementary materials for more details).

Discussion

The proposed reinforcement learning framework provides a simple and versatile method for optimizing synchronization in pulse-coupled oscillator networks to maximize the degree and resilience of synchronization. Given that maintaining synchronization in the presence of message corruption and communication delays is crucial for numerous systems and processes, the results are expected to have broad applications in biological and engineered systems. For example, biological systems such as cardiac pacemakers and neuron networks can be effectively modeled by PCOs (61). In addition, the proposed method are well positioned to synchronize clocks in wireless sensor networks (41, 42, 53) and coordinate motions in robot networks (32, 34). Furthermore, to the best of our knowledge, this paper is the first to use a distributed reinforcement learning approach to optimize synchronization under noncontinuous pulse-based interactions, which is different from continuous-time smooth interactions in the Kuramoto model (62). It has direct ramifications for the deployment of reinforcement learning in general multiagent systems to optimize dynamical processes in general networks.

Future works include minimizing the synchronization time for PCOs and incorporating mechanisms to allow oscillators to adapt their intrinsic frequencies to deal with large frequency heterogeneity. Further improvements in the design of reward values, along with expanded action space (e.g. allow adjustments of both oscillator phases and frequencies), may lead to better control of both synchronization and desynchronization in networks of sensors and robots with discrete-time.

Materials and methods

Reinforcement learning protocol

State, action, and Q-value

Since the oscillator’s state and action values evolve on continuous intervals, parameterization and approximation decisions must be made to implement RL. We parameterize the continuous state interval into $P+1$ evenly spaced parameters, $s_0,s_1,\dots ,s_P$⁠. The state parameter $s_p$ corresponds to the phase value $\frac{2\pi p}{P}$ in ${\mathbb{S}}^1$⁠. Additionally, we discretize the actions with $A+1$ values, $a_0,a_1,\dots ,a_A$⁠, such that the actions are limited to phase changes that keep the oscillator within its current cycle. The possible actions that can be taken by an oscillator for $s_p$ can be expressed as $a_k=-s_p+\frac{2\pi k}{A}$ for $k=0,1,\dots ,A$⁠. We represent the value of each state-action pair with a $(P+1)\times (A+1)$ matrix $Q(s,a)$⁠. Each element of $Q(s,a)$ estimates the amount of expected reward by taking the action a at state s.

We implement episodic RL using an on-policy temporal-difference RL technique (see, e.g. 63). Off-policy RL techniques, such as Q-learning, tend to perform worse when there is a need to approximate continuous states and action spaces based on (64). A policy π for our MDP consists of a set of actions, one action for each state parameter $s_p$⁠, such that $\pi (s_p)=a_p$⁠. This policy represents a straight-line approximation of the PRF for an oscillator $F(\varphi )$⁠. For on-policy learning, we choose a policy before each episode. To avoid confusion with phase $2\pi$⁠, we denote policies with subscripts, such that the policy used for episode t is ${\pi }_t$⁠. The choice for a policy is based on the current state-action value estimates $Q(s,a)$⁠. We use a soft-max, or Boltzmann, distribution to choose an action. Initially, the values for all state-action pairs are set to zero. Thus, the initial policy is equally likely to choose any action.

Reward design

Our goal is to have the oscillators synchronize their phases. Due to the dynamics of PCO networks, the choice of reward and how to update the values of state-action pairs are both critical. Since the length of the containing arc measures how well the network is synchronized, we reward actions that decrease the length of the containing arc and penalize actions that increase that length. Thus, the reward will include the decrease in the length of the containing arc, $\mathrm{\Delta }\mathrm{\Lambda }(\varphi )=\mathrm{\Lambda }(\varphi {)}^{-}-\mathrm{\Lambda }(\varphi {)}^{+}$⁠, where $\mathrm{\Lambda }(\varphi {)}^{-}$ and $\mathrm{\Lambda }(\varphi {)}^{+}$ are the lengths of the containing arc before and after the action was taken, respectively.

Individual oscillators in a PCO network do not know the true state of the network. However, each oscillator can approximately estimate the state of the network by keeping track of the phase differences between itself and other oscillators when a pulse is received. Therefore, we can use the oscillator’s estimated state to approximate the change in the containing arc and, thus, determine the reward value.

Multiple actions can result in the same decrease in the containing arc. To encourage efficient synchronization, we penalize the oscillator by $f(a_k)=\frac{a_k^2}{2\pi }$ based on the magnitude of the action $a_k$⁠. Therefore, when an oscillator receives its kth pulse and takes an action $a_k$ that decreases the (estimated) containing arc by an amount $\mathrm{\Delta }{\mathrm{\Lambda }}_k(\varphi )$⁠, the total reward is given by

where $w_{\mathrm{\Lambda }}$ and $w_a$ are positive weights.

During a training episode t, we let the network evolve for a fixed amount of time following a given policy ${\pi }_t$⁠. When an oscillator receives a pulse, it uses ${\pi }_t$ to determine its action, i.e. phase adjustment, based on a given coupling strength l and its current phase ${\varphi }_i$⁠.

Let us denote the two closest state parameters to ${\varphi }_i$ for the kth received pulse as $s_{L,k}$ and $s_{H,k}$⁠, respectively, where $s_{L,k}\le {\varphi }_i\le s_{H,k}$ holds. The corresponding actions from policy ${\pi }_t$ for those state parameters are denoted as $a_{L,k}$ and $a_{H,k}$⁠, respectively. To determine the action ${\psi }_i$ that the oscillator takes, we weight the actions of the two nearest state parameters based on the proximity of the current phase to those state variables. That is, the phase adjustment ${\psi }_i$ in Eq. (1) is

with ${\rho }_{L,k}=\frac{s_{H,k}-{\varphi }_i}{s_{H,k}-s_{L,k}}$ and ${\rho }_{H,k}=\frac{{\varphi }_i-s_{L,k}}{s_{H,k}-s_{L,k}}$⁠.

Before an oscillator adjusts its phase, it records its current state as $S_k={\varphi }_i$⁠, where k is an index for the number of received pulses during an episode. The action taken, based on the policy ${\pi }_t$⁠, is recorded as $A_k={\psi }_i$⁠. The oscillator calculates its reward $R_k$ based on Eq. (7) and then evolves freely until another pulse is received.

Q-value update

Once an episode for the network is complete, we use the resulting state-action-reward sequences to perform a batch update of the state-actions value matrix $Q(s,a)$⁠. The update is based on the Sarsa algorithm in (64), where the reward $R_k$ will apply to the state-action values for the two nearest state parameters to $S_k$ and their corresponding actions from the episode policy ${\pi }_t$⁠. The update for $s_{L,k}$ is

and the update for $s_{H,k}$ is

where α is the learning rate and γ is the discount rate. Here, $Q_{E,k+1}$ denotes the average estimated value of the next state-action pair (⁠$S_{k+1},A_{k+1}$⁠) and is calculated as

This update is performed for every state-action pair except for the final pair, which is the terminal state.

After all updates have been completed for each oscillator’s state-action-reward sequence, an episode of training is complete. With the updated state-action value matrix $Q(s,a)$⁠, a new policy is chosen for the next episode, and the process is repeated. After all episodes of training are complete, we use $Q(s,a)$ to determine the optimal policy ${\pi }^{*}$⁠, where ${\pi }^{*}(s_p)={\text{arg max}}_{a_i}Q(s_p,a_i)$ for each state parameter $s_p$⁠. We note that different oscillators in a network can have different optimal policies ${\pi }^{*}$⁠.

In the implementation, for each case, we let the network evolve for 15 cycles for each episode with initial phases randomly selected in $[0,2\pi )$⁠, and use a coupling strength $l=1.0$⁠. We parameterize the state with 101 evenly spaced parameters and discretize the policy actions into 201 evenly spaced values. For each experiment, we simulate a network for 100,000 episodes. We use $w_{\mathrm{\Lambda }}=\frac{N}{N-1}$ and $w_a=\frac{1}{l}$ for Eq. (7).

Acknowledgments

We thank Steven Strogatz for helpful discussions.

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

The work was supported in part by the National Science Foundation under Grant 2012075.

Authors’ Contribution

T.A. and Y. W. conceived the project. Z.C, T.A., and Y.W. performed the research. Z.C, T.A., Y.Z. and Y.W. discussed the results. Z.C, T.A., Y.Z. and Y.W. wrote the paper.

Data Availability

All data is included in the manuscript and/or supporting information.

References

Author notes