---
abstract: 'The property of desynchronization in an all-to-all network of homogeneous impulse-coupled oscillators is studied. Each impulse-coupled oscillator is modeled as a hybrid system with a single timer state that self-resets to zero when it reaches a threshold, at which event all other impulse-coupled oscillators adjust their timers following a common reset law. In this setting, desynchronization is considered as each impulse-coupled oscillator’s timer having equal separation between successive resets. We show that, for the considered model, desynchronization is an asymptotically stable property. For this purpose, we recast desynchronization as a set stabilization problem and employ Lyapunov stability tools for hybrid systems. Furthermore, several perturbations are considered showing that desynchronization is a robust property. Perturbations on both the continuous and discrete dynamics are considered. Numerical results are presented to illustrate the main contributions.'
author:
- 'Sean Phillips and Ricardo G. Sanfelice[^1]'
bibliography:
- 'Biblio.bib'
- 'RGS.bib'
- 'SAP.bib'
---

Introduction
============

Impulse-coupled oscillators are multi-agent systems with state variables consisting of timers that evolve continuously until a state-dependent event triggers an instantaneous update of their values. Networks of such oscillators have been employed to model the dynamics of a wide range of biological and engineering systems. In fact, impulse-coupled oscillators have been used to model groups of fireflies [@Mirollo.90.SIAMJAM.BiologicalOscillators], spiking neurons [@Pikovsky.ea.03.BiologicalOscillators; @gerstner2002spiking], muscle cells [@Peskin.75.BiologicalOscillators], wireless networks [@hong2010cooperative], and sensor networks [@Liu05adynamic]. With synchronization being a property of particular interest, such complex networks have been found to coordinate the values of their state variables by sharing information only at the times the events/impulses occur [@Mirollo.90.SIAMJAM.BiologicalOscillators; @Abbott1993].

The opposite of synchronization is [*desynchronization*]{}. In simple words, desynchronization in multi-agent systems is the notion that the agents’ periodic actions are separated “as far apart” as possible in time. Desynchronization is similar to clustering or splay-state configurations, and is sometimes referred in the literature as inhibited behavior [@mauroy:037122; @glass1988clocks]. For impulse-coupled oscillators, desynchronization is given as the behavior in which the separation between all of the timers impulses is equal [@Patel.Desync2007]. This behavior has been found to be present in communication schemes in fish [@Benda.Neuron] and in networks of spiking neurons [@Pfurtscheller19991842; @1997Natur.390.70S]. Desynchronization of oscillators has recently been shown to be of importance in the understanding of Parkinson’s disease [@Mabi.Moehlis2010; @Majtanik.Dolan2004], in the design of algorithms that limit the amount of overlapping data transfer and data loss in wireless digital networks [@hong2010cooperative], and in the design of round-robin scheduling schemes for sensor networks [@Liu05adynamic].

Motivated by the applications mentioned above and the lack of a full understanding of desynchronization in multi-agent systems, this paper pertains to the study of the dynamical properties of desynchronization in a network of impulse-coupled oscillators with an all-to-all communication graph. The uniqueness of the approach emerges from the use of hybrid systems tools, which not only conveniently capture the continuous and impulsive behavior in the networks of interest, but also are suitable for analytical study of asymptotic stability and robustness to perturbations.

More precisely, the dynamics of the proposed hybrid system capture the (linear) continuous evolution of the states as well their impulsive/discontinuous behavior due to state triggered events. Analysis of the asymptotic behavior of the trajectories (or solutions) to these systems is performed using the framework of hybrid systems introduced in [@teel2012hybrid; @Goebel.ea.09.CSM]. To this end, we recast the study of desynchronization as a set stabilization problem. Unlike synchronization, for which the set of points to stabilize is obvious, the complexity of desynchronization requires first to determine such a collection of points, which we refer to as the [*desynchronization set*]{}. We propose an algorithm to compute such set of points. Then, using Lyapunov stability theory for hybrid systems, we prove that the desynchronization set is asymptotically stable by defining a Lyapunov-like function as the distance between the state and (an inflated version of) the desynchronization set. In our context, asymptotic stability of the desynchronization set implies that the distance between the state and the desynchronization set converges to zero as the amount of time and the number of jumps get large. Using the proposed Lyapunov-like function and invoking an invariance principle, the basin of attraction is characterized and shown to be the entire state space minus a set of measure zero, which turns out to actually be an exact estimate of the basin of attraction. Furthermore, also exploiting the availability of a Lyapunov-like function, we analytically characterize the time for the solutions to reach a neighborhood of the desynchronization set. In particular, this characterization provides key insight for the design of algorithms used in applications in which desynchronization is crucial, such as wireless digital networks and sensor networks.

The asymptotic stability property of the desynchronization configuration is shown to be robust to several types of perturbations. The perturbations studied here include a generic perturbation in the form of an inflation of the dynamics of the proposed hybrid system model of the network of interest and several kinds of perturbations on the timer rates. Using the tools presented in [@teel2012hybrid; @Goebel.ea.09.CSM], we analytically characterize the effect of these perturbations on the already established asymptotic stability property of the desynchronization set. In particular, these perturbations capture situations where the agents in the network are heterogeneous due to having differing timer rates, threshold values, and update laws. To verify the analytical results, we simulate networks of impulse-coupled oscillators under several classes of perturbations. Specifically, we show numerical results when perturbations affect the update laws and the timer rates.

The remainder of this paper is organized as follows. Section \[sec:hs\] is devoted to hybrid modeling of networks of impulse-coupled oscillators. Section \[sec:AN\] introduces an algorithm to determine the desynchronization set. Section \[sec:lyapunov\] presents the stability results while the time to convergence is characterized in Section \[sec:timetoconverge\]. The robustness results are in Section \[sec:robustness\]. Section \[sec:numerics\] presents numerical results illustrating our results. Final remarks are given in Section \[sec:conclusion\].

[**Notation**]{}

Hybrid System Model of Impulse-Coupled Oscillators {#sec:hs}
==================================================

Mathematical Model {#sec:hybridmodel}
------------------

In this paper, we consider a model of $N$ impulse-coupled oscillators. Each impulse-coupled oscillator has a continuous state ($\tau_i$ for the $i$-th oscillator) defining its internal timer. Once the timer of any oscillator reaches a threshold (${\bar{\tau}}$), it triggers an impulse and is reset to zero. At such an event, all the other impulse-coupled oscillators rescale their timer by a factor given by $(1+ \varepsilon)$ times the value of their timer, where $\varepsilon \in (-1,0)$.[^2] Figure \[fig:example\] shows a trajectory of two impulse-coupled oscillators [with states]{} ${\tau_1}$ and ${\tau_2}$. In this figure, the dark red circles indicate when a timer state has reached the threshold and, thus, resets to zero. The light green circles indicate when an oscillator is externally reset and, hence, decreases its timer by $(1+\varepsilon)$ times its current state.

According to this outline of the model, the dynamics of the impulse-coupled oscillators involve impulses and timer resets, which are treated as true discrete events and instantaneous updates, while the smooth evolution of the timers before/after these events define the continuous dynamics. We follow the hybrid formalism of [@teel2012hybrid; @Goebel.ea.09.CSM], where a hybrid system is given by four objects $(C,f,D,G)$ defining its *data*:

-   *Flow set:* a set $C \subset {\mathbb{R}}^{N}$ specifying the points where flows are possible (or continuous evolution).

-   *Flow map:* a single-valued map $f: {\mathbb{R}}^{N} \to {\mathbb{R}}^{N}$ defining the flows.

-   *Jump set:* a set $D \subset {\mathbb{R}}^{N}$ specifying the points where jumps are possible (or discrete evolution).

-   *Jump map:* a set-valued map $G: {\mathbb{R}}^{N} \rightrightarrows{\mathbb{R}}^{N}$ defining the jumps.

A hybrid system capturing the dynamics of the impulse-coupled oscillators is denoted as ${{{\mathcal}H}}_N := (C,f,D,G)$ and can be written in the compact form $${{{\mathcal}H}}_N: \qquad \tau \in {\mathbb{R}}^{N} \qquad \left\{ \begin{array}{llll} \dot{\tau} &=& f(\tau) &\quad \tau \in C \\ \tau^{+} &\in& G(\tau) & \quad \tau \in D \end{array}\right. , \label{eqn:HS}$$ where $N \in {\mathbb{N}}\setminus \{0,1\}$ is the number of impulse-coupled oscillators. The state of ${{{\mathcal}H}}_{N}$ is given by The flow and jump sets are defined to constrain the evolution of the timers. The flow set is defined by $$\startmodif C := P_N, \stopmodif
\label{eqn:flawiest}$$ where $I := \{1,2, \ldots , N\}$ and ${\bar{\tau}}> 0$ is the threshold. During flows, an internal clock gradually increases based on the homogeneous rate, $\omega$. Then, the flow map is defined as with $\omega > 0$ [[ defining the natural frequency of each impulse-coupled oscillator]{}]{}. The impulsive events are captured by a jump [[ set $D$ and a jump]{}]{} map $G$. Jumps occur when the state is in the jump set $D$ defined as $$D := \left\{ \tau \in P_N : \ \exists i \in I \ \mbox{s.t.} \ \tau_i = {\bar{\tau}}\right\} . 
\label{eqn:D}$$ From such points, the $i$-th timer is reset to zero and forces a jump of all other timers. Such discrete dynamics are captured by the following jump map: for each $\tau \in D$ define $
G(\tau) = \left[ g_1(\tau) \ \ g_2(\tau)\ \ \ldots \ \ g_N(\tau)  \right]^\top,
$ where, for each $i \in I$, $$g_i (\tau) = \left\{ \begin{array}{l} 0 \qquad \qquad \ \ \ \ \mbox{if } \tau_{i} = {\bar{\tau}}, \tau_r < {\bar{\tau}}\ \ \forall r \in I\setminus \{i\} \\ 
\{0 , \tau_{i}(1+\varepsilon) \} \ \mbox{if } \tau_{i} = {\bar{\tau}}\ \exists r \in I \setminus \{i\} \  \mbox{s.t.} \  \tau_r = {\bar{\tau}}\\ 
(1+\varepsilon)\tau_{i} \qquad \  \mbox{if } \tau_i < {\bar{\tau}}\ \exists r \in I \setminus \{i\} \  \mbox{s.t.} \  \tau_r = {\bar{\tau}}\end{array} \right. \label{eqn:gi}$$ with parameters $\varepsilon \in (-1,0)$ and ${\bar{\tau}}> 0$; for $\tau \in D$, $g_i$ is not empty. When a jump is triggered, the state $\tau_i$ jumps according to the $i$-th component of the jump map $g_i$. When a state reaches the threshold ${\bar{\tau}}$, it is reset to zero only when all other states are less than that threshold; otherwise, if multiple timers reach the threshold simultaneously, the jump map is set valued to indicate that either $g_i(\tau) = 0$ or $g_i(\tau) = (1+\varepsilon)\tau_{i}$ is possible. This is to ensure that the jump map satisfies the regularity conditions outlined in Section \[sec:ModelForAnalysis\].[^3]

For example, consider the case $N=2$ the hybrid system ${{{\mathcal}H}}_N = (C,f,D,G)$ has state given by $$\tau=\left[ \begin{array}{c} \tau_{1} \\ \tau_{2}\end{array} \right] \in P_2 := [0,{\bar{\tau}}]\times[0,{\bar{\tau}}] .$$ The states $\tau_{1}$ and $\tau_{2}$ are the timers for both of the oscillators. The hybrid system ${{{\mathcal}H}}_2$ has the following data:

$${{{\mathcal}H}}_2 = \left\{\begin{array}{ll} C = P_2, & \qquad f(\tau) = \left[ \begin{array}{c} 1
\\ 1 \end{array}\right] \forall \tau \in C, \\ \noalign{\medskip} 
D = \left\{ \tau \in P_2 \ : \ \exists i \in \{1,2\} \ s.t. \ \tau_{i} = {\bar{\tau}}\right\}, & \qquad
G(\tau) = \left[ \begin{array}{c} g_1(\tau) \\ g_2(\tau) \end{array} \right] \forall \tau \in D , \end{array} \right.$$ where the functions $g_1$ and $g_2$ are defined as $$g_1 (\tau) = \left\{ \begin{array}{ll}  0 & \mbox{if } \tau_{1} = {\bar{\tau}}, \tau_2 < {\bar{\tau}}\\ \{0 , \tau_{1}(1+\varepsilon) \} & \mbox{if } \tau_{1} = {\bar{\tau}}, \tau_2 = {\bar{\tau}}\\
(1+\varepsilon)\tau_{1} & \mbox{if } \tau_1 < {\bar{\tau}}, \tau_2 = {\bar{\tau}}\end{array} \right.
\qquad \qquad  g_2 (\tau) = \left\{ \begin{array}{ll}  0 & \mbox{if } \tau_{2} = {\bar{\tau}}, \tau_1 < {\bar{\tau}}\\ \{0 , \tau_{2}(1+\varepsilon) \} & \mbox{if } \tau_{2} = {\bar{\tau}}, \tau_1 = {\bar{\tau}}\\
(1+\varepsilon)\tau_{2} & \mbox{if } \tau_2 < {\bar{\tau}}, \tau_1 = {\bar{\tau}}\end{array} \right. .$$

Basic Properties of ${{{\mathcal}H}}_{N}$ {#sec:BasicCond}
-----------------------------------------

\[sec:ModelForAnalysis\]

### Hybrid Basic Conditions

To apply analysis tools for hybrid systems in [@teel2012hybrid], which will be summarized in Section \[sec:analysis\], the data of the hybrid system ${{{\mathcal}H}}_{N}$ must meet certain mild conditions. These conditions, referred to as the [*hybrid basic conditions*]{}, are as follows:

1.  $C$ and $D$ are closed sets in ${\mathbb{R}}^N$.

2.  $f: {\mathbb{R}}^N\to{\mathbb{R}}^N$ is continuous on $C$.

3.  $G :{\mathbb{R}}^N{\rightrightarrows}{\mathbb{R}}^N$ is an outer semicontinuous[^4] set-valued mapping, locally bounded on $D$, and such that $G(x)$ is nonempty for each $x \in D$.

\[lem:BasicConds\] ${{{\mathcal}H}}_{N}$ satisfies the hybrid basic conditions.

Note that satisfying the hybrid basic conditions [[ implies that ${{{\mathcal}H}}_N$ is well-posed [@teel2012hybrid Theorem 6.30], which automatically gives robustness to vanishing state disturbances; see [@teel2012hybrid; @Goebel.ea.09.CSM]. Section \[sec:robustness\] considers different types of perturbations that ${{{\mathcal}H}}_N$ can withstand. ]{}]{}

### Solutions to ${{{\mathcal}H}}_N$

\[lem:solutions2HS\] From every point in $C \cup D$, there exists a solution and every maximal solution to ${{{\mathcal}H}}_{N}$ is complete and bounded.

Due to the jump map $G$, if the elements of the solution are initially equal (denote this set as ${\mathcal{S}}:= \{\tau \in P_N : \exists i,r \in I, i \neq r, \tau_{i} = \tau_r \}$) it is possible for them to remain equal for all time. Furthermore, it is also possible for solutions to be initialized on the jump set such that one element is at the threshold and another is equal to zero then after the jump they will be equal, e.g. let $\tau_{1} = {\bar{\tau}}$, $\tau_{2} = 0$ then $\tau_{1}^{+} = \tau_{2}^{+} = 0$. We denote this set as ${\mathcal{G}}:= \{\tau \in D\setminus {\mathcal{S}}: \exists i,r \in I, i \neq r, \tau_{i} = 0, \tau_r = {\bar{\tau}}\}$. The next result considers solutions initialized on the set ${\mathcal{X}}:= {\mathcal{S}}\cup {\mathcal{G}}$.

\[lem:defX\] For each $\tau(0,0) \in {\mathcal{X}}$, there exists a solution $\tau$ to ${{{\mathcal}H}}_{N}$ from $\tau(0,0)$ such that, for some $M \in \{0,1\}$, $\tau(t,j) \in {\mathcal{S}}$ for all $t+j \geq M$, $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$.

Consider a solution $\tau$ to the hybrid system ${{{\mathcal}H}}_N$ with initial condition $\tau(0,0) \in \startmodif {\mathcal{S}}\stopmodif$. Due to the flow map for each state being equal, $\tau$ remains in ${\mathcal{S}}$ during flows. Furthermore, at points $\tau \in {\mathcal{S}}\cap D$, the jump map $G$ is set valued by the definition of $g_i$ in . From these points, $G(\tau) \cap {\mathcal{S}}\neq \emptyset$. In fact, for each $\tau(0,0) \in {\mathcal{S}}$, there exists at least one solution such that $\tau(t,j) \in {\mathcal{S}}$ for all $t + j \geq 0$, with $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$. Consider the case of solutions initialized at $\tau(0,0) \in \startmodif {\mathcal{G}}\stopmodif$ (Note that $\tau(0,0) \in D$). It follows that for some $r \in I$, $\tau_r(0,0) = {\bar{\tau}}$ and $g_r(\tau(0,0)) = 0$. Therefore, after the initial jump, we have that $G(\tau(0,0)) \cap {\mathcal{S}}\neq \emptyset,$ by which using previous arguments implies that $\tau(t,j) \in {\mathcal{S}}$ for all $t + j \geq 1$.

Furthermore, there is a distinct ordering to the jumps. If $\tau$ is such that $\tau_i \neq \tau_r$ for all $i\neq r$ then the ordering of each $\tau_i$ is preserved after $N$ jumps. More specifically, we have the following result.

\[lem:ordering\] For every solution $\tau$ to ${{{\mathcal}H}}_N$ with $\tau(0,0) \notin {\mathcal{X}}$, if at $(t_{j},j) \in {\mathop{\rm dom}\nolimits}\tau$ we have for some sequence of nonrepeated elements $\{i_m\}^N_{m = 1}$ of $I$ (that is, a reordering of the elements of the set $I = \{1,2,\ldots,N\}$) then, after $N$ jumps, it follows that

Let $\tau$ be a solution to ${{{\mathcal}H}}_{N}$ from $P_{N}\setminus {\mathcal{X}}$. There exists a sequence $i_{k}$ of distinct elements with $i_{k} \in I$ for each $k \in I$, such that $0 \leq \tau_{i_{1}}(t,j) < \tau_{i_{2}}(t,j) < \ldots < \tau_{i_{N}}(t,j) \leq {\bar{\tau}}$ over $[t_{0},t_{1}]\times \{0\}$. After the jump at $(t,j) = (t_{1},0)$ we have $0 = \tau_{i_{N}}(t,j+1) < \tau_{i_{1}}(t,j+1) < \tau_{i_{2}}(t,j+1) < \ldots <  \tau_{i_{N-1}}(t,j+1) < {\bar{\tau}}$. Continuing this way for each jump, it follows that after $N-1$ more jumps, the solution is such that $0 \leq \tau_{i_{1}}(t_{N},j+N) < \tau_{i_{2}}(t_{N},j+N) < \ldots < \tau_{i_{N}}(t_{N},j+N) \leq {\bar{\tau}}$ and the order at time $(t,j)$ is preserved.

Using these properties of solutions to ${{{\mathcal}H}}_N$, the next section defines the set to which these solutions converge and establishes its stability properties.

Dynamical Properties of ${{{\mathcal}H}}_{N}$ {#sec:analysis}
=============================================

Our goal is to show that the desynchronization configuration [[ of ${{{\mathcal}H}}_N$, which is defined in Section \[sec:AN\]]{}]{}, is asymptotically stable. We recall from [@teel2012hybrid; @Goebel.ea.09.CSM] the following definition of asymptotic stability for general hybrid systems with state $x \in {\mathbb{R}}^n$.

\[def:AS\] A closed set ${\mathcal{A}}\subset {\mathbb{R}}^n$ is said to be

-   [*stable*]{} if for each $\varepsilon>0$ there exists $\delta>0$ such that each solution $x$ with $|x(0,0)|_{{\mathcal{A}}}\leq \delta$ satisfies $|x(t,j)|_{{\mathcal{A}}} \leq \varepsilon$ for all $(t,j)\in{\mathop{\rm dom}\nolimits}x$;

-   [*attractive*]{} if there exists $\mu > 0$ such that every maximal solution $x$ with $|x(0,0)|_{{\mathcal{A}}}\leq \mu$ is complete and satisfies\
    $\lim_{(t,j) \in {\mathop{\rm dom}\nolimits}x, t+j\to\infty} |x(t,j)|_{{\mathcal{A}}}=0$;

-   [*asymptotically stable*]{} if stable and attractive;

-   [*weakly globally asymptotically stable*]{} if ${\mathcal{A}}$ is stable and if, for every initial condition, there exists a maximal solution that is complete and satisfies $\lim_{(t,j) \in {\mathop{\rm dom}\nolimits}x, t+j\to\infty} |x(t,j)|_{{\mathcal{A}}}=0$.

The set of points from where the attractivity property holds is the basin of attraction and excludes all points where the system trajectories may never converge to ${\mathcal{A}}$. In fact, it will be established in Section \[sec:lyapunov\] that the basin of attraction for asymptotic stability of desynchronization of ${{{\mathcal}H}}_N$ does not include any point $\tau$ such that any two or more timers are equal or become equal after a jump, which is the set ${\mathcal{X}}$ defined in Lemma \[lem:defX\].

Construction of the set ${\mathcal{A}}$ for ${{{\mathcal}H}}_N$ {#sec:AN}
---------------------------------------------------------------

In this section, we identify the set of points corresponding to the impulse-coupled oscillators being desynchronized, namely, we define the [*desynchronization set*]{}. We define desynchronization as the behavior in which the separation between all of the timers’ impulses is equal (and nonzero), see Figure \[fig:example\]. More specifically desynchronization is defined as follows:

\[def:desynch\] A solution $\tau$ to ${{{\mathcal}H}}_N$ is desynchronized if there exists $\Delta > 0$ and a sequence of non-repeated elements $\{i_m\}^N_{m = 1}$ of $I$ (that is, a reordering of the elements of the set $I = \{1,2, \ldots, N \}$) such that $\lim_{j \to \infty} (t_j^{i_m} - t_j^{i_{m+1}} )= \Delta$ for all $m \in \{1,2, \ldots, N-1\}$ and $\lim_{j \to \infty} (t_j^{N} - t_j^{i_{1}}) = \Delta,$ where $\{t_j^{i_m}\}_{j = 0}^\infty$ is the sequence of jump times of the state $\tau_{i_m}$.

In fact, this separation between impulses leads to an ordered sequence of impulse times with equal separation. The desynchronization set ${\mathcal{A}}$ for the hybrid system ${{{\mathcal}H}}_N$ captures such a behavior and is parameterized by $\varepsilon$, the threshold ${\bar{\tau}}$, and the number of impulse-coupled oscillators $N$.

To define this set, first we provide some basic intuition about the dynamics of ${{{\mathcal}H}}_N$ when desynchronized. The set ${\mathcal{A}}$ must be forward invariant and such that trajectories staying in it satisfy the property in Definition \[def:desynch\]. Due to the definition of the flow map $f$, there exist sets in the form of “lines" $\ell_k$, each of them in the direction ${\mathbf 1}$, which is the direction of the flow map, intersecting the jump set at a point which, for the $k$-th line, we denote as ${\widetilde}\tau^{k}$. We define the desynchronization set as the union of sets $\ell_{k}$ collecting points $\tau = {\widetilde}\tau^k + {\mathbf 1}s \in P_N$ parameterized by $s \in {\mathbb{R}}$.

To identify ${\widetilde}\tau^k$, consider a point ${\widetilde}\tau^{k} \in D\setminus {\mathcal{X}}$ with components satisfying ${\widetilde}\tau_{1}^{k} = {\bar{\tau}}> {\widetilde}\tau_{2}^{k} > {\widetilde}\tau_{3}^{k} > ... > {\widetilde}\tau_{N}^{k}$. Due to Definition \[def:desynch\], it must be true that the difference between jump times are constant. This means that there must be some correlation between $\Delta$ and the difference between, in this case, $\tau_{1}^{k}$ and $\tau_{2}^{k}$. Moreover, there must be a correlation between $\tau_{1}^{k}$ and all other states at jumps. It follows that this point belongs to ${\mathcal{A}}$ only if the distance between the expiring timer (${\widetilde}\tau_{1}^{k}$) and each of its other components (${\widetilde}\tau_{i}^{k}$, $i \in I \setminus \{1\}$) is equal to the distance between the value after the jump of the timer expiring next (${\widetilde}\tau_{2}^{k}\null^{+}$) and the value after the jump of its other components (${\widetilde}\tau_{i}^{k}\null^{+}$, $i \in I \setminus \{2\}$), respectively. This property ensures that, when in the desynchronization set, the relative distance between the leading timer and each of the other timers is equal, before and after jumps. More precisely, $$\begin{aligned}
{\widetilde}\tau_{1}^{k} - \widetilde{\tau}_{i}^{k} = {\widetilde}\tau_{2}^{k}\null^{+} - {\widetilde}\tau_{\mbox{\scriptsize next}(i)}^{k}\hspace{-0.21in}\null^{+} \qquad
\qquad \forall \ i \in I \setminus \{1\} \label{eqn:ANcondition},\end{aligned}$$ where ${\widetilde}\tau^{k}\null^{+} = G({\widetilde}\tau^{k})$ and next$(i) = i + 1$ if $i + 1 \leq N$ and $1$ otherwise.[^5] Since ${\mathcal{X}}$ contains all points such that at least two or more timers are the same, we can consider the case when one component of ${\widetilde}\tau^k$ is equal to ${\bar{\tau}}$ at a time. For each such case, we have $(N - 1)!$ possible permutations of the other components and $N$ possible timer components equal to ${\bar{\tau}}$, leading to $N!$ total possible sets $\ell_{k}$.

For the $N$ case, the algorithm above results in the system of equations $\Gamma \tau_s = b$, where and $
b = \bar{\tau} {\bf 1},
$ where $\tau_s$ is the state ${\widetilde}\tau^{k}$ sorted into decreasing order. It can be shown that for any $\varepsilon \in (-1,0)$, a solution $\tau_{s}$ exists (see Lemma \[eqn:tauSsolution\]). Then, $\tau_s$ needs to be unsorted and becomes ${\widetilde}\tau^k$ in the definition of the set $\ell_k$.

The solution to $\Gamma\tau_s = b$ is the result of a single case of $\tau \in D \setminus {\mathcal{X}}$. As indicated above, to get a full definition of the set ${\mathcal{A}}$, the $N!$ sets $\ell_k$ should be computed. For arbitrary $N$, the set ${\mathcal{A}}$ is given as a collection of sets $\ell_{k}$ given by $${\mathcal{A}}= \bigcup_{k =1}^{N!}\ell_{k} \label{eqn:AN},$$ where, for each $k \in \{1,2,\dots,N!\}$, $\ell_k := \{\tau  :\tau =  {\widetilde}\tau^k + {\mathbf 1}s \in P_N, s \in {\mathbb{R}}\}.$

Lyapunov Stability {#sec:lyapunov}
------------------

Lyapunov theory for hybrid systems is employed to show that the set of points ${\mathcal{A}}$ is asymptotically stable. Our candidate Lyapunov-like function, which is defined below and uses the distance function, is built by observing that there exist points where the distance to ${\mathcal{A}}$ may increase during flows. This is due to the sets $\ell_{k}$ being a subset $P_{N}$. To avoid this issue, we define $${\widetilde}{\mathcal{A}}= \bigcup_{k=1}^{N!} {\widetilde}\ell_k \supset {\mathcal{A}}$$ where ${\widetilde}\ell_k$ is the extension of $\ell_k$ given by $$\begin{aligned}
{\widetilde}\ell_k = \left\{\tau \in {\mathbb{R}}^{N} : \tau = {\widetilde}\tau^k + {\mathbf 1}s, s \in {\mathbb{R}}\right\}. \label{eqn:wtellk}\end{aligned}$$ Then, with this extended version of ${\mathcal{A}}$, the proposed candidate Lyapunov-like function for asymptotic stability of ${\mathcal{A}}$ for ${{{\mathcal}H}}_N$ is given by the locally Lipschitz function $$\label{eqn:lyapunov}
V(\tau) = \min\{|\tau|_{{\widetilde}\ell_1},|\tau|_{{\widetilde}\ell_2}, \ldots ,|\tau|_{{\widetilde}\ell_k}, \ldots ,|\tau|_{{\widetilde}\ell_N!}\} \quad \forall \ \tau \in P_N \setminus {\mathcal{X}}$$ where, for some $k$, $|\tau|_{{\widetilde}\ell_k}$ is the distance between the point $\tau$ and the set ${\widetilde}\ell_k$.[^6] The following theorem establishes asymptotic stability of ${\mathcal{A}}$ for ${{{\mathcal}H}}_N$. We show that the change in $V$ during flows is zero and that at jumps we have a strict decrease of $V$; namely, $V(G(\tau)) - V(\tau) = -|{\varepsilon}| V(\tau)$. A key step in the proof is in using [@teel2012hybrid Theorem 8.2] on a restricted version of ${{{\mathcal}H}}_N$.

\[thm:stability\] For every $N \in {\mathbb{N}}, N > 1$, ${\bar{\tau}}> 0,\omega > 0$, and $\varepsilon \in (-1,0)$, the hybrid system ${{{\mathcal}H}}_N$ is such that the compact set ${\mathcal{A}}$ is

Let the set ${\mathcal{X}}_{v}$ define the $v$-inflation of ${\mathcal{X}}$ (defined in Lemma \[lem:defX\]), that is, the open set[^7] ${\mathcal{X}}_{v} := \{\tau \in {\mathbb{R}}^N : |\tau|_{{\mathcal{X}}} < v\}$, where $v \in (0,v^{*})$ and $v^{*} = \min_{x\in{\mathcal{X}}, y\in{\widetilde}{\mathcal{A}}} |x-y|$. Given any $v \in (0,v^{*})$, we now consider a restricted hybrid system ${\widetilde}{{{\mathcal}H}}_{N} = (f,{\widetilde}{C},G,{\widetilde}{D})$, where ${\widetilde}{C} := C \setminus {\mathcal{X}}_{v}$ and ${\widetilde}{D} := D \setminus {\mathcal{X}}_{v}$, which are closed. We establish that ${\widetilde}{\mathcal{A}}$ is an asymptotically stable set for ${\widetilde}{{{\mathcal}H}}_N$.

Note that the continuous function $V$, given by , is defined as the minimum distance from $\tau$ to ${\widetilde}{\mathcal{A}}$, where ${\widetilde}{\mathcal{A}}$ is the union of $N!$ sets ${\widetilde}\ell_{k}$ in . To determine the change of $V$ during flows[^8], we consider the relationship between the flow map and the sets ${\widetilde}\ell_{k}$. The inner product between a vector pointing in the direction of the set ${\widetilde}\ell_{k}$ and the flow map on ${\widetilde}{C}$ satisfies $${\mathbf 1}^{\top}f(\tau) = {\mathbf 1}^{\top}(\omega{\mathbf 1}) = \omega N=|{\mathbf 1}||\omega {\mathbf 1}| = |{\mathbf 1}||f(\tau)|\cos \theta$$, which is only true if $\theta$ is zero. Therefore, the direction of the flow map and of the vector defining ${\widetilde}\ell_{k}$ are parallel, implying that the distance to the set ${\widetilde}{\mathcal{A}}$ is constant during flows.

The change in $V$ during jumps is given by $V(G(\tau)) - V(\tau)$ for $\tau \in {\widetilde}D \setminus {\widetilde}{\mathcal{A}}$. Due to the fact that we can rearrange the components of $\tau \in P_N \setminus {\mathcal{X}}$, without loss of generality, we consider a single jump condition, namely, we consider $\tau$ such that ${\bar{\tau}}= \tau_{1} > \tau_{2} > \ldots > \tau_{N-1} > \tau_{N}$. Using the formulation in Section \[sec:AN\] and Lemma \[eqn:tauSsolution\], the elements of the vector ${\widetilde}\tau^{k}$ associated with ${\widetilde}\ell_{k}$ for this case of $\tau$ are given by ${\widetilde}\tau_{i}^{k} = \frac{\sum_{p=0}^{N-i}({\varepsilon}+1)^p}{\sum_{p=0}^{N-1}({\varepsilon}+1)^p}{\bar{\tau}}$, which by Lemma \[lem:consum1\] is equal to $\frac{({\varepsilon}+1)^{N - i +1} - 1}{({\varepsilon}+1)^{N} - 1}{\bar{\tau}}$. After the jump, $G(\tau)$ is single valued and is such that its elements are ordered as follows: $g_{2}(\tau) > g_{3}(\tau) > \ldots > g_{N}(\tau) > g_{1}(\tau) = 0.$ Specifically, the jump map is $G(\tau) = [0,(1+{\varepsilon})\tau_{2},\ldots, (1+{\varepsilon})\tau_{N}]^{\top}$. Then, the formulation in Section \[sec:AN\] and [Lemma \[eqn:tauSsolution\]]{} leads to a case of ${\widetilde}\tau^{k}$ denoted as ${\widetilde}\tau^{k'}$. By [Lemma \[lem:consum1\]]{}, the elements of the vector ${\widetilde}\tau^{k'}$ are given by ${\widetilde}\tau^{k'}_{1} = \frac{{\varepsilon}}{({\varepsilon}+1)^{N} - 1}{\bar{\tau}}$ and ${\widetilde}\tau^{k'}_{i} = \frac{({\varepsilon}+1)^{N - i +2} - 1}{({\varepsilon}+1)^{N} - 1}{\bar{\tau}}$ for $i > 1$. Due to the ordering of $\tau$ and $G(\tau)$, ${\widetilde}\tau^{k'}$ is a one-element shifted (to the right) version of ${\widetilde}\tau^{k}$.

From the definition of ${\widetilde}\tau^{k}$ above, $V$ at $\tau$ reduces to $$V(\tau) = |\tau|_{{\widetilde}\ell_{k}} = \left|({\widetilde}\tau^{k} - \tau) - \frac{1}{N}(({\widetilde}\tau^{k} - \tau)^{\top}{\mathbf 1}){\mathbf 1}\right|$$ for some $k$. Note that $$({\widetilde}\tau^{k} - \tau)^{\top}{\mathbf 1}= \sum_{i = 1}^{N} {\widetilde}\tau^{k}_{i} - \sum_{i=1}^{N} \tau_{i}$$ reduces to $\sum_{i = 2}^{N} {\widetilde}\tau^{k}_{i} - \sum_{i=2}^{N} \tau_{i}$ since $\tau_{1} = {\widetilde}\tau_{1}^{k} = {\bar{\tau}}$. Using [Lemmas \[lem:consum1\] and \[lem:consum2\]]{}, it follows that $$\sum_{i = 2}^{N} {\widetilde}\tau^{k}_{i}= \frac{\sum_{i=2}^N \sum_{p=0}^{N-i}({\varepsilon}+1)^p}{\sum_{p=0}^{N-1}({\varepsilon}+1)^p}{\bar{\tau}}= \frac{(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}.$$ Then, the first element of the vector inside the norm in the expression of $V(\tau)$ is given as $$\begin{aligned}
&({\widetilde}\tau^{k}_{1} - \tau_{1}) - \frac{1}{N}\left(\frac{(({\varepsilon}+1)^{N} - 1 )- N{\varepsilon}}{{\varepsilon}(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}- \sum_{i=2}^{N} \tau_{i}\right) 
\\ & \qquad\qquad\qquad\qquad= -\frac{(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}N(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}+ \frac{1}{N}\sum_{i=2}^{N} \tau_{i},\end{aligned}$$ while the elements with $m \in \{2,3,\ldots,N\}$ are given by

After the jump at $\tau$, since $G(\tau)$ is single valued, $V(G(\tau))$ is given by $$|G(\tau)|_{{\widetilde}\ell_{k'}} = \left|({\widetilde}\tau^{k'} - G(\tau)) - \frac{1}{N}(({\widetilde}\tau^{k'} - G(\tau))^{\top}{\mathbf 1}){\mathbf 1}\right|.$$ Note that $({\widetilde}\tau^{k'} - G(\tau))^{\top}{\mathbf 1}= \sum_{i = 1}^{N} {\widetilde}\tau^{k'}_{i} - \sum_{i = 1}^{N} g_{i}(\tau)$ reduces to $\sum_{i = 1}^{N} {\widetilde}\tau^{k'}_{i} - \sum_{i = 2}^{N} (1+ {\varepsilon})\tau_i$, since $g_{1}(\tau) = 0$ and $g_{i}(\tau) = (1+{\varepsilon})\tau_{i}$ for $i > 1$. Using [Lemmas \[lem:consum1\] and \[lem:consum2\]]{}, it follows that $$\begin{aligned}
\sum_{i = 1}^{N} {\widetilde}\tau^{k'}_{i} &= \frac{\sum_{i=1}^N \sum_{p=0}^{N-i}({\varepsilon}+1)^p}{\sum_{p=0}^{N-1}({\varepsilon}+1)^p}{\bar{\tau}}\\ &= \frac{({\varepsilon}+1)(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}\end{aligned}$$ which leads to $$({\widetilde}\tau^{k'} - G(\tau))^{\top}{\mathbf 1}= \frac{({\varepsilon}+1)(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}- \sum_{i = 2}^{N} (1+{\varepsilon})\tau_{i}.$$ The first element inside the norm in $V(G(\tau))$ is given by $$\begin{aligned}
&({\widetilde}\tau_{1}^{k'} - g_1(\tau)) - \frac{1}{N}\left(\frac{({\varepsilon}+1)(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}\right. \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \left. - \sum_{i = 2}^{N} (1+{\varepsilon})\tau_{i}\right) \\
& = \frac{{\varepsilon}}{({\varepsilon}+1)^{N} - 1}{\bar{\tau}}- \frac{({\varepsilon}+1)(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}N(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}\\ 
& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad + \frac{1}{N}\sum_{i = 2}^{N} (1+{\varepsilon})\tau_{i} \\
 &= (1+{\varepsilon})\left(-\frac{(({\varepsilon}+1)^{N} - 1) - N{\varepsilon}}{{\varepsilon}N(({\varepsilon}+1)^{N} - 1)}{\bar{\tau}}+ \frac{1}{N}\sum_{i = 2}^{N} \tau_{i}\right).\end{aligned}$$ For each element $m > 1$, it follows that

Combining the expressions for each of the elements inside the norm of $V(G(\tau))$, it follows that $V(G(\tau)) = (1+{\varepsilon}) V(\tau)$.

Then, the change during jumps is given by $V(G(\tau)) - V(\tau) = {\varepsilon}V(\tau)$ where ${\varepsilon}\in (-1,0)$. With the property of $V$ during flows established above, the change of $V$ along solutions is bounded during flows and jumps by the nonpositive functions $u_{{\widetilde}{C}}$ and $u_{{\widetilde}{D}}$, respectively, defined as follows: $u_{{\widetilde}{C}} (z)= 0$ for each $z \in {\widetilde}{C}$ and $u_{{\widetilde}{C}} (z)= -\infty$ otherwise; $u_{{\widetilde}{D}}(z) = {\varepsilon}V(z)$ for each $z \in {\widetilde}{D}$ and $u_{{\widetilde}{D}}(z) = -\infty$ otherwise. Using Lemma \[lem:BasicConds\], the fact that ${\widetilde}{C}$ and ${\widetilde}{D}$ are closed, and the fact that every maximal solution to ${\widetilde}{{{\mathcal}H}}$ is bounded and complete, by [@teel2012hybrid Theorem 8.2], every maximal solution to ${\widetilde}{{{\mathcal}H}}_N$ approaches the largest weakly invariant subset of $L_{V}(r') \cap {\widetilde}{C} \cap [L_{u_{{\widetilde}{C}}}(0) \cup (L_{u_{{\widetilde}{D}}}(0) \cap G(L_{u_{{\widetilde}{C}}}(0)))] = L_{V}(r') \cap {\widetilde}{C}$ for $r' \in V({\widetilde}{C})$. Since every maximal solution jumps an infinite number of times, the largest invariant set is given for $r' = 0$ due to the fact that $V(G(\tau)) - V(\tau) = {\varepsilon}V(\tau) < 0$ if $r' > 0$. Then, the largest invariant set is given by $L_{V}(0) \cap {\widetilde}{C} = {\widetilde}{\mathcal{A}}\cap {\widetilde}{C}$ which is identically equal to ${\mathcal{A}}$. Hence, the set ${\mathcal{A}}$ is attractive. Stability is guaranteed from the fact that $V$ is nonincreasing during flows and strictly decreasing during jumps. Then, the set ${\widetilde}{\mathcal{A}}$ is asymptotically stable for the hybrid system ${\widetilde}{{{\mathcal}H}}_{N}$. We have that ${\mathcal{A}}$ is (strongly) forward invariant and from Theorem 3.4 we know that ${\mathcal{A}}$ is uniformly attractive from a neighborhood of itself. Then by Proposition 7.5 in [@teel2012hybrid], it follows that ${\mathcal{A}}$ is asymptotically stable.

Note that the set of solutions to ${\widetilde}{{{\mathcal}H}}_{N}$ coincides with the set of solutions to ${{{\mathcal}H}}_{N}$ from $P_{N} \setminus {\mathcal{X}}_{v}$. Therefore, the set ${\mathcal{A}}$ is asymptotically stable for ${{{\mathcal}H}}_N$ with basin of attraction ${{\mathcal}B}_{{\mathcal{A}}} = P_{N} \setminus {\mathcal{X}}_{v}$. [[ Since $v$ is arbitrary, it follows that the basin of attraction is equal to $P_{N} \stopmodif \setminus {\mathcal{X}}$.]{}]{}

Note that the jump map $G$, at points $\tau \in {\mathcal{X}}$, is set valued by definition of $g_i$ in . From these points there exist solutions to ${{{\mathcal}H}}_N$ that jump out of ${\mathcal{X}}$. In fact, consider the case $\tau \in {\mathcal{X}}$. We have that $\tau_{i} = \tau_r$ for some $i,r \in I$. Then, after the jump it follows that $g_{i}(\tau) \in \{0, (1+{\varepsilon}){\bar{\tau}}\}$ and $g_{r}(\tau) \in \{0, (1+{\varepsilon}){\bar{\tau}}\},$ and there exist $g_{i}$ and $g_{r}$ such that $g_{i} = g_{r}$ or $g_{i} \neq g_{r}$. Since for every point in ${\mathcal{X}}$ there exists a solution that converges to ${\mathcal{A}}$ and also a solution that stays in ${\mathcal{X}}$, ${\mathcal{X}}$ is weakly forward invariant.[^9]

Characterization of Time of Convergence {#sec:timetoconverge}
---------------------------------------

In this section, we characterize the time to converge to a neighborhood of ${\mathcal{A}}$. The proposed (upper bound) of the time to converge depends on the initial distance to the set ${\widetilde}{\mathcal{A}}$ and the parameters of the hybrid system $(\varepsilon,{\bar{\tau}})$.

\[thm:timetoconverge\] For every $N \in {\mathbb{N}}$, $N > 1$, and every $c_1, c_2$ such that $\overline{c} > c_2 > c_1 > 0$ with $\overline{c} = \max_{x \in {\mathcal{X}}} |x|_{\widetilde{\mathcal{A}}}$, every maximal solution to ${{{\mathcal}H}}_N$ with initial condition $\tau(0,0) \in (P_{N} \setminus {\mathcal{X}}) \cap {\widetilde}L_{V}(c_2)$ is such that

Let $\tau_0 = \tau(0,0)$ and pick a maximal solution $\tau$ to ${{{\mathcal}H}}_{N}$ from $\tau_{0}$. At every jump time $(t_{j},j) \in {\mathop{\rm dom}\nolimits}\tau$, define $\bar{g}_1 = \tau(t_1,1)$, $\bar{g}_2 = \tau(t_2,2), \ldots, \bar{g}_J = \tau(t_J,J)$, for some $J \in {\mathbb{N}}$. From Theorem \[thm:stability\], we have that there is no change in the Lyapunov function during flows. Furthermore, we have that for each $\tau \in D \setminus {\mathcal{A}}$ the difference $V(G(\tau)) - V(\tau) = \varepsilon V(\tau)$ with $\varepsilon \in (-1,0)$. Since, for every $j,$ $\tau(t_{j},j) \in D$, we have which implies At the next jump, we have Proceeding in this way, after $J$ jumps we have From $V(\bar g_{J}) = (1+\varepsilon)^{J}V(\tau_{0})$, we want to find $J$ so that $V(\bar g_{J}) \leq c_{1}$ when $V(\tau_{0}) \leq c_{2}$. Considering the worst cast for $V(\tau_{0})$, we want $(1+\varepsilon)^{J}c_{2} \leq c_{1}$, which implies $\frac{c_{2}}{c_{1}} \leq \left(\frac{1}{1+\varepsilon} \right)^{J}$, and therefore $J = \left\lceil\frac{\log \frac{c_{2}}{c_{1}}}{ \log \frac{1}{1+\varepsilon}}\right\rceil > 0.$ For each $j$, the time between jumps satisfies $t_1 - t_0 \leq \frac{{\bar{\tau}}}{\omega}, 
t_2 - t_1 \leq \frac{{\bar{\tau}}}{\omega}, 
\ldots, 
t_j - t_{j-1} \leq \frac{{\bar{\tau}}}{\omega}.$ Then, we have that after $J$ jumps, $
\sum_{j=1}^{J} t_{j} - t_{j-1} \leq J\frac{{\bar{\tau}}}{\omega}.
$ With $t_{0} = 0$, the expression reduces to $t_{J} \leq J\frac{{\bar{\tau}}}{\omega}  = \left\lceil\frac{\log \frac{c_{2}}{c_{1}}}{ \log \frac{1}{1+\varepsilon}}\right\rceil\frac{{\bar{\tau}}}{\omega}.$ Then, after $t+j \geq t_{J} +J$, the solution is at least $c_{1}$ close to the set ${\widetilde}{\mathcal{A}}$. Defining $M = t_{J} + J$ we then have

Figure \[fig:timetoconverge\] shows the time to converge (divided by $\frac{{\bar{\tau}}}{\omega}+1$) versus $\varepsilon$ with constant $c_{2} = 0.99{\bar{\tau}}$ and varying values of $c_{1}$. As the figure indicates, the time to converge decreases as $|{\varepsilon}|$ increases, which confirms the intuition that the larger the jump the faster oscillators desynchronize.

Robustness Analysis {#sec:robustness}
-------------------

[[ Lemma \[lem:BasicConds\] establishes that the hybrid model of $N$ impulse-coupled oscillators satisfies the hybrid basic conditions. In light of this property, the asymptotic stability property of ${\mathcal{A}}$ for ${{{\mathcal}H}}_N$ is preserved under certain perturbations; i.e., asymptotic stability is robust [@teel2012hybrid]. In the next sections, we consider a perturbed version of ${{{\mathcal}H}}_N$ and present robust stability results. In particular, we consider generic perturbations to ${{{\mathcal}H}}_N$, and two different cases of perturbations only on the timer rates to allow for heterogeneous timers.]{}]{}

### Robustness to Generic Perturbations

We start by revisiting the definition of perturbed hybrid systems in [@teel2012hybrid].

.1cm

Given a hybrid system ${{{\mathcal}H}}$ and a function $\rho: {\mathbb{R}}^N \to {\mathbb{R}}_{\geq 0}$, the $\rho$-perturbation of ${{{\mathcal}H}}$, denoted ${{{\mathcal}H}}_\rho$, is the hybrid system $$\left\{\begin{array}{cc}x \in C_\rho & \quad \dot{x} \in F_\rho(x) \\
x \in D_\rho & \quad x^+ \in G_\rho(x) \\
 \end{array} \right.$$ where

.1cm

Using this definition, we can deduce a generic perturbed hybrid system modeling $N$ impulse-coupled oscillators. Then, for the hybrid system ${{{\mathcal}H}}_N$, we denote ${{{\mathcal}H}}_{N,\rho}$ as the $\rho$-perturbation of ${{{\mathcal}H}}_N$. Given the perturbation function $\rho : {\mathbb{R}}^N \to {\mathbb{R}}_{\geq 0}$, the perturbed flow map is given by where the perturbed flow set $C_\rho$ is given by For example, if $N = 2$ and $\rho(\tau) = \bar\rho > 0$ for all $\tau \in {\mathbb{R}}^N$, which would correspond to constant perturbations on the lower value and threshold, then $C_\rho = C + \rho {\mathbb{B}}$. The perturbed jump map and jump set are defined as where $g_{i,\rho}$ is the $i$-th component of $G_\rho$. The following result establishes that the hybrid system ${{{\mathcal}H}}_N$ is robust to small perturbations.

(robustness of asymptotic stability)\[thm:robustofAS\] If $\rho : {\mathbb{R}}^N \to {\mathbb{R}}_{\geq 0}$ is continuous and positive on ${\mathbb{R}}^{N} \setminus {\mathcal{A}}$, then ${\mathcal{A}}$ is semiglobally practically robustly ${\mathcal}{KL}$ asymptotically stable with basin of attraction $B_{{\mathcal{A}}} = P_N \setminus {\mathcal{X}}$, i.e., for every compact set $K \subset B_{{\mathcal{A}}}$ and every $\alpha > 0$, there exists $\delta \in (0,1)$ such that every maximal solution $\tau$ to ${{{\mathcal}H}}_{N,\delta\rho}$ from $K$ satisfies $|\tau(t,j)|_{{\mathcal{A}}} \leq \beta(|\tau(0,0)|_{{\mathcal{A}}},t+j) + \alpha$ for all $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$.

From Lemma \[lem:BasicConds\], the hybrid system ${{{\mathcal}H}}_N$ satisfies the hybrid basic conditions. Therefore, by [@teel2012hybrid Theorem 6.8] ${{{\mathcal}H}}_N$ is nominally well-posed and, moreover, by [@teel2012hybrid Proposition 6.28] is well-posed. From the proof of Theorem \[thm:stability\], we know that the set ${\mathcal{A}}$ is an asymptotically stable compact set for the hybrid system ${{{\mathcal}H}}_N$ with basin of attraction $B_{{\mathcal{A}}}$. Since by Lemma \[lem:solutions2HS\], every maximal solution is complete, then [@teel2012hybrid Theorem 7.20] implies that ${\mathcal{A}}$ is semiglobally practically robustly ${\mathcal}{KL}$ asymptotically stable.

Section \[sec:jumpperturbs1\] showcases of ${{{\mathcal}H}}_{N}$ with $\rho$-perturbations on the jump map.

### Robustness to Heterogeneous Timer Rates

We consider the case when the continuous dynamic rates are perturbed in the form of for a given solution $\tau$. For example, consider the perturbation of the flow map given by $$\begin{aligned}
f(\tau) = \omega {\mathbf 1}+ \Delta\omega \label{eqn:fdelta}\end{aligned}$$ where $\Delta\omega \in {\mathbb{R}}^n$ is a constant defining a perturbation from the natural frequencies of the impulse-coupled oscillators. Then for some $k$, during flows, along a solution $\tau$ such that over $[t_{j},t_{j+1}]\times\{j\}$ satisfies $V(\tau(t,j)) = |\tau(t,j)|_{{\widetilde}\ell_{k}}$, it follows that $c$ reduces to $c(t,j) = \left(\frac{r_{\ell_{k}}^{\top}(\tau(t,j))(\frac{1}{N}{\bf \underline{1} - \bf I})}{|\tau(t,j)|_{\ell_{k}}}\right) \Delta\omega.$[^10] Furthermore, the norm of the hybrid arc $c$ can be bounded by a constant $\bar{c}$ given by $$\begin{aligned}
\bar{c} = \left|\left(\frac{1}{N}\underline{\bf 1} - {\bf I}\right)\Delta\omega\right| \label{eqn:ctbound}.\end{aligned}$$ Building from this example, the following result provides properties of the distance to ${\widetilde}{\mathcal{A}}$ from solutions $\tau$ to ${{{\mathcal}H}}_N$ under generic perturbations on $f$ (not necessarily as in ).

\[thm:vanishingC\] Suppose that the perturbation on the flow map of ${{{\mathcal}H}}_N$ is such that a perturbed solution $\tau$ satisfies, for each $j$ such that $\{t : (t,j) \in {\mathop{\rm dom}\nolimits}\tau\}$ has more than one point, $\frac{d}{dt} |\tau(t,j)|_{{\widetilde}{\mathcal{A}}} = c(t,j)$ for all $t \in \{t : (t,j) \in {\mathop{\rm dom}\nolimits}\tau\}$ [[ and $\tau(t,j) \in P_N \setminus {\mathcal{X}}$ for all $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$]{}]{}, for some hybrid arc $c$ with ${\mathop{\rm dom}\nolimits}c = {\mathop{\rm dom}\nolimits}\tau$. Then, the following hold:

-   The asymptotic value of $|\tau(t,j)|_{{\widetilde}{\mathcal{A}}}$ satisfies $$\begin{aligned}
    \lim_{t+j \to \infty}|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq \lim_{t+j \to \infty} \sum_{i = 0}^{j} (1+{\varepsilon})^{j-i} \int_{t_i}^{t_{i+1}} c(t,j) dt\end{aligned}$$

-   If there exists $\bar{c} > 0$ such that $|c(t,j)| \leq \bar{c}$ for each $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$ then $$\begin{aligned}
    \lim_{t+j \to \infty}|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq \frac{\bar{c}{\bar{\tau}}}{|{\varepsilon}|\omega}. \label{thmeqn:distupperbound}\end{aligned}$$

-   If ${\widetilde}{j} : {\mathbb{R}}_{\geq 0} \to {\mathbb{N}}$ is a function that [[ chooses the appropriate minimum $j$ such that $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$]{}]{} for each time $t$ and $t \mapsto c(t,{\widetilde}{j}(t))$ is absolutely integrable, i.e., $\exists B$ such that then $$\begin{aligned}
    \lim_{t + j \to \infty} |\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq \frac{B}{{\varepsilon}}. \label{eqn:limboundbeps}\end{aligned}$$

Consider a maximal solution $\tau$ to ${{{\mathcal}H}}_N$ with initial condition $\tau(0,0) \in P_N \setminus {\mathcal{X}}$. This proof uses the function $V$ from the proof of Theorem \[thm:stability\]. With $V$ equal to the distance from $\tau$ to the set ${\widetilde}{\mathcal{A}}$, then, for each $\tau \in D \setminus {\mathcal{X}}$, we have that $V(G(\tau)) - V(\tau) = {\varepsilon}V(\tau)$. Using the fact that $V(\tau) = |\tau|_{{\widetilde}{\mathcal{A}}}$ and the fact that, $G$ along the solution is single valued, it follows that $|\tau|_{{\widetilde}{\mathcal{A}}}$ after a jump can be equivalently written as By assumption, in between jumps, the distance to the set ${\widetilde}{\mathcal{A}}$ is such that $\frac{d}{dt}|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} = c(t,j)$, which implies that at $t_{j+1}$ the distance to the desynchronization set is given by $$|\tau(t_{j+1},j)|_{{\widetilde}{\mathcal{A}}} = \int_{t_{j}}^{t_{j+1}} c(s,j) ds + |\tau(t_{j},j)|_{{\widetilde}{\mathcal{A}}}.$$ It follows that Then, proceeding in this way, we obtain $$\begin{aligned}
|\tau(t_{j},j)|_{{\widetilde}{\mathcal{A}}} &= (1+{\varepsilon})^{j}|\tau(0,0)|_{{\widetilde}{\mathcal{A}}} \\ & \qquad  \qquad \qquad +\sum_{i = 0}^{j-1}(1+{\varepsilon})^{j-i}\int_{t_{i}}^{t_{i+1}} c(s,i)ds. \end{aligned}$$ For the case of generic $t_{j+1} \geq t \geq t_j$, we have that $$\begin{aligned}
|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} = (1+{\varepsilon})^{j}|\tau(0,0)|_{{\widetilde}{\mathcal{A}}} + \sum_{i = 0}^{j}(1+{\varepsilon})^{j-i}\int_{t_{i}}^t c(s,i)ds.\end{aligned}$$ Since, we know that as either $t$ or $j$ goes to infinity, $j$ or $t$ go to infinity as well, respectively. The expression reduces to If $c(t,j) \leq \bar{c}$, it follows that

Lastly, since this hybrid system has the property that for any maximal solution $\tau$ with $(t,j) \in {\mathop{\rm dom}\nolimits}\tau$, if $t$ approaches $\infty$ then the parameter $j$ also approaches $\infty$, the expression given by $\lim_{t + j \to \infty} |\tau(t,j)|_{{\widetilde}{\mathcal{A}}}$ can be simplified. To do this, we know that the series $\sum_{i=0}^j (1+{\varepsilon})^{j -i} = \frac{(1+{\varepsilon})^{j+1}-1}{{\varepsilon}}$ approaches $\frac{1}{|{\varepsilon}|}$ as $j \to \infty$. Since $1+{\varepsilon}>0$ for ${\varepsilon}\in (-1,0)$, the series is absolutely convergent and its partial sum $s_j = \sum_{i=0}^j (1+{\varepsilon})^{j -i}$ is such that $\{s_j\}^\infty_{j=m}$ is a nondecreasing sequence (for each $m$). This implies that $s_j \leq 1/|{\varepsilon}|$ for all $j$ and for each $m$. Then, it follows that $(1+{\varepsilon})^{j-i} \leq \frac{1}{|{\varepsilon}|}$ for every $j,i \in {\mathbb{N}}$. Since the expression is a function of $j$ only and, for complete solutions, $t$ is such that as $t \to \infty$, then $j \to \infty$, we obtain

Numerical Analysis {#sec:numerics}
==================

This section presents numerical results obtained from simulating ${{{\mathcal}H}}_N$. First, we present results for the nominal case of ${{{\mathcal}H}}_N$ given by . Then, we present results for ${{{\mathcal}H}}_N$ under different types of perturbations. The Hybrid Equations (HyEQ) Toolbox in [@Sanfelice.ea.13.HSCC] was used to compute the trajectories.

Nominal Case
------------

The possible solutions to the hybrid system ${{{\mathcal}H}}_N$ fall into four categories: always desynchronized, asymptotically desynchronized, never desynchronized, and initially synchronized. The parameters used in these simulations are ${\bar{\tau}}= 1$ and $\varepsilon = -0.2$.

A solution of ${{{\mathcal}H}}_N$ that starts in $P_N \setminus ({\mathcal{X}}\cup {\mathcal{A}})$ asymptotically converges to ${\mathcal{A}}$, as Theorem \[thm:timetoconverge\] indicates. show solutions to both ${{{\mathcal}H}}_{2}$ and ${{{\mathcal}H}}_{3}$ converging to their respective desynchronization sets.

 \
\[\]\[\]\[.7\][$\tau_1,\tau_{2}$]{} \[\]\[\]\[.7\][$t$ \[seconds\]]{} \[\]\[\]\[.8\]

For ${{{\mathcal}H}}_2$, if $\tau(0,0) = [0,0.1]^\top$, then the initial sublevel set is ${\widetilde}L_{V}(c_{2})$ with $c_{2} = 0.24$. Using Theorem \[thm:timetoconverge\], the time to converge to the sublevel set ${\widetilde}L_{V}(c_{1})$ with $c_{1} = 0.1$ leads to $M = 7.84$. Figure \[fig:H2notinX2\] shows a solution to the system for 10 seconds of flow time. From the figure, it can be seen that $V(\tau(t,j)) \approx 0.1$ at $(t,j) = (3,4)$. Then, the property guaranteed by Theorem \[thm:timetoconverge\], namely, $V(\tau(t,j)) \leq c_{1}$ for each $(t,j)$ such that $t + j \geq M$, is satisfied. Figure \[fig:H3notinX3\], shows a solution and the distance of this solution to ${\mathcal{A}}$. Notice that the initial sub level set is ${\widetilde}{L}_{V}(c_{2})$ with $c_{2} =0.32$. From Theorem \[thm:timetoconverge\] it follows that the time to converge to ${\widetilde}{L}_{V}(c_{1})$ with $c_{1} = 0.1$ is given by $M =10.14$, which is actually already satisfied at $(t,j) = (2.2,4)$. show solutions to ${{{\mathcal}H}}_N$ that asymptotically desynchronize for $N \in \{7,10\}$.

Perturbed Case
--------------

In this section, we present numerical results to validate the statements in Section \[sec:robustness\].

### Simulations of ${{{\mathcal}H}}_N$ with perturbed jumps {#sec:jumpperturbs1}

$\bullet$ [**Perturbation of the threshold in the jump set:**]{} We replace the jump set $D$ by $D_{\rho} := \{\tau : \exists i \in I \ s.t. \ \tau_i = {\bar{\tau}}+ \rho_i\}$ where $\rho_i \in [0, \bar{\rho}_i]$, $\bar{\rho}_i > 0$ for each $i \in I$. To avoid maximal solutions that are not complete, the flow set $C$ is replaced by $C_\rho := [0,{\bar{\tau}}+\rho_1]\times[0,{\bar{\tau}}+\rho_2] \times \ldots \times [0,{\bar{\tau}}+\rho_N]$. Furthermore, the components of the jump map are also replaced by $$g_{\rho_i} (\tau) = \left\{ \begin{array}{ll} 0 \qquad \qquad \ \ \ \ &\mbox{if } \tau_{i} = {\bar{\tau}}+\rho_i, \tau_r < {\bar{\tau}}+\rho_j \ \ \forall j \in I\setminus \{i\} \\ 
\{0 , \tau_{i}(1+\varepsilon) \} \ &\mbox{if } \tau_{i} = {\bar{\tau}}+\rho_i \ \exists j \in I \setminus \{i\} \  \mbox{s.t.} \  \tau_r = {\bar{\tau}}+\rho_j \\ 
(1+\varepsilon)\tau_{i} \qquad \  &\mbox{if } \tau_i < {\bar{\tau}}+\rho_i \ \exists j \in I \setminus \{i\} \  \mbox{s.t.} \  \tau_r = {\bar{\tau}}+\rho_j\end{array} \right. \label{eqn:gi_perturbed} .$$

$\bullet$ [**Perturbations on the reset component of the jump map:**]{} Under the effect of the perturbations considered in this case, instead of reseting $\tau_i$ to zero, the perturbed jump resets $\tau_i$ to a value $\rho_i \in {\mathbb{R}}_{\geq 0}$, for each $i \in I$. The perturbed hybrid system has the following data: and where, for each $i \in I$, the perturbed jump map is given by $$g_i (\tau) = \left\{ \begin{array}{l} \rho_i \qquad \qquad \ \ \ \ \mbox{if } \tau_{i} = {\bar{\tau}}, \tau_r < {\bar{\tau}}\ \ \forall j \in I\setminus \{i\} \\ 
\{\rho_i , \tau_{i}(1+\varepsilon) \} \ \mbox{if } \tau_{i} = {\bar{\tau}}\ \exists j \in I \setminus \{i\} \  \mbox{s.t.} \  \tau_r = {\bar{\tau}}\\ 
(1+\varepsilon)\tau_{i} \qquad \  \mbox{if } \tau_i < {\bar{\tau}}\ \exists j \in I \setminus \{i\} \  \mbox{s.t.} \  \tau_r = {\bar{\tau}}\end{array} \right. .\label{eqn:gi_resetperturbation}$$

This case of perturbations exemplifies Theorem \[thm:robustofAS\] with $\rho$ affecting only the jump map of ${{{\mathcal}H}}_N$. Figures \[figs:ResetPerturb\_di\_equal\_di=0.2\] and \[figs:ResetPertub\_di\_notequal\_di\] show several simulations to this perturbation of ${{{\mathcal}H}}_N$. All of the simulations in this section use parameters $\omega = 1$, ${\bar{\tau}}= 3$, ${\varepsilon}= -0.3$, and $N = 2$.

The first case of the perturbed jump map $G_\rho$ considered is for $\rho_1 = \rho_2 = 0.02$. Figure \[ResetPerturb\_di\_equal\_t1t2plot\_di=0.2\] shows a solution to the perturbed ${{{\mathcal}H}}_2$ from the initial condition $\tau(0,0) = [2.4,2.3]^\top$ on the $(\tau_1,\tau_2)$-plane. Notice that for $\tau \in D$ such that $\tau_i = {\bar{\tau}}$ the jump map resets $\tau_i$ to $\rho_i$ (red dashed line) and not to $0$ as in the unperturbed case. The solution for this case approaches a region around ${\widetilde}{\mathcal{A}}$, as Theorem \[thm:robustofAS\] guarantees. Figure \[ResetPerturb\_di\_equal\_timeplot\_di=0.2\] shows the distance to the set ${\widetilde}{\mathcal{A}}$ over time for 10 solutions of the perturbed system ${{{\mathcal}H}}_2$ with initial conditions $\tau(0,0) \in P_2\setminus {\mathcal{X}}_{2}$. This figure shows that solutions approach a distance of about $0.12$ after 25 seconds.

Now, consider the case where $\rho_1 \neq \rho_2$. Figure \[figs:ResetPertub\_di\_notequal\_di\] shows the distance to ${\widetilde}{\mathcal{A}}$ for two sets of solutions with different values for $\rho_1$ and $\rho_2$. More specifically, Figure \[ResetPertub\_di\_NotEqual\_timeplot\_di=\[0.15,0.25\]\] shows the case of $\rho_1 = 0.15$ and $\rho_2 = 0.25$. For this case, it can be seen that the solutions converge after $\approx 28$ seconds of flow time and, after that time, satisfy $|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq 0.25$. Figure \[ResetPertub\_di\_NotEqual\_timeplot\_di=\[0.01,0.02\]\] shows the case of $\rho_1 = 0.02$ and $\rho_2 = 0.01$. For this case, this figure shows that, after $\approx 28$ seconds of flow time, the solutions satisfy $|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq 0.04$. These simulations validate Theorem \[thm:robustofAS\] with $\rho$ affecting only the jump map, verifying that the smaller the size of the perturbation the smaller the steady-state value of the distance to ${\widetilde}{\mathcal{A}}$.

the component $(1+{\varepsilon})\tau_i$ of the jump map is perturbed, namely, we use $\tau_i^+ = (1+{\varepsilon}) \tau_i + \rho_i(\tau_i)$, where $\rho_i : {\mathbb{R}}_{\geq 0} \to P_N \setminus {\mathcal{X}}$ is a continuous function. The perturbed jump map $G_{\rho}$ has components $g_{\rho i}$ that are given as $g_{i}$ in but with $\tau_{i}(1+{\varepsilon}) + \rho_{i}(\tau_{i})$ replacing $\tau_{i}(1+{\varepsilon})$.

Consider the case $\rho_i(\tau_{i}) = {\widetilde}\rho_i \tau_i$ with ${\widetilde}\rho_i \in (0,|{\varepsilon}|)$ and let ${\widetilde}{\varepsilon}_i = {\varepsilon}+ {\widetilde}\rho_i \in (-1,0)$. Then $\tau_i^+$ reduces to $\tau_i^+ = (1+{\widetilde}{\varepsilon}_i)\tau_i$ and the jump map $g_{\rho i}$ is given by with ${\widetilde}{\varepsilon}_{i}$ in place of ${\varepsilon}$. This type of perturbation is used to verify Theorem \[thm:robustofAS\] with $\rho$ affecting only the “bump” portion of the jump map. and \[figs:ResetBumpPerturb\_di\_notequal\_di\] show simulations to ${{{\mathcal}H}}_N$ with the parameters $\omega = 1$, ${\bar{\tau}}= 3$, ${\varepsilon}= -0.3$, and $N = 2$.

Consider the case of ${{{\mathcal}H}}_2$ with $G_\rho$ when ${\widetilde}\rho_1 = {\widetilde}\rho_2 = 0.1$, leading to ${\widetilde}{\varepsilon}_1 = {\widetilde}{\varepsilon}_2$ = 0.2. shows a solution on the $(\tau_1,\tau_2)$-plane for this case with initial condition $\tau(0,0) = [0.1,0.2]^\top$. Notice that the solution approaches a region around ${\mathcal{A}}$ (green line), as Theorem \[thm:robustofAS\] guarantees. Figure \[ResetBumpPertub\_di\_equal\_timeplot\_di=\[0.1,0.1\]\] shows the distance to the set ${\widetilde}{\mathcal{A}}$ over time for 10 solutions with initial conditions $\tau(0,0) \in C$. It shows that solutions approach a distance to ${\widetilde}{\mathcal{A}}$ of $\approx0.09$ after $\approx40$ seconds of flow time.

Next, we consider the case of $G_{\rho}$ with ${\widetilde}{\varepsilon}_1 \neq {\widetilde}{\varepsilon}_2$. Figure \[ResetBumpPertub\_di\_NotEqual\_timeplot\_di=\[0.15,0.1\]\] shows the distance to ${\widetilde}{\mathcal{A}}$ for 10 solutions with perturbations given by ${\widetilde}\rho_1 = 0.15$ and ${\widetilde}\rho_2 = 0.1$. For this case, the distance to ${\widetilde}{\mathcal{A}}$ satisfies $|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq 0.3$ after $\approx40$ seconds of flow time. Figure \[fig:ResetBumpPertub\_di\_NotEqual\_timeplot\_di=\[0.02,0.01\]\] shows simulation results with ${\widetilde}\rho_1 = 0.02$ and ${\widetilde}\rho_2 = 0.01$. Notice that the smaller the value of the perturbation is, the closer the solutions get to the set ${\widetilde}{\mathcal{A}}$. For this case, after $\approx30$ seconds of flow time, the distance to ${\widetilde}{\mathcal{A}}$ satisfies $|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq 0.06$. These simulations validate Theorem \[thm:robustofAS\] with $\rho$ affecting only the jump map, verifying that the smaller the size of the perturbation the smaller the steady-state value of the distance to ${\widetilde}{\mathcal{A}}$ would be.

### Perturbations on the Flow Map {#sec:flowpertrubs}

In this section, we consider a class of perturbations on the flow map. More precisely, consider the case when there exists a function $(t,j) \mapsto c(t,j)$ such that $c(t,j) \leq \bar{c}$ with $\bar{c}$ as in . Then, from Theorem \[thm:vanishingC\] with , we know that $$\begin{aligned}
\lim_{t+j \to \infty} |\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq \left|\frac{\bar{c}{\bar{\tau}}}{{\varepsilon}\omega}\right| \leq \left|\frac{\left|(\frac{1}{N}\underline{\bf 1} - {\bf I})\Delta\omega\right|{\bar{\tau}}}{{\varepsilon}\omega}\right|. \label{eqn:flowpert2}\end{aligned}$$

Figure  shows a simulation so as to verify this property. The parameters of this simulation are $N = 2$, $\omega = 1$, ${\varepsilon}= -0.3$, ${\bar{\tau}}= 4$, and $\Delta\omega = [0.120,0.134]^\top$. It follows from that $\overline{c} = 0.0105$. Then, from , it follows that $\lim_{t+j \to \infty}|\tau(t,j)|_{{\widetilde}{\mathcal{A}}} \leq 0.1047$. Specifically, Figure \[fig:c=constant\_tau1tau2plot\] shows a solution on the $(\tau_1,\tau_2)$-plane of the perturbed hybrid system ${{{\mathcal}H}}_2$ with initial condition $\tau(0,0) = [0,0.01]^\top$. This figure shows the solution (blue line) converging to a region around ${\widetilde}{\mathcal{A}}$ (between dash-dotted lines about ${\mathcal{A}}$ in green). Figure \[fig:c=constant\_disttimetraj\] shows the distance to the set ${\widetilde}{\mathcal{A}}$ of 10 solutions with initial conditions $\tau(0,0) \in C$ with a dashed line denoting the upper bound on the distance in . Notice that all solutions are within this bound after approximately 15 seconds of flow time and stay within this region afterwards.

Conclusion {#sec:conclusion}
==========

We have shown that desynchronization in a class of impulse-coupled oscillators is an asymptotically stable and robust property. These properties are established within a solid framework for modeling and analysis of hybrid systems, which is amenable for the study of synchronization and desynchronization in other impulse-coupled oscillators in the literature. The main difficulty in applying these tools lies on the construction of a Lyapunov-like quantity certifying asymptotic stability. As we show here, invariance principles can be exploited to relax the conditions that those functions have to satisfy, so as to characterize convergence, stability, and robustness in the class of systems under study. Future directions of research include the study of nonlinear reset maps, such as those capturing the phase-response curve of spiking neurons, as well as impulse-coupled oscillators connected via general graphs.

Appendix {#sec:appendix}
========

The following result derives the solution to $\Gamma\tau_{s} = b$ with $\Gamma$ given in and $b = {\bar{\tau}}{\mathbf 1}$ via Gaussian elimination.

\[eqn:tauSsolution\] For each ${\varepsilon}\in (-1,0)$, the solution $\tau_s$ to $\Gamma\tau_{s} = b$ with $\Gamma$ given in and $b = {\bar{\tau}}{\mathbf 1}$ is such that its elements, denoted as $\tau_s^k$ for each $k \in \{1,2,\ldots,N\}$, are given by $\tau_{s}^{k} = \frac{\sum_{i=0}^{N-k}({\varepsilon}+1)^i}{\sum_{i=0}^{N-1}({\varepsilon}+1)^i}{\bar{\tau}}$.

The $N \times N$ matrix in and the $N \times 1$ matrix $b = {\bar{\tau}}{\mathbf 1}$ leads to the augmented matrix $[\Gamma|b]$ given by To solve for $\tau_{s}^{k}$, we apply the Gauss-Jordan elimination technique to to remove the elements $-({\varepsilon}+1)$ above the diagonal. Starting from the $N$-th row to remove the $-({\varepsilon}+1)$ component in the $N-1$ row, and continuing up to the second row, gives Denoting the augmented matrix in as $[\Gamma'|b']$, with $\tau_{s}^{1} = {\bar{\tau}}$ and $\tau_{s}^{2} = \frac{\sum_{i=0}^{N-2}({\varepsilon}+1)^i}{\sum_{i=0}^{N-1}({\varepsilon}+1)^i}{\bar{\tau}}$, the solution for each element of $\tau_{s}^{k}$ with $k > 2$ can be derived from as $\Gamma'_{k,2}\tau^{2}_{s} + \tau^{k}_{s} = b'_{k}$ where $\Gamma'_{k,2}$ denotes the $(k,2)$ entry of $\Gamma'$. Noting that $\tau_{s}^{1}$ can be rewritten as $\tau_{s}^{1} = \frac{\sum_{i=0}^{N-1}({\varepsilon}+1)^i}{\sum_{i=0}^{N-1}({\varepsilon}+1)^i}{\bar{\tau}}$ leads to $\tau_{s}^{k} = \frac{\sum_{i=0}^{N-k}({\varepsilon}+1)^i}{\sum_{i=0}^{N-1}({\varepsilon}+1)^i}{\bar{\tau}}$.

\[lem:consum1\] For each $x \neq 1$, and $m,n \in {\mathbb{N}}$ such that $n-1 \geq m$, the finite sum $\sum_{i=m}^{n-1} x^i$ satisfies $
\sum_{i=m}^{n-1} x^i = \frac{x^n - x^m}{x-1}.
$

\[lem:consum2\] For each $x \neq 1$, [[ and each $m,N \in {\mathbb{N}}$ such that $N \geq m$]{}]{}, the finite sum $\sum_{n=m}^N\sum_{i=0}^{N-n}x^i$ satisfies

[^1]: Department of Computer Engineering, University of California, Santa Cruz, CA 95064. Email: [seaphill,ricardo@ucsc.edu]{}. This research has been partially supported by the National Science Foundation under CAREER Grant no. ECS-1150306 and by the Air Force Office of Scientific Research under Grant no. FA9550-12-1-0366.

[^2]: Cf. the model for synchronization in [@Mirollo.90.SIAMJAM.BiologicalOscillators] where $\varepsilon > 0$.

[^3]: In [@mauroy:037122], a more general flow map and a jump map incrementing $\tau_i$ by ${\varepsilon}> 0$ are considered.

[^4]: A set-valued mapping $G : {\mathbb{R}}^N{\rightrightarrows}{\mathbb{R}}^N$ is [*outer semicontinuous*]{} if its graph $\{(x,y): x \in {\mathbb{R}}^N, y \in G(x)\}$ is closed, see [@teel2012hybrid Lemma 5.10] and [@RockafellarWets98].

[^5]: Note that $G$ is single valued at each ${\widetilde}\tau^k \notin {\mathcal{X}}$.

[^6]: The set ${\widetilde}\ell_k$ can be described as a straight line in ${\mathbb{R}}^n$ passing through a point ${\widetilde}\tau^k$ and with slope ${\mathbf 1}$. Then, $|\tau|_{{\widetilde}\ell_k}$ can be written as the general point-to-line distance $|({\widetilde}\tau^k - \tau) - 1/N (({\widetilde}\tau^k - \tau)^\top{\mathbf 1}){\mathbf 1}|$.

[^7]: The set ${\mathcal{X}}_{v}$ is open since every point $\tau \in {\mathcal{X}}_{v}$ is an interior point of ${\mathcal{X}}_{v}$.

[^8]: Its derivative can be computed using Clarke’s generalized gradient [@Clarke90].

[^9]: For example, consider the case $N = 2$. If $\tau(0,0) = [{\bar{\tau}},{\bar{\tau}}]^\top \in D$, then there are nonunique solutions due to the jump map begin set valued. It follows that after the jump, each $\tau_i$ can be mapped to any point in $\{0,\tau_i(1+{\varepsilon})\}$, which leads to any of the following four options of the states $(\tau_1,\tau_2)$ after such a jump: $(0,0),(0,{\bar{\tau}}(1+{\varepsilon})),({\bar{\tau}}(1+{\varepsilon}),0)$ or $({\bar{\tau}}(1+{\varepsilon}),{\bar{\tau}}(1+{\varepsilon}))$. If the state is mapped to either $(0,0)$ or $({\bar{\tau}}(1+{\varepsilon}),{\bar{\tau}}(1+{\varepsilon}))$, then it remains in ${\mathcal{X}}_2$. Conversely, if any of the other options are chosen, then $({\tau_1},{\tau_2})$ leaves ${\mathcal{X}}_2$ and converges to ${\mathcal{A}}$ asymptotically.

[^10]: Let $r_{\ell_{k}}(\tau)$ be the vector defined by the minimum distance from $\tau$ to the line $\ell_{k}$. Then, it follows that $V(\tau) = (r_{\ell_{k}}^{\top}(\tau)r_{\ell_{k}}(\tau))^{\frac{1}{2}}$. To determine its change during flows, note that on $C\setminus ({\mathcal{X}}\cup {\mathcal{A}})$ the gradient is given by $
    \nabla V(\tau) = \frac{\partial}{\partial \tau} \left( r_{\ell_{k}}^{\top}(\tau)r_{\ell_{k}}(\tau) \right)^{\frac{1}{2}}
    = \frac{\left(r_{\ell_{k}}^{\top}(\tau) \frac{\partial}{\partial \tau} r_{\ell_{k}}(\tau) \right)}{|\tau|_{\ell_{k}}} 
    $ where each $j$-th entry of $\frac{\partial}{\partial \tau} r_{\ell_{k}}(\tau)$ is given by $
    \frac{\partial}{\partial \tau} r^{j}_{\ell_{k}}(\tau) = \frac{\partial}{\partial \tau}\left(({\widetilde}\tau_j\null\hspace{-.08cm}^{k} - \tau_{j}) - \frac{1}{N}\sum_{i=1}^N({\widetilde}\tau_i\null\hspace{-.08cm}^k - \tau_i)^\top \right)
    = \left[\frac{1}{N}, \frac{1}{N},\ldots, \frac{1}{N}, -1 + \frac{1}{N}, \frac{1}{N}, \ldots, \frac{1}{N} \right]
    $ – the term $-1 + \frac{1}{N}$ corresponds to the $j$-th element of the vector. It follows that $\frac{\partial}{\partial \tau} r_{\ell_{k}}(\tau) = \frac{1}{N}{\bf \underline{1}} - {\bf I}$. Then, for each $\tau \in C \setminus {\mathcal{X}}$, $
    \langle \nabla V(\tau), f(\tau) \rangle = \left(\frac{r_{\ell_{k}}^{\top}(\tau)(\frac{1}{N}{\bf \underline{1} - \bf I})}{|\tau|_{\ell_{k}}}\right)f(\tau)
    $.
