Robust Almost Global Splay State Stabilization of Pulse Coupled Oscillators

TL;DR

This paper presents a new method for stabilizing the evenly distributed phase configuration of pulse coupled oscillators, ensuring robustness and almost global stability through a novel Lyapunov approach.

Contribution

It introduces a novel phase response function design and Lyapunov function for almost global stabilization of splay states in pulse coupled oscillators.

Findings

Achieves almost global asymptotic stability of splay states

Demonstrates robustness to frequency perturbations

Validates effectiveness through an example

Abstract

This technical note deals with the problem of asymptotically stabilizing the splay state configuration of a network of identical pulse coupled oscillators through the design of the their phase response function. The network of pulse coupled oscillators is modeled as a hybrid system. The design of the phase response function is performed to achieve almost global asymptotic stability of a set wherein oscillators' phases are evenly distributed on the unit circle. To establish such a result, a novel Lyapunov function is proposed. Robustness with respect to frequency perturbation is assessed. Finally, the effectiveness of the proposed methodology is shown in an example.

Equations

$$\mathcal{H}\left\{\begin{array}[]{ccll}\dot{x}&=&f(x)&\quad x\in C\\ x^{+}&\in&G(x)&\quad x\in D\end{array}\right.$$

$$\begin{array}[]{lr}C\coloneqq[0,2\pi]^{n},&f(x)\coloneqq{\bf 1}_{n}\omega\qquad\forall x\in C\end{array}$$

$$g_{i}(x)\coloneqq\left\{\begin{array}[]{ll}0&\mbox{if}\,x_{i}=2\pi,\forall j\in\mathcal{V}\setminus\{i\}\,x_{j}\neq 2\pi\\ x_{i}+Q(x_{i})&\mbox{if}\,x_{i}\neq 2\pi,\exists j\in\mathcal{V}\setminus\{i\}\,x_{j}=2\pi\\ \{0,x_{i}+Q(x_{i})\}&\mbox{if}\,x_{i}=2\pi,\exists j\in\mathcal{V}\setminus\{i\}\,x_{j}=2\pi\end{array}\right.$$

$$\mathcal{A}=\left\{x\in C:y=\mathrm{sort}(x),d(y_i,y_{{\mathrm{next}}_{\mathcal{V}}(i)})=\frac{2\pi }{n}\forall i\in \mathcal{V}\right\}$$

$$\mathcal{X}=\{x\in C:\exists i,j\in \mathcal{V}\text{mbox}s.t.i\ne j,\text{mbox}andd(x_i,x_j)=0\}$$

$$\underset{i\in \mathcal{V}}{\max }d(x_i,x_{{\mathrm{next}}_{\mathcal{V}}(i)})=\frac{2\pi }{n}$$

$${\mathcal{Q}}_1:=\left\{(w_1,w_2)\in {\mathbb{R}}^2:w_1\in \left[0,\frac{2\pi (n-1)}{n}\right]\wedge w_2=0\right\}$$

$$\begin{array}{rl}{\mathcal{Q}}_2& :=\{(w_1,w_2)\in {\mathbb{R}}^2:\\ & w_1\in \left(\frac{2\pi (n-1)}{n},2\pi \right]\wedge \frac{2\pi (n-1)}{n}-w_1<w_2<0\}\end{array}$$

$$\mathcal{H}_{d}\left\{\begin{array}[]{ll}\begin{array}[]{rcl}\dot{x}&=&\omega{\bf 1}_{n}+d\\ \end{array}&x\in C\\ \begin{array}[]{rcl}x^{+}&\in&G(x)\\ \end{array}&x\in D\end{array}\right.$$

$$Q(z)=\{\begin{array}{ll}0& \text{mbox}ifz\in [0,\frac{4\pi }{3}]\\ -\frac{7}{10}(z-\frac{4\pi }{3})& \text{mbox}ifz\in [\frac{4\pi }{3},2\pi ]\end{array}$$

Full text

I Introduction

Pulse coupled oscillators are oscillators connected together to form a network of systems interacting through the exchange of impulsive messages. In particular, each time an oscillator completes a cycle it sends a pulse (fires) to the other oscillators. Such oscillators, in response to the received pulse, reset their phase according to a phase response function; see [2, 8]. Despite their simplicity, pulse coupled oscillators may manifest complex collective behaviors. For such a reason, their adoption as modeling and design tools in complex biological and engineered systems has provided interesting results in different fields; see, e.g., [8, 11, 12] just to cite a few. Concerning engineered systems, the most investigated behaviors in pulse coupled oscillators are synchronization and desynchronization. When a network of pulse coupled oscillators evolve in synchronization, oscillators evolve with the same phase each other. When a network of pulse coupled oscillators reaches desynchronization, oscillators have different phases and the time between consecutive firings is constant [11]. Desynchronization has relevant applications in engineered systems as it provides a simple, robust, and decentralized way to generate round-robin schedulers; [11]. In particular, as pointed out in [11], a network of pulsed coupled oscillators can manifest two possible kinds of desynchrony: splay state (or strict desynchronization) and weak desynchronization. In the splay state, oscillators’ phases are evenly distributed around the unit circle ([18]) implying that the time between consecutive firing instants is constant. Instead, in the weak desynchronization, there is no equal spacing between oscillators’ phases and equal separation of firing times is achieved by making oscillators decrease their phase at each firing time. This kind of desynchronized behavior has been studied in [13] via hybrid systems tools. It is worthwhile to observe that although the occurrence of both behaviors ensures constant firing intervals, since in the splay state configuration no interaction between oscillators is needed, such a configuration reveals to be more robust to pulse lossess and therefore preferable in some applications; see [11].

In this technical note, we address the problem of stabilizing the splay state configuration of a network of pulse coupled oscillators connected in an all-to-all topology. Inspired by [9, 10, 13], we adopt a hybrid systems approach to tackle the considered problem. In particular, we model the network of pulse coupled oscillators as a hybrid system $\mathcal{H}$ in the framework of [7]. Then, building on this model, similarly as in [13], we recast the splay state stabilization problem as a set stabilization problem for $\mathcal{H}$. As a second step, by relying on the notion of shortest containing arc ([4, 10]), we show that the splay state configuration coincides with the (unique) configuration for which the length of the shortest containing arc is maximized. Building on such a property, we define a Lyapunov function related to the length of the shortest containing arc. Finally, thanks to the invariance principle for hybrid systems [16], inspired by [11], we provide sufficient conditions on the phase response function to ensure almost global asymptotic stability of the splay state configuration. Furthermore, by relying on some regularity properties satisfied by the data of $\mathcal{H}$, we show that the splay state configuration is structurally robust to frequency perturbations, making the results presented in this paper appealing in practice. Finally, the theoretical findings are illustrated via a numerical example.

Contribution

The contribution of this paper are as follows. Firstly, we formalize the splay state stabilization problem for pulse coupled oscillators in the hybrid systems framework in [7], which allows to obtain formal robustness guarantees about the perturbed behavior of the network. In this sense, this work can be seen as complementary to [13], where weak desynchronization is considered. Secondly, we propose sufficient conditions on the phase response function ensuring asymptotic stability of the splay state configuration. Thirdly, for the first time, we make use of the notion of shortest containing arc in a desynchronization setting to build a Lyapunov function. Observe that, as opposed to synchronization where the relationship between the length of the shortest containing arc and the synchronous configuration is obvious, in the case of splay state, such a relationship is nontrivial and needs to be worked out.

A preliminary version of this work was presented in [5]. On the one hand, this paper aims at completing the preliminary results presented in [5] by providing further technical clarifications, refinements, and proofs. On the other hand, a thorough comparison with existing results as well as the analysis in the presence of frequency perturbations are presented.

The remainder of the technical note is organized as follows. Sections I.A concerns preliminary notions/results. Section I.B gives some basic notions on hybrid systems. Section II presents the system under consideration, the hybrid modeling and the problem we solve. Section III is dedicated to the main results. Section IV is devoted to a numerical example. Numerical simulations of hybrid systems are performed in Matlab via the Hybrid Equations Toolbox [15].

Notation: The set $\mathbb{N}$ is the set of strictly positive integers, $\mathbb{N}_{0}=\mathbb{N}\cup\{0\}$, $\mathbb{R}_{\geq 0}$ represents the set of non-negative real scalars, and $\mathbb{C}$ represents the set of complex numbers. We denote $\iota=\sqrt{-1}$ as the imaginary unit and $\mathbb{S}^{1}=\{z\in\mathbb{C}\colon|z|=1\}$ the unit circle, where $|z|$ is the modulus of $z$. Given a set $S$, $\mathcal{C}(S)$ denotes the cardinality of $S$. For a vector $x\in\mathbb{R}^{n}$, $x_{i}$ denotes its $i$-th entry and $|x|$ the Euclidean norm. The symbol ${\bf 1}_{n}$ denotes the vector of $\mathbb{R}^{n}$ whose entries are equal to one. Given two vectors $x,y$, we denote $(x,y)=[x^{\prime}\,\,y^{\prime}]^{\prime}$. The set $\mathbb{B}$ is the closed unit ball, of suitable dimensions, in the Euclidean space. Given a vector $x\in\mathbb{R}^{n}$ and a closed set $\mathcal{A}$, the distance of $x$ to $\mathcal{A}$ is defined as $|x|_{\mathcal{A}}=\inf_{y\in{\mathcal{A}}}|x-y|$. Given $\mathcal{I}\subset[0,2\pi]$, we denote $\Psi(\mathcal{I})=\{z\in\mathbb{S}^{1}\colon z=e^{i\omega}\quad\omega\in\mathcal{I}\}$. Given a set $S$, $m(S)$ stands for the Lebesgue measure of $S$. Let $n\in\mathbb{N}$, we denote $\mathcal{V}=\{1,2,\dots,n\}$. Given $i\in\mathcal{V}$, we denote $\operatorname{next}_{\mathcal{V}}(i)=i+1$ if $i<n$ while $\operatorname{next}_{\mathcal{V}}(i)=1$ if $i=n$. Given a vector $x\in\mathbb{R}^{n}$, we denote $\operatorname{sort}(x)$ the vector whose entries correspond to the entries of $x$ sorted from the smallest to the largest value. Given a function $f\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, we denote $\operatorname{rge}f$ the range of $f$. Let $F\colon\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{m}$ be a set-valued mapping, we denote $\operatorname{gph}F=\{(x,y)\in\mathbb{R}^{n}\times\mathbb{R}^{m}\colon y\in F(x)\}$. Given a function $f\colon\operatorname{dom}f\rightarrow\mathbb{R}$ and a real scalar $\mu$, $f^{-1}=\{x\in\operatorname{dom}f\colon f(x)=\mu\}$ is the $\mu$-level set of $f$.

I-A Preliminaries

Definition 1**.**

Given $\theta_{1},\theta_{2}\in\mathbb{R}$, $\operatorname{d}(\theta_{1},\theta_{2})$ represents the geodesic distance between the point $x_{1}=e^{\iota\theta_{1}}$ and $x_{2}=e^{\iota\theta_{2}}$ on the unit circle $\mathbb{S}^{1}$.

In the sequel, with a slight abuse of notation, given $\theta_{1},\theta_{2}\in\mathbb{R}$, we refer to $\operatorname{d}(\theta_{1},\theta_{2})$ as the geodesic distance between $\theta_{1}$ and $\theta_{2}$.

Definition 2**.**

A subset $\alpha\subset\mathbb{S}^{1}$ is said to be an arc if it is closed and connected. Given $\theta\in\mathbb{R}^{n}$ and an arc $\alpha\subset\mathbb{S}^{1}$. We say that $\alpha$ contains $\theta$ if for each $i=1,2,\dots,n$, $x_{i}=e^{\iota\theta_{i}}\in\alpha$. In particular, given $\theta\in\mathbb{R}^{n}$, we denote $\Gamma(\theta)\subset\mathbb{S}^{1}$ the set of all the arcs containing $\theta$. Given an arc $\alpha\subset\mathbb{S}^{1}$, we denote $\ell(\alpha)\in[0,2\pi)$ the length of the arc $\alpha$.

Definition 3**.**

Given $\theta\in\mathbb{R}^{n}$, we denote $\gamma(\theta)$ the length of the shortest arc containing $\theta$. In particular111It can be readily shown that for any $\theta\in\mathbb{R}^{n}$, $\mathcal{C}(\Gamma(\theta))=n$. Therefore, $\gamma$ is well defined., $\gamma(\theta)=\min_{\omega\in\Gamma(\theta)}\ell(\omega)$.

Now, consider the following result which will be useful in the sequel.

Lemma 1**.**

For each $\theta\in\mathbb{R}^{n}$, one has $\gamma(\theta)\leq 2\pi-\frac{2\pi}{n}$.

Proof.

Without loss of generality, assume by contradiction that there exists $\hat{\theta}\in[0,2\pi]^{n}$ with $0=\hat{\theta}_{1}\leq\hat{\theta}_{2}\leq\dots\leq\hat{\theta}_{n}$, such that $\gamma(\hat{\theta})>2\pi-\frac{2\pi}{n}$. Since $2\pi-\frac{2\pi}{n}\geq\pi$, as argued in [9], one has that $\gamma(\hat{\theta})=2\pi-\underset{i\in\mathcal{V}}{\max}\operatorname{d}(\hat{\theta}_{i},\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)})$, hence for each $i\in\mathcal{V}$, $\operatorname{d}(\hat{\theta}_{i},\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)})<\frac{2\pi}{n}$. On the other hand, since by assumption $\gamma(\hat{\theta})>2\pi-\frac{2\pi}{n}\geq\pi$, it can be easily noticed that necessarily for each $i\in\mathcal{V}\setminus\{n\}$, $\operatorname{d}(\hat{\theta}_{i},\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)})=|\hat{\theta}_{i}-\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)}|=\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)}-\hat{\theta}_{i}$ and $\operatorname{d}(\hat{\theta}_{n},\hat{\theta}_{1})=2\pi-\hat{\theta}_{n}$. Thus, $\hat{\theta}_{n}=\sum_{i\in\mathcal{V}\setminus\{n\}}(\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)}-\hat{\theta}_{i})<2\pi\frac{n-1}{n}$ but this implies that $\operatorname{d}(\hat{\theta}_{n},\hat{\theta}_{1})=2\pi-\hat{\theta}_{n}>\frac{2\pi}{n}$ contradicting the fact that for each $i\in\mathcal{V}$, $\operatorname{d}(\hat{\theta}_{i},\hat{\theta}_{\operatorname{next}_{\mathcal{V}}(i)})<\frac{2\pi}{n}$. ∎

I-B Hybrid systems

We consider hybrid systems with state $x\in\mathbb{R}^{n}$ of the form

[TABLE]

that we represent by the shorthand notation $\mathcal{H}=(C,f,D,G)$. A solution to $\mathcal{H}$ is any hybrid arc defined over a subset of $\mathbb{R}_{\geq 0}\times\mathbb{N}_{0}$ that satisfies the dynamics of $\mathcal{H}$. A solution to a hybrid system is said to be complete if its domain is unbounded and maximal if it is not the truncation of another solution. Given a set $S$, we denote $\mathcal{S}_{\mathcal{H}}(S)$ the set of all maximal solutions $\phi$ to $\mathcal{H}$ with $\phi(0,0)\in S$. Given $\mathcal{H}=(C,f,D,G)$, we say that $\mathcal{H}$ satisfies the hybrid basic conditions ([7]) if: $C$ and $D$ are closed in $\mathbb{R}^{n}$; $f\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is continuous on $C$, $G\colon\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is an outer semicontinuous set-valued map 222A set-valued map $M\colon\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is outer semicontinuous if its graph is closed; see [14]. nonempty and locally bounded on $D$.

Definition 4** ($(\tau,\varepsilon)$-closeness of hybrid arcs [7]).**

Given $\tau\geq 0$ and $\varepsilon>0$, two hybrid arcs $\xi_{1}$ and $\xi_{2}$ are $(\tau,\varepsilon)$-close if

• (a)

for all $(t,j)\in\operatorname{dom}\xi_{1}$ with $t+j\leq\tau$ there exists $s$ such that $(s,j)\in\operatorname{dom}\xi_{2}$, $|t-s|<\varepsilon$, and $|\xi_{1}(t,j)-\xi_{2}(s,j)|\!<\!\varepsilon$;

• (b)

for all $(t,j)\in\operatorname{dom}\xi_{2}$ with $t+j\leq\tau$ there exists $s$ such that $(s,j)\in\operatorname{dom}\xi_{1}$, $|t-s|<\varepsilon$, and $|\xi_{2}(t,j)-\xi_{1}(s,j)|\!<\!\varepsilon$.

II Problem statement

II-A System description

In this paper, we analyze an all-to-all network of $n\in\mathbb{N}\setminus\{1\}$ identical pulse-coupled oscillators (PCO). Each oscillator is characterized by a phase variable $x_{i}\in[0,2\pi]$ for each $i\in\mathcal{V}\coloneqq\{1,2,\dots,n\}$, that evolves continuously from [math] to $2\pi$ according to integrate-and-fire dynamics, i.e., $\dot{x}_{i}=\omega$, where $\omega\in\mathbb{R}_{>0}$ is the natural frequency of the oscillators; when, $x_{i}$ reaches $2\pi$, the oscillator $i$ fires, i.e., emits a pulse and resets its phase to zero. Such a pulse is instantaneously broadcast through the network to the other oscillators $x_{j}$ with $j\in\mathcal{V}\setminus\{i\}$ that in turn reset their phases according to a phase response function $Q$ [2], i.e., $x_{j}^{+}=x_{j}+Q(x_{j})$, where $x^{+}_{j}$ denotes the value of $x_{j}$ right after a reset.

The problem studied in this technical note consists of designing the function $Q$ in a way such that the configuration in which the $n$ phases are evenly distributed on the unit circle (splay state [18]) is asymptotically stable. In the sequel, we refer to this problem as splay state stabilization problem.

Remark 1**.**

Assuming an all-to-all communiction between oscillators allows to provide a simple yet robust solution to the considered problem, which is therefore suitable for practical implementations. This assumption is common in the literature on desynchronization of PCOs; see, e.g., [11, 13]. Indeed, although this assumption simplifies the analysis of the resulting network, it is often used in engineered systems, like wireless network [11] and in the study of some biological systems such as cardiac pacemaker cells [12] and populations of spiking neurons [6]. Further insights about the complications encountered in the case of more general topologies are illustrated throughout the paper.

II-B Hybrid modeling and dynamical properties

Due to the hybrid behavior of a network of PCOs, as in [9, 13], in this paper we model such a network as a hybrid system $\mathcal{H}$ with state $x=(x_{1},x_{2},\dots,x_{n})\in\mathbb{R}^{n}$, following the formalism in [7]. To this end, we need to define $C,f,D,G$ such that $\mathcal{H}=(C,f,D,G)$ provides a model of the PCOs network dynamics.

In particular, we define the flow map and the flow set as follows

[TABLE]

While the jump set can be defined as follows $D=\{x\in[0,2\pi]^{n}\colon\underset{i\in\mathcal{V}}{\max}\,x_{i}=2\pi\}$ so to enforce a jump whenever at one oscillator fires. Finally, to define the jump map $G$ of $\mathcal{H}$, we follow the idea in [13]. In particular, for each $x\in D$ we set $G(x)=(g_{1}(x),g_{2}(x),\dots,g_{n}(x))$ where for each $i\in\mathcal{V}$

[TABLE]

Notice that whenever one single oscillator fires, its state is reset to zero while the phase of each other oscillator is reset according to the function $Q$. If multiple oscillators fire, then the jump map is set-valued allowing to obtain a hybrid system satisfying the hybrid basic conditions. Moreover, to make $\mathcal{H}$ an accurate description of the evolution of the network of PCOs previously described, we consider the following standing assumption; see [19].

Assumption 1**.**

The function $Q\colon[0,2\pi]\rightarrow\mathbb{R}$ is such that $\operatorname{rge}(z\mapsto z+Q(z))\subset[0,2\pi]$ $\triangle$

The above assumption ensures that $G(D)\subset C\cup D=C$ avoiding the existence of solutions ending in finite time due to jumping outside $C\cup D$.

Remark 2**.**

Although the dynamics of the hybrid system $\mathcal{H}$ and the problem setting may look similar to the ones studied in [13], the problem addressed in this paper cannot be solved by employing the methodology proposed in [13]. Indeed, the desyncronized behavior studied in [13] (weak desynchronization) is defined as the behavior in which the separation between all of the firing events is equal (and nonzero). Obviously, the occurrence of this behavior does not ensure even spacing between phases; see the top figure in Figure 1. In fact, the reset rule adopted in [13] is such that the firing of a node makes all the other oscillators decrease their phase, hence preventing from achieving the splay state configuration333The simulation shown in the top of Figure 1 is performed by taking the parameter $\varepsilon$ in [13] equal to $-0.2$ and by picking an initial condition in the desynchronizaiton set given in [13].; see Figure 1. Therefore, further developments are needed to solve the splay state stabilization problem analyzed in this paper.

Remark 3**.**

In real-time decentralized implementations, whenever multiple oscillators fire (usually) each of the firing oscillators resets its phase to zero overlooking the incoming pulse stimulus received from other firing oscillators. This behavior is well known in the literature of PCOs and prevents to achieve convergence towards the splay state for every initial condition; see, e.g., [11]. It is interesting to observe that, although in the case of simultaneous firings our model gives rise to multiple solutions, the behavior described here above is still captured by our model. Indeed, if two or more oscillators fire simultaneously the definition of $G$ allows to consider (among the others) the solutions to $\mathcal{H}$ for which after one jump the phase of two or more oscillators are all equal to zero, preventing $\mathcal{H}$ from approaching the splay state configuration.

Now consider the following assumption.

Assumption 2**.**

The function $Q\colon[0,2\pi]\rightarrow\mathbb{R}$ satisfies the following properties:

• (A.1)

$Q$* is continuous on $[0,2\pi]$*

• (A.2)

$Q$* is such that $Q(2\pi)\neq 0$ $\triangle$*

Proposition 1**.**

Let Assumption 2 hold. Then, the following properties hold:

• (a)

$\mathcal{H}$* satisfies the hybrid basic conditions;*

• (b)

For each initial condition $\xi\in C\cup D$, there exists at least one nontrivial solution to $\mathcal{H}$. In particular, every $\phi\in\mathcal{S}_{\mathcal{H}}(C\cup D)$ is complete, non Zeno, and $\sup_{j}\operatorname{dom}\phi=\infty$.

Proof.

Point (a) is straightforward to prove. Indeed, $C$ and $D$ are closed in $\mathbb{R}^{n}$ and $f$ is continuous on $C$. Moreover, the jump map $G$ is obviously bounded and nonempty on $D$ and, similarly as shown in [13], is outer semicontinuous on $D$ since $Q$ is continuous444This latter property can be proven by directly using the definition of outer semicontinuity [7, Definition 5.9].. Hence, (a) holds. To prove (b), observe that since $\mathcal{H}$ satisfies the hybrid basic conditions, existence of nontrivial solutions from $C\cup D$ and completeness of each $\phi\in\mathcal{S}_{\mathcal{H}}(C\cup D)$ follows directly from the application of 555In particular, it can be verified that item (VC) in [7, Proposition 6.10] holds and that due to $G(D)\subset C\cup D$ and $C\cup D$ compact, items (b) and (c) in [7, Proposition 6.10] are ruled out.[7, Proposition 6.10]. Moreover, $\sup_{j}\operatorname{dom}\phi=\infty$ since $\phi$ is complete and the length of each flow interval is at most $\frac{2\pi}{\omega}$. To prove that each $\phi\in\mathcal{S}_{\mathcal{H}}(C\cup D)$ is non Zeno, we make use of [16, Lemma 2.7]. In particular, first observe that due to (A.2), since $Q(2\pi)\neq 0$, one has $G(D)\cap D=\emptyset$. Pick any $\phi\in\mathcal{S}_{\mathcal{H}}$ and let $E=\bigcup_{j\in\mathbb{N}_{0}}[t_{j},t_{j+1}]\times\{j\}$ be the domain of $\phi$. Then, since every maximal solution to $\mathcal{H}$ is precompact666A solution $\phi$ is precompact if it is bounded and complete. and, as shown earlier, $\mathcal{H}$ satisfies the hybrid basic conditions, thanks to [16, Lemma 2.7], it follows that there exists a positive real scalar $\gamma$ (dependent on the solution $\phi$) such that for each $j\geq 1$, $t_{j+1}-t_{j}\geq\gamma$ and this rules out the existence of Zeno solutions, concluding the proof. ∎

Remark 4**.**

The fact that $\mathcal{H}$ satisfies the hybrid basic conditions entails two main advantages. The first advantage is that working with a hybrid system fulfilling the hybrid basic conditions allows to exploit a large numbers of results available in the literature of hybrid systems such as the invariance principle in [16]. The second advantage is that the satisfaction of the hybrid basic conditions ensures that $\mathcal{H}$ is structurally robust to small perturbations, this aspect will be illustrated in Section III.b.

II-C Splay state stabilization

In this section, we formalize the splay state stabilization problem as a set stabilization problem for $\mathcal{H}$. In this paper, we consider the following notions of stability for the hybrid system $\mathcal{H}$:

Definition 5**.**

([7, Definition 7.1]) Let ${\cal A}\subset\mathbb{R}^{n}$ be a compact set. The set ${\cal A}$ is

• •

stable for $\mathcal{H}$ if for every $\epsilon>0$ there exists $\delta>0$ such that every solution to $\mathcal{H}$ with $|\phi(0,0)|_{\mathcal{A}}\leq\delta$ satisfies $|\phi(t,j)|_{\mathcal{A}}\leq\epsilon$ for all $(t,j)\in\operatorname{dom}\phi$;

• •

locally attractive for $\mathcal{H}$ if every maximal solution to $\mathcal{H}$ is complete and there exists $\mu\!>\!0$ such that for every maximal solution $\phi$ to $\mathcal{H}$ with $|\phi(0,0)|_{\mathcal{A}}\leq\mu$ one has $\lim_{t+j\rightarrow\!\infty\!}|\phi(t,j)|_{\mathcal{A}}\!=\!0$

• •

locally asymptotically stable (LAS) for $\mathcal{H}$, if it is both stable and locally attractive for $\mathcal{H}$.

Definition 6**.**

([7, Definition 7.3]) Let ${\cal A}\subset\mathbb{R}^{n}$ be a compact set locally asymptotically stable for $\mathcal{H}$. The basin of attraction of $\mathcal{A}$, denoted by $\mathcal{B}_{\mathcal{A}}$, is the set of points such that every $\phi\in\mathcal{S}_{\mathcal{H}}(\mathcal{B}_{\mathcal{A}})$ is bounded, complete, and $\lim_{t+j\rightarrow\infty}|\phi(t,j)|_{\mathcal{A}}=0$.

Remark 5**.**

Observe that the basin of attraction $\mathcal{B}_{\mathcal{A}}$ always contains points in $\mathbb{R}^{n}\setminus(C\cup D)$ since no solution exists from such points, implying that the completeness, boundedness, and convergence requirements are vacuously verified.

Definition 7**.**

(Almost global asymptotic stability) Let ${\cal A}\subset\mathbb{R}^{n}$ be a compact set asymptotically stable for $\mathcal{H}$. The set $\mathcal{A}$ is said to be almost globally asymptotically stable for $\mathcal{A}$ if $m(\mathbb{R}^{n}\setminus\mathcal{B}_{\mathcal{A}})=0$.

Remark 6**.**

Almost global asymptotic stability ensures that convergence towards the attractor $\mathcal{A}$ occurs for almost all initial conditions. Due to its relevance in practical applications, this property has been extensively used in the study of synchronization problems; see [4].

To make use of the above stability definition, analogously to [13], we show that the splay state configuration can be captured by defining an adequate compact set $\mathcal{A}\subset C$ such that whenever the state of $\mathcal{H}$ belongs to $\mathcal{A}$ then, the system is in the splay state configuration. To this aim, we first consider the following definition:

Definition 8**.**

(Splay state) Let $\phi$ be a solution to $\mathcal{H}$ and let for each $(t,j)\in\operatorname{dom}\phi$, $\psi(t,j)\coloneqq\operatorname{sort}(\phi(t,j))$. We say that $\phi$ is a splay state solution if there exists a strictly positive constant $c$, such that for each $(t,j)\in\operatorname{dom}\phi$ and each $i\in\mathcal{V}$, one has $\operatorname{d}(\psi_{i}(t,j),\psi_{\operatorname{next}_{\mathcal{V}}(i)}(t,j))=c$.

In particular, it can be easily argued that the constant $c$ in the above definition has to be equal to $\frac{2\pi}{n}$ to guarantee even spacing between the phases. This observation leads to the following definition of the set $\mathcal{A}$ characterizing the splay state configuration.

[TABLE]

Notice that $\mathcal{A}$ is compact since it is bounded (included in $C$ that is bounded) and closed since the function $\operatorname{sort}\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is continuous777Continuity of the mapping $x\mapsto\operatorname{sort}(x)$ can be easily shown by noticing that for each $x\in\mathbb{R}^{n}$ and any sequence $x_{k}\rightarrow x$, one has that $\operatorname{sort}(x_{k})\rightarrow\operatorname{sort}(x)$..

Remark 7**.**

Clearly, due to the difference between the splay state and the desynchronized behavior studied in [13], the set $\mathcal{A}$ in (3) is different from the desynchronization set defined in [13]. However, similarly to [13], the set $\mathcal{A}$ in (3) also enjoys a simple geometric characterization. In particular, it can be easily shown that the set $\mathcal{A}$ corresponds to the union of $n!$ segments888The fact that the set $\mathcal{A}$ can be represented as the union of $n!$ segments can be easily seen from the definition of the set $\mathcal{A}$ by observing that: for each $i\in\mathcal{V}$, the constraint on the geodesic distance $\operatorname{d}(\cdot)$ given in (3) gives rise to a segment in the space $[0,2\pi]^{n}$; the $n$ constraints considered in (3) can be described equivalently by considering all the orderings of the components of the vector $x$, which are $n!$. Thus, to determine such segments, for each of the considered orderings, it suffices to pick two points belonging to the set $\mathcal{A}$. in $[0,2\pi]^{n}$, which in turn can be written as the intersection of a set of lines $\mathcal{L}$ in $\mathbb{R}^{n}$ with the box $[0,2\pi]^{n}$. Nonetheless, in contrast to [13], in this paper this characterization is not useful to derive sufficient conditions for asymptotic stability of the set $\mathcal{A}$. This aspect is shown in Section III.

As already mentioned earlier, one cannot expect to make the set $\mathcal{A}$ attractive for every initial condition $\xi\in C$. Indeed, let

[TABLE]

Then, for each $\xi\in\mathcal{X}$ and for any choice of the function $Q$, there exists a solution $\phi$ to $\mathcal{H}$ with $\phi(0,0)=\xi$ which has identically synchronized components, i.e., $\phi$ does not approach $\mathcal{A}$. On other other hand, since the set $\mathcal{X}$ has measure zero, we solve the following problem:

Problem 1**.**

Design the function $Q$ such that the set $\mathcal{A}$ defined in (3) is almost globally asymptotically stable for the hybrid system $\mathcal{H}$.

III Main Results

III-A Nominal behavior

This first result relates the set $\mathcal{A}$ to another compact set whose role is key in the derivation of the main result presented in this paper.

Proposition 2**.**

Let $\mathcal{A}_{c}=\left\{x\in[0,2\pi]^{n}\colon\gamma(x)=2\pi\frac{n-1}{n}\right\}$, then $\mathcal{A}=\mathcal{A}_{c}.$

Proof.

First observe that for $n=2$ the result is trivial, so in the sequel we assume $n>2$. We prove the above claim in two parts. First we show that $\mathcal{A}\subset\mathcal{A}_{c}$. The above set inclusion is trivial to show. In particular, pick $x=(x_{1},x_{2},\dots,x_{n})\in\mathcal{A}$ and assume without loss of generality that $0=x_{1}<x_{2}<\dots<x_{n}$. Since for each $i\in\mathcal{V}$, $\operatorname{d}(x_{i},x_{\operatorname{next}_{\mathcal{V}}(i)})=\frac{2\pi}{n}$, one has $x_{2}=\frac{2\pi}{n},x_{3}=\frac{4\pi}{n},\dots,x_{n}=2\pi\frac{n-1}{n}$. Hence, it can be easily argued that $\gamma(x)=\sum_{i\in\mathcal{V}\setminus\{n\}}\operatorname{d}(x_{i},x_{\operatorname{next}_{\mathcal{V}}(i)})=2\pi\frac{n-1}{n}$, showing that $\mathcal{A}\subset\mathcal{A}_{c}$.

Now we show that $\mathcal{A}_{c}\subset\mathcal{A}$. Pick $x\in\mathcal{A}_{c}$, then $\gamma(x)=2\pi\frac{n-1}{n}>\pi$. Without loss of generality, assume $0=x_{1}\leq x_{2}\leq\dots\leq x_{n}$. Since $\gamma(x)\geq\pi$, as shown in [9], one has $\gamma(x)=2\pi-\underset{i\in\mathcal{V}}{\max}\operatorname{d}(x_{i},x_{\operatorname{next}_{\mathcal{V}}(i)})$, which, due to $x\in\mathcal{A}_{c}$, implies

[TABLE]

Suppose by contradiction that $x\notin\mathcal{A}$. Then, there exists $i\in\mathcal{V}$ such that $\operatorname{d}(x_{i},x_{\operatorname{next}_{\mathcal{V}}(i)})\neq\frac{2\pi}{n}$. Assume without loss of generality $i=1$. Thanks to (5), it has to be $\operatorname{d}(x_{1},x_{2})<\frac{2\pi}{n}$, which yields $x_{2}\in[0,\frac{2\pi}{n})$. Therefore, since $\gamma(x)>\pi$, by following the same steps as in the proof of Lemma 1, one can easily show that necessarily $x_{n}\in(\pi,2\pi\frac{n-1}{n})$, i.e., $\operatorname{d}(x_{n},x_{1})>\frac{2\pi}{n}$ contradicting (5). ∎

The above result shows that the set $\mathcal{A}$ corresponds to the set in which the length of the shortest containing arc is maximized. This fact is crucial and allows to exploit the notion of containing arc, in a dual fashion with respect to [10], for the construction of a Lyapunov function candidate for $\mathcal{H}$.

Now, with the aim of solving Problem 1, inspired by [11], we consider the following assumption on the function $Q$.

Assumption 3**.**

The function $Q$ is such that, the function $z\mapsto z+Q(z)$ is injective and $\operatorname{gph}Q\subset\mathcal{Q}_{1}\cup\mathcal{Q}_{2}$ where

[TABLE]

$\triangle$**

Remark 8**.**

Furthermore, notice that the satisfaction of Assumption 3 guarantees the satisfaction of Assumption 1, so the latter does not need to be checked.

Essentially, the above assumption guarantees that at each firing, the oscillators whose phase belongs to $({2\pi\frac{n-1}{n}},2\pi]$ decrease their phase without exiting $({2\pi\frac{n-1}{n}},2\pi]$. Instead, the oscillators whose phase belongs to $[0,{2\pi\frac{n-1}{n}}]$ are unaffected by the firing event. Moreover, having assumed $z\mapsto z+Q(z)$ to be injective guarantees that the firing of one oscillator does not make the phase of two or more oscillators jump to the same value. Indeed, in this case convergence towards $\mathcal{A}$ is impossible. This aspect and other important properties inherited by $\mathcal{H}$ from Assumption 2 and Assumption 3 are stated in the following result.

Lemma 2**.**

. Let $\mathcal{X}$ be defined as in (4) and let $\mathcal{U}=C\setminus\mathcal{X}$. If Assumption 2 and Assumption 3 hold, then for every $\phi\in\mathcal{S}_{\mathcal{H}}(\mathcal{U})$ one has $\operatorname{rge}\phi\subset\mathcal{U}$.

Proof.

Due to the definition of the flow dynamics, solution cannot exit the set $\mathcal{U}$ by flowing. Therefore, one needs to rule out the existence of solutions entering the set $\mathcal{X}$ via a jump. To this end, pick a solution $\phi\in\mathcal{S}_{\mathcal{H}}(\mathcal{U})$ and denote $(t,k)\mapsto\phi(t,k)\coloneqq(\phi_{1}(t,k),\phi_{2}(t,k),\dots,\phi_{n}(t,k))$. Assume that $(t,k)\in\operatorname{dom}\phi$ is such that $(t,k+1)\in\operatorname{dom}\phi$ and that $\phi(t,k)\in\mathcal{U}$. Since, $\phi(t,k)\in\mathcal{U}$, it follows that for each $i\neq j\in\mathcal{V}$ one has $\operatorname{d}(\phi_{i}(t,k),\phi_{j}(t,k))>0$. Let $\bar{i}\in\mathcal{V}$ such that $\phi_{\bar{i}}(t,k)=2\pi$, then $\phi_{\bar{i}}(t,k+1)=0$. While, for each $i\neq\bar{i}\in\mathcal{V}$, $\phi_{i}(t,k+1)=\phi_{i}(t,k)+Q(\phi_{i}(t,k))$. Thus, since $\phi(t,k)\in\mathcal{U}$ and $z\mapsto z+Q(z)$ is injective, for each $i\neq j\in\mathcal{V}\setminus\{\bar{i}\}$ $\phi_{i}(t,k+1)\neq\phi_{j}(t,k+1)$. To conclude the proof, one needs to show that for each $i\in\mathcal{V}\setminus\{\bar{i}\}$ $\phi_{i}(t,k+1)\notin\{0,2\pi\}$. But this is verified. Indeed, since $\phi(t,j)\in\mathcal{U}$, either $\phi_{i}(t,k)\in(0,{2\pi\frac{n-1}{n}}]$ or $\phi_{i}(t,k)\in({2\pi\frac{n-1}{n}},2\pi)$. In the first case, $\phi_{i}(t,k+1)=\phi_{i}(t,k)\notin\{0,2\pi\}$, while in the second case $\phi_{i}(t,k+1)\in[\phi_{i}(t,k)-\frac{2\pi}{n},\phi_{i}(t,k))\subset({2\pi\frac{n-1}{n}},2\pi)$, concluding the proof. ∎

Remark 9**.**

It is worthwhile to note that the proof of the above claim strongly relies on the fact that the oscillators are connected through an all-to-all network. Indeed, when such an assumption cannot be fulfilled, even if the function $z\mapsto z+Q(z)$ is injective, maximal solutions to $\mathcal{H}$ may leave the set $\mathcal{U}$ in finite time preventing from approaching $\mathcal{A}$. This fact makes the solution of the splay state stabilization problem for general topologies very challenging and unachievable with the simple interaction mechanism considered in this paper. Such an aspect is acknowledged by the existing literature on PCOs such as [3].

Now we give the main result of this paper.

Theorem 1**.**

Let Assumption 2 and Assumption 3 hold. Then, the set $\mathcal{A}$ defined in (3) is almost globally asymptotically stable for $\mathcal{H}$ and the basin of attraction of $\mathcal{A}$ includes $\mathcal{U}$.

Proof.

For every $x\in\mathbb{R}^{n}$, define the following Lyapunov function candidate $V(x)={2\pi\frac{n-1}{n}}-\gamma(x)$. Observe that $V$ is continuous in $\mathbb{R}^{n}$ (see [5]) and thanks to Lemma 1 and Proposition 2 positive definite with respect to $\mathcal{A}$ on $C\cup D$. Define, $u_{C}(x)=0$ for each $x\in C$ and $u_{C}(x)=-\infty$ otherwise; $u_{D}(x)=\max_{g\in G(x)}\{V(g)-V(x)\}$ for each $x\in D$ and $u_{D}(x)=-\infty$ otherwise. Observe that since during flows there is no interaction among the oscillators, due to the definition of the flow map, $V$ is constant during flows, hence the growth of $V$ is bounded by $u_{C},u_{D}$; see [7]. Let, $U=\mathbb{R}^{n}\setminus\mathcal{X}$, we want to show that for each $x\in U$, one has that $u_{D}(x)\leq 0$. To this end, one needs to show that $\gamma(x)$ is non decreasing across jumps from $x\in U\cap D$. That is, one needs to prove that for each $x\in U\cap D,g\in G(x)$, $\gamma(g)-\gamma(x)\geq 0$. Let $\Theta=\Psi([{2\pi\frac{n-1}{n}},2\pi])$ and for each $x\in[0,2\pi]^{n}$ let $\alpha(x)\subset\mathbb{S}^{1}$ be the shortest arc containing $x$. Three different situations need to be analyzed; see Figure 2. Since $G(x)$ is single-valued for $x\in U\cap D$, we denote $g(x)$ as the unique element of $G(x)$.

[Case 1 ]: $\alpha(x)\subset\Theta$; Figure 2 (a). Pick $x\in U\cap D$, and observe that necessarily for each $i\in\mathcal{V}$, $x_{i}\in[{2\pi\frac{n-1}{n}},2\pi]$ and that there exits a unique $k\in\mathcal{V}$ such that $x_{i}=2\pi$. Then due to Assumption 3 and due to the definition of $G$, one has that $g_{k}(x)=0$ and for each $l\neq k\in\mathcal{V}$, $g_{l}(x)\leq 0$. Hence, it can be argued that $\gamma(g(x))\geq\gamma(x)$.

[Case 2]: $\alpha(x)\supset\Theta$; Figure 2 (b). Pick $x\in U\cap D$, and observe that in this case $\gamma(g(x))=\gamma(x).$

[Case 3]: $\alpha(x)$ has an endpoint $h$ in $\Theta$; Figure 2 (c). Pick $x\in U\cap D$, and notice that either $\alpha(x)$ is the counterclockwise arc connecting $h$ to the other endpoint or $\alpha(x)\subset\Theta$, being the latter already discussed above. Then due to Assumption 3 and due to the definition of $G$, one has that there exists $k\in\mathcal{V}$, such that $g_{k}(x)=0$. Moreover, for each $l\neq k\in\mathcal{V}$ such that $x_{l}\in[0,{2\pi\frac{n-1}{n}}]$, one has $g_{l}(x)=0$, while for each $s\neq k\in\mathcal{V}$ such that $x_{s}\in({2\pi\frac{n-1}{n}},2\pi]$, one has $g_{s}(x)<0$. Thus, it can be argued that $\gamma(g(x))\geq\gamma(x)$. Therefore, for each $x\in U$ one has that $u_{D}(x)\leq 0$. Hence, since $U$ is a neighborhood of $\mathcal{A}$, $V$ is continuous and positive definite with respect to $\mathcal{A}$ on $C\cup D$, $G(D)\subset C\cup D$, and $\mathcal{H}$ satisfies the hybrid basic conditions, by [7, Theorem 8.8] it follows that $\mathcal{A}$ is stable for $\mathcal{H}$. Now, we show that $\mathcal{A}$ is locally attractive for $\mathcal{H}$ implying that such a set is LAS. To this end, we make use of [7, Theorem 8.8 (b)]. In particular, we want to show that for each $r\in(0,{2\pi\frac{n-1}{n}}]$ the largest weakly invariant999A set $S\subset\mathbb{R}^{n}$ is said to weakly invariant if for every $\xi\in S$, $\tau>0$, there exists at least one complete solution $\phi\in\mathcal{S}_{\mathcal{H}}(S)$ with $\operatorname{rge}\phi\subset S$ (weak forward invariance) such that for some $(t^{\star},j^{\star})\in\operatorname{dom}\phi$ with $t^{\star}+j^{\star}\geq\tau$ one has $\phi(t^{\star},j^{\star})=\xi$ (weak backward invariance); see [7] for further details on weak invariance in hybrid systems. subset contained in $V^{-1}(r)\cap U\cap(u_{C}^{-1}(0)\cup(u_{D}^{-1}(0)\cap G(u_{D}^{-1}(0)))=V^{-1}(r)\cap U\cap C=V^{-1}(r)\cap\mathcal{U}$ is empty. To this end, it suffices to show that for every $r\in(0,{2\pi\frac{n-1}{n}}]$, every $\phi\in\mathcal{S}_{\mathcal{H}}(V^{-1}(r)\cap\mathcal{U})$ leaves the set $V^{-1}(r)$. Pick any $\phi\in\mathcal{S}_{\mathcal{H}}(V^{-1}(r)\cap\mathcal{U})$, without loss of generality, suppose that $\phi(0,0)\in D$, and for each $(t,j)\in\operatorname{dom}\phi$, define $\psi(t,j)\colon\operatorname{sort}(\phi(t,j))$. Notice that, thanks to Assumption 3, if $\alpha(\xi)$ has an endpoint in $\Psi(({2\pi\frac{n-1}{n}},2\pi])$, then $V(\phi(0,1))<V(\phi(0,0))=r$. Hence, the only case to analyze is when the arc $\alpha(\phi(0,0))$ has an endpoint $h$ belonging to $\mathbb{S}^{1}\setminus\Psi(({2\pi\frac{n-1}{n}},2\pi])$. Now observe that, since $r\in(0,{2\pi\frac{n-1}{n}}]$, $\phi(0,0)\notin\mathcal{X}$, during flows the geodesic distance between oscillators’ phase is unchanged and, thanks to Assumption 3, at jumps it never increases more than $\frac{2\pi}{n}$, it follows that for each $(t,j)\in\operatorname{dom}\phi$, there exists $i\in\mathcal{V}$ such that $\operatorname{d}(\psi(t,j)_{i}(t,j),\psi(t,j)_{\operatorname{next}_{\mathcal{V}}(i)}(t,j))<\frac{2\pi}{n}$. Therefore, since from Assumption 3 at each jump all oscillators whose phase is in $({2\pi\frac{n-1}{n}},2\pi)$ decrease their phases and the endpoint $h$ of $\alpha(\phi(t,j))$ by assumption does not belong to $\Psi(({2\pi\frac{n-1}{n}},2\pi))$, it can be easily observed that for each $(t,j)\in\operatorname{dom}\phi$, there exist $(t_{s},j_{s}),(t_{s},j_{s}+1)\in\operatorname{dom}\phi$ such that $\alpha(\phi(t_{s},j_{s}))$ has an endpoint in $\Psi(({2\pi\frac{n-1}{n}},2\pi))$, implying that $\gamma(\phi(t_{s},j_{s}+1))>\gamma(\phi(t_{s},j_{s}))=r$. Hence, according to [7, Theorem 8.8], $\mathcal{A}$ is locally asymptotically stable for $\mathcal{H}$. To conclude the proof, we need to show that the basin of attraction of $\mathcal{A}$ includes $\mathcal{U}$. Indeed, in that case one has $\mathcal{B}_{\mathcal{A}}=(\mathbb{R}^{n}\setminus(\mathcal{U}\cup\mathcal{X}))\cup\mathcal{U}=\mathbb{R}^{n}\setminus\mathcal{X}$, which guarantees that $\mathcal{A}$ is almost globally asymptotically stable since $\mathbb{R}^{n}\setminus\mathcal{B}_{\mathcal{A}}=\mathcal{X}$ and 101010The set $\mathcal{X}$ is the union of a finite number of isolated points and affine subspaces of $\mathbb{R}^{n}$ of dimension $n-1$, i.e., sets of Lebesgue measure zero. Therefore, by countable subadditivity of the Lebesgue measure, one has $m(\mathcal{X})=0$; see [17]. $m(\mathcal{X})=0$. To this end, from [7, Theorem 8.2], one that each $\phi\in\mathcal{S}_{\mathcal{H}}(\mathcal{U})$ converges111111Maximal solutions to $\mathcal{H}$ are complete, bounded (precompact) and do not leave $C$ due to $C\cup D\cup G(D)=C$ compact., for some $r\in V(C)$, to the largest and nonempty weakly invariant subset of $V^{-1}(r)\cap C$. Then since we shown that for each $r\in(0,{2\pi\frac{n-1}{n}}]$, any $\phi\in\mathcal{S}_{\mathcal{H}}(\mathcal{U})$ leaves $V^{-1}(r)$, if follows that each $\phi\in\mathcal{S}_{\mathcal{H}}(\mathcal{U})$ approaches $\mathcal{A}$ and this finishes the proof. ∎

Remark 10**.**

The above result makes use of the notion of the shortest containing arc to ensure almost global asymptotic stability of the set $\mathcal{A}$. Nonetheless, due to the geometric characterization of the set $\mathcal{A}$ illustrated in Remark 7, it may be tempting to assume that the Lyapunov analysis proposed in [13] building on the function $x\mapsto\widetilde{V}(x)=|x|_{\mathcal{L}}$ (where the set $\mathcal{L}$ is defined in Remark (7)) can be adapted to deal with the setting considered in this paper. To this end, it is worthwhile to observe that in [13], to assess asymptotic stability of the desynchronization set via the adoption of the proposed distance-like function, the authors make use of the inherent relationships between the desynchronization set and the considered phase response function. However, due to different phase response function and set to be stabilized, such relationships do not hold in our setting. Thus, the arguments in [13] cannot be directly employed to assess asymptotic stability of the set $\mathcal{A}$ via the use of the proposed distance-like function. In fact, numerical experiments show that such a distance-like function, although approaching asymptotically zero, may increase along the solutions of the hybrid system $\mathcal{H}$ when the phase response function is designed according to the prescriptions given by Assumption 2 and Assumption 3, making such a function unsuitable as a Lyapunov candidate in our setting. This aspect is made evident in Section IV through numerical simulations.

III-B Robustness analysis

In real-world settings, due to parameter variations or constructive imperfections, assuming identical natural frequency $\omega$ for each oscillator is unrealistic. Therefore, having an insight on the behavior of the hybrid system $\mathcal{H}$ in the presence of small perturbations on the oscillators’ natural frequency is an interesting aspect. To delve into this issue, consider the following hybrid model of the network of PCOs with frequency perturbations

[TABLE]

where $d\in\mathbb{R}^{n}$ represents the perturbation affecting the flow dynamics. To formally describe the deviation of $\mathcal{H}_{d}$ from the nominal behavior captured by $\mathcal{H}$, we make use of the notion of $(\tau,\varepsilon)$-closeness of hybrid arcs given in Section I-B. In particular, by relying on the regularity of the data of $\mathcal{H}$, we can establish the following result.121212In the published version, there is a typo in the proof of Proposition 3. Namely, $\mathcal{H}_{\rho}$ should be $\mathcal{H}^{\prime}_{\rho}$ to avoid confusion with $\mathcal{H}_{d}$ in (6).

Proposition 3**.**

For each $\varepsilon>0$ and $\tau\geq 0$, there exists $\delta>0$ such that for each measurable function $d\colon\mathbb{R}_{\geq 0}\rightarrow\delta\mathbb{B}$ and for each solution131313Given a measurable function $d\colon\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{n}$, $\phi_{d}$ is a solution to $\mathcal{H}_{d}$ if $(\phi,d^{\prime})$ with $d^{\prime}(t,j)=d(t)$ for each $(t,j)\in\operatorname{dom}\phi$ is a solution pair to $\mathcal{H}_{d}$; see [1] for further details. $\phi_{d}$ to (6) from $C$, there exists a solution $\phi$ to $\mathcal{H}$ from $C$ such that $\phi$ and $\phi_{d}$ are $(\tau,\varepsilon)$-close.

Proof.

Let $\rho$ be a positive real scalar, define $\mathcal{H}^{\prime}_{\rho}=(C,F_{\rho},D,G)$, where for each $x\in C$ $F_{\rho}(x)=f(x)+\rho\mathbb{B}$. Since $\mathcal{H}$ satisfies the hybrid basic conditions and is pre-forward complete141414Given a set $S$, a hybrid system $\mathcal{H}$ is said to be pre-forward complete from $S$ if every $\phi\in\mathcal{S}_{\mathcal{H}}(S)$ is either complete or bounded; see [7, Definition 6.12]. from the (compact) set $C$, from [7, Proposition 6.34] it follows that for every $\varepsilon>0$ and $\tau\geq 0$ there exists $\delta>0$ such that for every solution $\phi_{\delta}$ to $\mathcal{H}^{\prime}_{\delta}$ there exists a solution $\phi$ to $\mathcal{H}$ such that $\phi$ and $\phi_{\delta}$ are $(\tau,\varepsilon)$-close. To conclude, if suffices to notice that whenever $d(t)\in\delta\mathbb{B}$ for each $t\geq 0$, solutions to $\mathcal{H}_{d}$ are solutions to {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathcal{H}^{\prime}_{\delta}} and this concludes the proof. ∎

Essentially, the above result states that in the presence of small perturbations on the natural frequency, the evolution of the perturbed network of PCOs does not differ too much from the one of the unperturbed network. Therefore, one can expect that the convergence towards the splay state is practically preserved in the case of (small) frequency perturbations.

IV Numerical Example

IV-A Nominal behavior

In this example in which $n=3$, we want to show both the effectiveness of the proposed methodology and the differences between the dynamics of $\mathcal{H}$ and the results given in [13]. We use the following phase response function

[TABLE]

Figure 3 shows the evolution of a solution $\phi$ to $\mathcal{H}$ starting from an initial condition $\phi(0,0)\notin\mathcal{A}$. As expected, as $t+j$ goes to infinity, $\phi$ approaches the splay state configuration and the Lyapunov function $V$ is constant during flows, nonincreasing across jumps, and approaches [math].

Now, as mentioned earlier, we show that the Lypunov like function proposed in [13], i.e., $x\mapsto\widetilde{V}(x)=|x|_{\mathcal{L}}$, cannot be employed in the setting analyzed in this paper. To this end, in Figure 4 we report the evolution of $\widetilde{V}$ along the solution $\phi$ to $\mathcal{H}$ analyzed above. Simulations show that, although $\widetilde{V}$ approaches zero as $t+j$ goes to infinity, the considered function increases across some jumps and this makes it unsuitable as a Lyapunov candidate to assess asymptotic stability of the set $\mathcal{A}$ for the hybrid system $\mathcal{H}$.

IV-B Perturbed behavior

Now, we want to analyze the effect of small frequency perturbation on the dynamics of $\mathcal{H}$. In particular, to study the effect of uncorrelated periodic frequency perturbations, we select $d(t)=\epsilon(\sin(0.5t),\sin(0.5t+2/3\pi),\sin(0.5t+4/3\pi))$, where $\epsilon>0$ represents the maximum amplitude of the frequency variation. Figure 5 reports the evolution of $\mathcal{H}_{d}$ and of $\mathcal{H}$ for $\epsilon$ equal to $0.05$ and $0.03$, repsectively. As expected, the evolution of $\mathcal{H}_{d}$ does not differ too much from the evolution of $\mathcal{H}$. Figure 6 shows the evolution of the function $V$ along, respectively, the evolution of $\mathcal{H}$ and of $\mathcal{H}_{d}$ for the two considered values of $\epsilon$. Simulations show that the function $V$ in the presence of small frequency perturbations approaches a ball containing zero. In that simulation, one can remark that the larger the amplitude of the perturbation, the larger the radius of the ball approached by $V$.

V Conclusion

In this paper, we studied splay state almost global asymptotic stabilization in pulse coupled oscillators. The considered problem was turned into the stabilization problem of a compact set $\mathcal{A}$ wherein oscillators’ phases are evenly distributed on the unit circle. Sufficient conditions on the phase response function to guarantee almost global asymptotic stability of set $\mathcal{A}$ were provided. In particular, almost global asymptotic stability of the set $\mathcal{A}$ is assessed via the use of a novel Lyapunov function along with the invariance principle for hybrid systems in [16]. Moreover, the proposed approach was shown to be robust to small frequency perturbations naturally present in practice. The results presented in this technical note are promising and the framework adopted quite flexible to envision interesting extensions of the results presented here. Among them, we mention the extension to more complex network topologies.