# Singularly Perturbed Hybrid Systems for Analysis of Networks with Frequently Switching Graphs

Aneel Tanwani, Hyungbo Shim, Andrew Teel

► To cite this version:

Aneel Tanwani, Hyungbo Shim, Andrew Teel. Singularly Perturbed Hybrid Systems for Analysis of Networks with Frequently Switching Graphs. IEEE Transactions on Automatic Control, 2025, 70 (7), pp.4344 - 4359. 10.1109/TAC.2024.3523242. hal-04855233

HAL Id: hal-04855233

https://laas.hal.science/hal-04855233v1

Submitted on 24 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Singularly Perturbed Hybrid Systems for Analysis of Networks with Frequently Switching Graphs

Aneel Tanwani Hyungbo Shim Andrew R. Teel

**Abstract**—For a class of hybrid systems, where jumps occur frequently, we analyze the stability of system trajectories in view of singularly perturbed dynamics. The specific model we consider comprises an interconnection of two hybrid subsystems, a timer which triggers the jumps, and some discrete variables to determine the index of the jump maps. The flow equations of these variables are singularly perturbed differential equations and, in particular, a smaller value of the singular perturbation parameter leads to an increase in the frequency of the jump instants. For the limiting value of this parameter, we consider a decomposition which comprises a quasi-steady-state system modeled by a differential equation without any jumps and a boundary-layer system described by purely discrete dynamics. Under appropriate assumptions on the quasi-steady-state system and the boundary-layer system, we derive results showing practical stability of a compact attractor when the jumps occur sufficiently often. As an application of our results, we discuss the control design problem in a network of second-order continuous-time coupled oscillators, where each agent communicates the information about its position to some of its neighbors at discrete times. Using the results developed in this article, we show that if the union of the communication graphs being used for information exchange between agents is connected, then the oscillators achieve practical consensus.

**Index Terms**—Hybrid systems; Singular perturbation; Practical stability; Discrete-time switched systems; Heterogeneous agents; Nonlinear oscillators.

## I. INTRODUCTION

The tools developed for the analysis of nonlinear dynamical systems have a rather broad applicability in studying complex systems comprising multiscale dynamics or interconnection of several subsystems. In particular, research in the area of hybrid dynamical systems [9] has greatly benefited from such tools to provide rigorous mathematical solutions to engineering problems involving continuous and discrete dynamics. This article aims at making a new contribution to the research evolving at the intersection of nonlinear analysis and hybrid systems. Motivated by the applications from designing control laws for networked systems with switching topology of the

communication graph, we propose models of hybrid systems where the jump dynamics occur rather frequently. We look at such behavior from the lens of singularly perturbed dynamics and analyze stability in terms of the perturbation parameter.

A common approach for checking the stability of hybrid systems involves finding a Lyapunov function that satisfies some dissipation inequalities along the flow dynamics and jump dynamics over the respective flow and jump sets, see for example [9, Theorem 3.18]. For many applications, it may be hard to compute one such function that satisfies those conditions, and therefore an interesting research direction is to construct this function from the stability information available for the individual components of the hybrid system. Such constructions work under certain conditions which often reveal an interesting interplay between stability of continuous and discrete dynamics. Taking switched systems with reset maps as an example of hybrid systems, it is seen that if the continuous dynamics are stable but discrete dynamics are not, then we can obtain stability conditions by restricting the average frequency of switching times in the form of a lower bound on *dwell-time* on the switching signals (ADT) [21]. On the other hand, if the continuous dynamics are not stable but the discrete dynamics compensate for the growth of the state during continuous-time evolution, then the overall system is asymptotically stable if the number of discrete jumps over each bounded interval is large enough [10]. Lower bounds on the number of discrete jumps necessary for stability are captured by the notion of *reverse-average dwell time* [10]. Using a similar philosophy, some recent works [19] also derive relative upper bounds on the *activation times* for unstable dynamics in the system compared to the lower bounds on the dwell-time for stable dynamics. Another relevant research direction is the stability analysis of interconnected hybrid systems [17]. In these works, the individual subsystems are assumed to have input-to-state stability (ISS) and one can either construct a Lyapunov function for the interconnection, or use trajectory-based analysis for proving stability under appropriate assumptions on gain functions.

This article considers an interconnected hybrid system where one subsystem has stable continuous dynamics, and the other subsystem is modeled by a discrete-time switching system and possibly unstable continuous dynamics. The particular structure that we consider in our work is analyzed using singular perturbation methods. Singular perturbation techniques have been conventionally studied for analysis of ordinary differential equations with slow and fast dynamics [14]. The basic idea is to see the system as a perturbation of a nominal system that comprises the boundary-layer system

The work of A. Tanwani is supported by the project CYPHAI, financed by ANR-JST CREST program, with grant number ANR-20-JSTM-0001. The work of H. Shim was supported by the National Research Foundation of Korea grant funded by MSIT (No. RS-2022-00165417). The work of A.R. Teel was supported in part by AFOSR grant number FA9550-21-1-0452.

A preliminary version of some results presented in this article appeared in the conference paper [28].

A. Tanwani is with CNRS - LAAS, University of Toulouse, France, Email address: aneel.tanwani@cnrs.fr. H. Shim is with ASRI, Department of Electrical and Computer Engineering, Seoul National University, Korea, hshim@snu.ac.kr. A.R. Teel is with Department of Electrical and Computer Engineering, University of California, Santa Barbara, USA, teel@ece.ucsb.edu.
(describing the limit of fast dynamics) and the quasi-steady-state system (described by slow system with limiting value of fast dynamics). This sort of decomposition into slow and fast dynamics has been used for analyzing robustness of input-to-state stability to slowly varying parameters and rapidly varying signals in [29]. More recently, we have seen the use of singular perturbation methods for a class of hybrid systems where fast evolution appears in continuous-time dynamics only [26], [31]. The solutions resulting from a continuous-time boundary-layer system generate an average vector field which approximates the slow subsystem. Using the stability of the average dynamics, with appropriate hypotheses on system data, the authors show semiglobal practical stability for original system. Similarly, singularly perturbed hybrid systems studied in [32] consider the perturbation parameter in continuous dynamics only but use Lyapunov function based methods for analysis.

For the class of interconnected hybrid systems analyzed in this paper, we analyze stability with respect to a parameter in system description which governs the frequency of jump dynamics and plays the role of a singular perturbation parameter. In particular, as the parameter gets smaller, the frequency of jumps increases. This viewpoint was adopted by the authors in [28] and this paper generalizes that framework in several directions with new proof techniques. For such systems, it is natural to consider the decomposition of dynamics in the so-called quasi-steady-state (QSS) system and a boundary-layer (BL) system. With appropriate stability assumptions on these subsystems, we can analyze practical stability of the interconnected system for small enough value of the perturbation parameter. An interesting element of this decomposition is that our boundary-layer system is described by a discrete-time switching system, and the stability of this system is formulated in terms of Lyapunov functions which is consistent with several studies in the literature [1], [7], [8], [23], [25]. For the stability of the overall interconnected hybrid system, we basically consider two proof techniques for our main results. The first approach highlights the underlying principle and is based on the construction of an appropriate Lyapunov function for the interconnection between slow continuous dynamics and fast discrete dynamics. Such constructions provide a generalization of the approach adopted in [28] and are also inspired by the construction provided in [30]. The second approach is motivated by the application studied in this paper, where one cannot readily check the stability of a boundary-layer system with a Lyapunov function. So, differently from [28], we do not construct a Lyapunov function for the interconnection in this case, and instead carry out trajectory-based analysis to get asymptotic bounds on the trajectories.

We see the use of singular perturbation methods for analysis in different applications. One such application of these methods, which is also studied in this article, is the synchronization problem in networked systems. For multi-agent systems with diffusive coupling and continuous-time heterogeneous agents [15], [22], it is observed that the agents achieve practical convergence if the coupling strength is large enough. One can view the coupling parameter as a singular perturbation parameter [16]. More recently, we have also seen the use

of discontinuous coupling for synchronization of networked system in [5]. These discontinuities in the couplings can be interpreted as changes in the graph topology that describes the connections between agents. Consensus problems with time-varying graphs have gathered significant interest in the literature [4], [12] and the techniques based on graph-theoretic methods for analysis of multi-agent systems can be found in [3], [20]. In particular, our work provides analysis of networked systems where the discontinuity in coupling arises due to time-sampled information exchange between neighbors, and the average value of the sampling interval plays the role of the singular perturbation parameter. The change in graph topology at communication times is conveniently modeled using the framework of discrete-time switched systems, and thus motivates the class of systems in (3) which is later used for our case study presented in Section IV. Under natural assumptions on the connectivity of the switching communication graph and dynamics of the agents, we show practical synchronization of the agents’ trajectories under sufficiently frequent communication.

## II. SYSTEM CLASS

In this paper, we consider hybrid dynamical systems described by ordinary differential inclusions over a flow set, and difference inclusions over a jump set. The peculiar feature of our models is that they are singularly perturbed, that is, there is a parameter $\epsilon > 0$ which affects the qualitative behavior of the trajectories and we analyze stability of the system for small enough values of this parameter.

### A. Timer with frequent jumps

An important modeling aspect for us is to consider a timer with the following dynamics

$$ \begin{aligned} \epsilon\dot{\tau} &\in [\sigma_1, \sigma_2], & \tau &\in [0, N_0] \\ \tau^+ = \tau - 1, & & \tau &\in [1, N_0]. \end{aligned} \quad (1) $$

The evolution of the timer $\tau$, determined by the parameters $N_0 \in \mathbb{N}$, $0 < \epsilon \le 1$, and $0 < \sigma_1 \le \sigma_2$, describes the interplay between the flow and the jump. In our work, we will consider hybrid systems where the jumps are triggered only by the timer $\tau$. We are particularly interested in the case where the parameter $\epsilon$ is very small, and we refer to $\epsilon$ as the singular perturbation parameter.

In the literature on hybrid systems, we find different setups where time-triggered jumps play a crucial role in the stability of the hybrid system. For example, the notion of (average) dwell-time [11], [21] is used for putting an upper bound on how frequently the jumps can occur, whereas the notion of reverse average dwell-time [10] is used for imposing a lower bound on the frequency of jumps. The timer in (1) has the property that it simultaneously puts a lower and upper bound on the number of jumps over an interval $(s, t]$, that is,

$$ 1 + \frac{\sigma_1(t-s)}{\epsilon} - N_0 \le J_{(s,t]} \le N_0 + \frac{\sigma_2(t-s)}{\epsilon} \quad (2) $$

where $J_{(s,t]}$ denotes the number of jumps in the variable $\tau$ over the interval $(s, t]$. It readily follows that, with such a timer,
there is no accumulation time for jumps (i.e., there are finitely many jumps in any finite continuous-time interval) and the domain of every complete solution has unbounded continuous-time element.

## B. Overall model

We are interested in modeling the evolution of two variables $(x, y) \in \mathbb{R}^n \times \mathbb{R}^m$. The dynamics of $x$ variable are modeled by an ordinary differential equation. For the evolution of $y$-dynamics, in addition to a flow map described by a differential equation, we have a collection of jump maps to define the jump dynamics. The times at which these jumps occur is modeled by the timer $\tau$ in (1) (that is, the jumps occur whenever $\tau$ jumps), and the jump map for $y$ is indexed by a variable $p$ belonging to a finite set $\mathcal{P}$.

The state variables for the overall system are $(x, y, \tau, p) \in \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R} \times \mathcal{P}$, and we consider a closed set $C_{xy}$ contained in an open domain $\mathcal{D}_{xy} \subset \mathbb{R}^n \times \mathbb{R}^m$. The state variables evolve inside the flow set

$$C := C_{xy} \times [0, N_0] \times \mathcal{P},$$

and go through an instantaneous change inside the jump set

$$\mathcal{D} := C_{xy} \times [1, N_0] \times \mathcal{P}.$$

The underlying equations that describe the evolution of the state variables are:

$$\left\{
\begin{array}{l}
\dot{x} = f_x(x, y, \epsilon), \\
\dot{y} = f_y(x, y, \epsilon), \\
\epsilon \dot{\tau} \in [\sigma_1, \sigma_2], \\
\dot{p} = 0,
\end{array}
\right.
\quad \text{for } (x, y, \tau, p) \in \mathcal{C}, \quad (3a)$$

$$\left\{
\begin{array}{l}
x^+ = x, \\
y^+ = g_y(x, y, p), \\
\tau^+ = \tau - 1, \\
p^+ \in G_p(p),
\end{array}
\right.
\quad \text{for } (x, y, \tau, p) \in \mathcal{D}. \quad (3b)$$

For the timer dynamics, the positive constants $\sigma_1$ and $\sigma_2$ are such that $\sigma_2 \ge \sigma_1 > 0$. We assume that the vector fields $f_x : \mathcal{D}_{xy} \times \mathbb{R}_{\ge 0} \to \mathbb{R}^n$, $f_y : \mathcal{D}_{xy} \times \mathbb{R}_{\ge 0} \to \mathbb{R}^m$ and the jump maps $g_y : \mathcal{D}_{xy} \times \mathcal{P} \to \mathbb{R}^m$ are continuous. The set-valued map $G_p : \mathcal{P} \rightrightarrows \mathcal{P}$ allows non-deterministic behavior for the variable $p$ at jumps. The solutions of system (3) are defined over a subset of $\mathbb{R}_{\ge 0} \times \mathbb{N}$, called *hybrid time domain* [9, Chap. 2]. In what follows, we use the notation $(x_0, y_0) := (x(0, 0), y(0, 0))$ to denote the initial condition for $(x, y)$-dynamics.

For the problems studied in this paper, the parameter $\epsilon$ appearing in (3) is taken to be small, which implies that the jumps occur frequently, and thus, the evolution of $y$ relies heavily on the jump dynamics. Since the jumps correspond to rapid evolution, we call $x$ the slow variable and $y$ the fast variable.

## C. System decomposition

We are interested in studying the behavior of system (3) when the parameter $\epsilon$ is sufficiently small and takes values close to zero. It turns out that the evolution of trajectories

can be approximated by two different dynamical systems. The first system is called the *boundary-layer system*, which is given by

$$\begin{align}
x^+ &= x \nonumber \\
y^+ &= g_y(x, y, p) \tag{4} \\
p^+ &\in G_p(p) \nonumber
\end{align}$$

for $(x, y) \in C_{xy}$, and $p \in \mathcal{P}$. Recalling the expression of the timer dynamics, we observe that a small value of $\epsilon$ corresponds to many jumps in the trajectories of system (3) within a small continuous-time interval. This means that, from an initial condition $(x(0, 0), y(0, 0))$, while the solution $x(\cdot)$ remains close to $x(0, 0)$, the solution $y(\cdot)$ evolves like a solution to the discrete-time system (4). In order to define an equilibrium point for the discrete-time system (4), we let $C_x$ denote the natural projection of $C_{xy}$ on $\mathbb{R}^n$, that is, $C_x := \{x \in \mathbb{R}^n : \exists y \text{ s.t. } (x, y) \in C_{xy}\}$; and we stipulate that there is a continuously differentiable function $h : \mathbb{R}^n \to \mathbb{R}^m$ such that

$$g_y(x, h(x), p) = h(x), \quad \forall x \in \bar{C}_x, \quad \forall p \in \mathcal{P}, \quad (5)$$

where $\bar{C}_x$ denotes the closure of $C_x$. We assume that $(x, h(x)) \in C_{xy}$ for all $x \in C_x$ in order for the problem to be well-posed. Imagine that $y(t, j)$ rapidly approaches $h(x(0, 0))$ after some jumps, and remains near $h(x(t, j))$ as time goes on. Then, the behavior of $x(t, j)$ can be approximated by the *quasi-steady-state* system defined as

$$\dot{x} = f_x(x, h(x), 0), \quad x \in \bar{C}_x, \quad (6)$$

which is the evolution of $x(\cdot)$ on the slow manifold $\{(x, y) \in \mathbb{R}^n \times \mathbb{R}^m : y = h(x)\}$, when $\epsilon = 0$. We associate a compact forward invariant set $A_x$ with (6) so that if $x(0) \in A_x$, then for every solution of (6), we have $x(t) \in A_x$.

## D. Stability notions

The primary focus of this paper is the stability analysis of system (3). In particular, we are interested in developing bounds on the $(x, y)$ trajectories, when the initial conditions belong to a given compact set and the perturbation parameter $\epsilon$ tends to 0. This stability notion is formalized in the following definition:

**Definition 1.** For the hybrid system (3), with boundary-layer system (4) and quasi-steady-state system (6), the compact set $A := \{(x, y) \in C_{xy} : x \in A_x, y = h(x)\} \times [0, N_0] \times \mathcal{P}$ is two-time-scale semiglobally practically stable, or TSPS for short, if for a given compact set $I_{xy} \subset C_{xy}$, there are a positive number $\epsilon^*$, class $KL$ functions $\beta_x$ and $\beta_y$, and class $K$ functions $\chi_x$ and $\chi_y$ such that, for any $0 < \epsilon \le \epsilon^*$ and for any $(x_0, y_0) \in I_{xy}$, the solution $(x(t, j), y(t, j))$ to (3) satisfies the following,

$$|x(t, j)|_{A_x} \le \beta_x(|x_0|_{A_x}; t) + \chi_x(\epsilon), \quad (7a)$$

$$|y(t, j) - h(x(t, j))| \le \beta_y(|y_0 - h(x_0)|, t/\epsilon) + \chi_y(\epsilon), \quad (7b)$$

$$\text{for all } (t, j) \in \text{dom}(x, y, \tau, p). \quad \triangleleft$$

For the results developed in this paper, we impose some stability assumptions on the boundary-layer system and the
quasi-steady-state system. This will allow us to establish the property stated in Definition 1 for the coupled system (3).

**Assumption 1.** The quasi-steady-state system (6) admits an asymptotically stable compact attractor $A_x \subset C_x$. △

**Assumption 2.** For every compact set $K_x \subset C_x$, the dynamical system (4), restricted to the set $\mathcal{D}_{\text{bl}} := \{(x, y, p) \in C_{xy} \times \mathcal{P} : x \in K_x\}$, admits a globally asymptotically stable compact set<fnref n="1" />

$$A_y := \{(x, y) \in C_{xy} : x \in K_x, y = h(x)\} \times \mathcal{P}, \quad (8)$$

where $h$ satisfies (5). △

### III. MAIN RESULTS

The main contribution of this article is to show that Assumption 1 and Assumption 2 together allow us to establish stability of system (3) in the sense of Definition 1 with small enough values of $\epsilon$.

**Theorem 1.** For the hybrid system (3), suppose that the quasi-steady-state system (6) satisfies Assumption 1 and the boundary-layer system (4) satisfies Assumption 2. Then, the set $A$ is TSPS.

To prove this result, we need to derive the inequalities (7a) and (7b). We will develop two different approaches for derivation of these inequalities, and hence two different techniques for proving Theorem 1. The difference in the two approaches arises due to the information that is available to us for checking Assumption 2, which is consequently used for the derivation of (7). In the first approach, we have a true Lyapunov function for the boundary-layer system which decreases along the dynamics (4) at each (jump) instant. In the second approach, we only have relaxed Lyapunov functions which decrease along the dynamics (4) after multiple (jump) instants. Since these two types of Lyapunov functions describe possibly different behavior of the boundary-layer system (4), the derivation of the inequalities in (7b) is also very different from each other. On the other hand, the derivation of (7a) in two approaches is similar.

#### A. Analysis with true Lyapunov function for boundary-layer system

Let us consider the following variants of Assumption 1 and Assumption 2 in terms of Lyapunov functions. These Lyapunov functions will then be used for the derivation of inequalities appearing in (7).

**Assumption 1a.** There exists a continuously differentiable function $V_x : \mathcal{D}_x = \{x : \exists y \text{ s.t. } (x, y) \in \mathcal{D}_{xy}\} \to \mathbb{R}$ such that, for all $x \in \bar{C}_x$,

$$\underline{\alpha}_x(|x|_{A_x}) \le V_x(x) \le \bar{\alpha}_x(|x|_{A_x}), \quad (9)$$

$$\nabla_x V_x(x) \cdot f_x(x, h(x), 0) \le -\alpha_x(V_x(x)), \quad (10)$$

where $\underline{\alpha}_x$ and $\bar{\alpha}_x$ are class $K_\infty$ functions, and $\alpha_x$ is a class $K$ function.

Following Assumption 1a, we use [18, Lemma 4.4] to get a class $KL$ function $\beta_x^0$ from $\alpha_x$ appearing in (10). Let<fnref n="2" />

$$\eta_x(s) := - \int_1^s \frac{dr}{\min\{r, \alpha_x(r)\}}, \quad \text{for } s > 0,$$

and

$$\beta_x^0(s, t) := \begin{cases} \eta_x^{-1}(\eta_x(s) + t), & \text{for } s > 0, \\ 0, & \text{for } s = 0. \end{cases} \quad (11)$$

For the boundary-layer system (4), the existence of true Lyapunov function is formalized as follows:

**Assumption 2a.** For every compact set $K_x \subset C_x$, there exist a continuously differentiable functions $V_y : \mathcal{D}_{xy} \times \mathcal{P} \to \mathbb{R}$, class $K_\infty$ functions $\underline{\alpha}_y$ and $\bar{\alpha}_y$, and a non-negative constant $\gamma_c < 1$ that satisfy the following items:

(BL-C1) for every $(x, y, p) \in \{(x, y) \in C_{xy} : x \in K_x\} \times \mathcal{P}$,

$$\underline{\alpha}_y(|y - h(x)|) \le V_y(x, y, p) \le \bar{\alpha}_y(|y - h(x)|),$$

(BL-C2) for every $(x, y, p) \in \{(x, y) \in C_{xy} : x \in K_x\} \times \mathcal{P}$, and for every $q \in G_p(p)$,

$$V_y(x, g_y(x, y, p), q) \le \gamma_c V_y(x, y, p).$$

**Remark 1.** For a fixed compact set $K_x$, it can be checked that Assumption 2a implies Assumption 2 as one would do so for proving stability using Lyapunov functions. Also, from converse Lyapunov theorems for discrete-time switching systems, one can show that Assumption 2a is also necessary for Assumption 2 to hold due to continuity assumption on the mappings $g_y$, closedness of the set $C_{xy}$, and the compactness of the attractor $A_y$, see for instance [8] and [13] for one possible demonstration of such a result. By similar reasoning, one can also check that Assumption 1 and Assumption 1a are equivalent. △

In the light of Remark 1, Theorem 1 is qualitatively equivalent to Theorem 2 stated next, but here, we provide more information about the functions appearing on the right-hand side of the inequalities in (7), using the information in Assumption 1a and Assumption 2a.

**Theorem 2.** For the hybrid system (3), suppose that the quasi-steady-state system (6) satisfies Assumption 1a and the boundary-layer system (4) satisfies Assumption 2a. Then, the set $A$ is TSPS, and the estimate (7a) holds with

$$\beta_x(s, t) := \underline{\alpha}_x^{-1} \left( \max \left\{ \beta_x^0(2\beta_x^0(\bar{\alpha}_x(s), 1), \max\{0, t-1\}), \beta_x^0(\bar{\alpha}_x(s), t) \right\} \right) \quad (12)$$

and (7b) holds with

$$\beta_y(s, \tau) = \underline{\alpha}_y^{-1} \left( 2\gamma_c^{\frac{\sigma_1}{N_0}\tau} \bar{\alpha}_y(s) \right). \quad (13)$$

The expressions for class $K$ functions $\chi_x$ and $\chi_y$ are derived in the proof of Theorem 2, which is carried out in

<footnote n="1">¹Here, we understand asymptotic stability of the set $A_y$ in the sense of [9, Definition 3.6] by regarding system (4) as a hybrid system with purely discrete dynamics and the jump set equal to $\mathcal{D}_{\text{bl}}$.</footnote>

<footnote n="2">²One can also take $\eta_x(s) = -\int_1^s \frac{dr}{\alpha_x(r)}$, for $s > 0$, if it holds that $\lim_{s \to 0^+} \eta_x(s) = +\infty$.</footnote>
Section V. An element of independent interest in the proof is the construction of a Lyapunov function for (3) to show boundedness of trajectories, whilst using some of the ideas from [30].

## B. Analysis with relaxed Lyapunov function for boundary-layer system

In Theorem 2, we made use of Assumption 2a which stipulates the existence of Lyapunov functions that decrease at each jump instant along the dynamics of the boundary-layer system. For certain applications, such as the one considered in Section IV, we may only have Lyapunov functions which decrease after a certain number of jumps. The existence of such relaxed Lyapunov functions is formulated in forthcoming Assumption 2b, and it will be used along with Assumption 1a to provide a different technique for the derivation of inequalities in (7). Notably, the functions on the right-hand side of (7b) depend on the information contained in Assumption 2b.

For convenience, let us introduce the following notation: for each $(x, y, p) \in C_{xy} \times \mathcal{P}$, let

$$ (g_y(x, y, p), G_p(p))^r := \{(y, q) : (y, q) = (y_r, p_r), (y_0, p_0) = (y, p), (y_{i+1}, p_{i+1}) \in (g_y(x, y_i, p_i), G_p(p_i)), 0 \le i \le r-1\}. $$

**Assumption 2b.** For every compact set $K_x \subset C_x$, there exist a positive integer $n$, continuously differentiable functions $V_y: \mathcal{Q}_{xy} \times \mathcal{P} \to \mathbb{R}$, class $K_\infty$ functions $\underline{\alpha}_y$ and $\overline{\alpha}_y$, and a non-negative constant $\gamma_m < 1$ that satisfy the following items:

(BL-M1) for every $(x, y, p) \in \{(x, y) \in C_{xy} : x \in K_x\} \times \mathcal{P}$,

$$ \underline{\alpha}_y(|y - h(x)|) \le V_y(x, y, p) \le \overline{\alpha}_y(|y - h(x)|), $$

(BL-M2) there exists a set $P \subset \mathcal{P}$, such that,

* every solution of $(p \in \mathcal{P}, p^+ \in G_p(p))$ reaches $P$ in at most $n$ steps,

* for every $(x, y, p) \in \{(x, y) \in C_{xy} : x \in K_x\} \times P$, there exists $r \in \{1, \dots, n\}$, such that for every $(g, q) \in (g_y(x, y, p), G_p(p))^r$, we have that $q \in P$ and

$$ V_y(x, g, q) \le \gamma_m V_y(x, y, p). \quad (14) $$

**Remark 2.** With continuous functions $g_y$, we can also show that Assumption 2b implies Assumption 2. In fact, this implication can be proven using the arguments given in the proof of Theorem 3 that appears in Section VI. Moreover, the reader may check that if Assumption 2a holds, then we can construct the functions stipulated in Assumption 2b from $V_y$ satisfying (BL-C1), (BL-C2), and the mappings $g_y$. Thus, taking Remark 1 into account, we observe that the three statements appearing in Assumption 2, Assumption 2a and Assumption 2b are all equivalent under the given regularity assumptions on the data of system (4). We refer the interested readers to [1], [6], [8], [23] for further details along the lines of analyzing stability of discrete-time switched systems using Lyapunov functions.

$\triangleleft$

**Theorem 3.** For the hybrid system (3), suppose that the boundary-layer system (4) satisfies Assumption 2b and the quasi-steady-state system (6) satisfies Assumption 1a. Then,

the set $\mathcal{A}$ is TSPS, and the estimate (7a) holds with $\beta_x(s, t)$ given in (12), and (7b) holds with

$$ \beta_y(s, \tau) := C_g \circ \underline{\alpha}_y^{-1} \left( 2 \gamma_m^{\frac{\sigma_1}{N_0 n} \tau} \overline{\alpha}_y(s) \right), \quad (15) $$

where $\gamma_m$ appears in (14), and $C_g$ is a class $\mathcal{K}$ function.

The proof of Theorem 3 is carried out in Section VI, and it contains the expressions for corresponding class $\mathcal{K}$ functions $\chi_x$ and $\chi_y$. The main idea of the proof is to get the bounds on the trajectories by looking at the perturbation of boundary-layer system (4) and the function $C_g$ in (15) essentially captures the size of perturbations.

## C. Boundary-layer system with constrained switching

A natural motivation for Assumption 2b is to analyze boundary-layer systems with constrained switching where we do not have a common Lyapunov function. To explain how Assumption 2b can be used in such cases, we consider a finite discrete set $\mathcal{Q}$, and a family of continuous mappings $g_y(\cdot, \cdot, q) : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^m$, with $q \in \mathcal{Q}$. The mappings $g_y(\cdot, \cdot, q)$ indexed by $q \in \mathcal{Q}$ would be used in the description of the discrete dynamics, and we are particularly interested in the case where the transition between the jump maps respects some pre-specified constraints. To describe these constraints, we consider a set $\mathcal{S}$ which comprises strings (also called words) of finite length formed by the alphabet $\mathcal{Q}$. For a string $s \in \mathcal{S}$, we use the notation $s_r \in \mathcal{Q}$ to denote the r-th letter of that string $s$, with $1 \le r \le |s|$, and $|s|$ denoting the total number of letters in $s$. We also consider the set-valued mapping $\text{succ}: \mathcal{S} \stackrel{\sim}{\rightleftharpoons} \mathcal{S}$, where $\text{succ}(s) \subset \mathcal{S}$ denotes a nonempty set of admissible successors of $s$. We emphasize that the set $\mathcal{Q}$ is used as index for different jump maps, whereas the set $\mathcal{S}$ is used to encode the constraints on switching sequence which respect certain assumptions related to stability. We assume that $\mathcal{S}$ is a finite set, and let $n := \max_{s \in \mathcal{S}} |s|$.

Now, we want to model our discrete dynamics with constraints on the switching sequence in the form of a conventional one-step difference equation by using certain discrete variables. Indeed, we suppose the jump map $g_y(\cdot, \cdot, q)$ is selected by the rule $q = s_r$, where the string $s \in \mathcal{S}$ and the indexing variable $r \in \mathbb{N}$ satisfy the following<fnref n="3" />:

$$ s^+ \in S_1(s, r) = \begin{cases} \{s\}, & \text{if } r < |s|, \\ \text{succ}(s), & \text{otherwise,} \end{cases} \quad (16a) $$

$$ r^+ \in S_2(s, r) = \begin{cases} \{r+1\}, & \text{if } r < |s|, \\ \{1\}, & \text{otherwise,} \end{cases} \quad (16b) $$

and $\text{succ}(s) \subset \mathcal{S}$ is a nonempty set of successors of $s \in \mathcal{S}$. We can now apply Assumption 2b to the boundary-layer system (4) with $p = (s, r)$, $\mathcal{P} = \mathcal{S} \times \{1, \dots, n\}$, and choosing $P = \mathcal{S} \times \{1\}$. This will provide us the sufficient conditions for asymptotic stability of switched systems with multiple Lyapunov functions in the same spirit as [1].

<footnote n="3">³The cardinality of the variable $s$ in (16a) may vary, so a precise description of the discrete dynamics for $(s, r)$ would require us to embed the variables in an Euclidean space of fixed dimension. This can be done by considering the variable $\bar{s}$ of fixed length $n$ which is obtained from $s$ by appending letters not in $\mathcal{Q}$.</footnote>
IV. CASE STUDY: NETWORK OF OSCILLATORS

As an application of Theorem 3, we study the problem of synchronization of coupled nonlinear oscillators. In particular, we consider N continuous-time controlled oscillators, which are connected to each other through an undirected graph. The control action for each oscillator is chosen as a function of the states of their neighbors at some discrete times. If this control is updated sufficiently often then we show that all the trajectories of the oscillators converge to each other in practical sense.

Each of the N oscillators (the agent) has a two dimensional state $(\omega_i, \xi_i) \in \mathbb{R}^2$, for $i = 1, \dots, N$, whose dynamic equation is written in the form of Liénard system as

$$ 
\begin{aligned}
\dot{\omega}_i &= -\omega_i + \xi_i \\
\dot{\xi}_i &= \phi_i(\omega_i, \xi_i) := (1 - f_i(\omega_i))(-\omega_i + \xi_i) - g_i(\omega_i)
\end{aligned}
\qquad (17) $$

where the functions $f_i$ and $g_i$ are assumed to be locally Lipschitz. For example, if the agent $i$ has $f_i(\omega_i) = \mu_i(\omega_i^2 - 1)$ and $g_i(\omega_i) = \omega_i$ with a constant $\mu_i > 0$, then the agent $i$ is a well-known Van der Pol oscillator [14]. If we look at an individual oscillator (17), it may have a stable or unstable limit cycle. In what follows, we consider a network of such oscillators, where each oscillator exchanges information about its state $\xi$ only with its neighbors at some discrete times, and consequently, $\xi$ gets reset whenever some new information is received from the neighbors. Our objective is to show that, with an appropriate reset rule, if the information between the neighbors is exchanged sufficiently often then all the oscillators in the network practically synchronize and converge practically to the same limit cycle.

A. Connectivity graphs

We consider a set of N vertices $V = \{1, \dots, N\}$, and let $Q$ be the index set for different collections of edges associated to these nodes, that is, for each $q \in Q$, $\mathcal{E}_q$ is a set of *unordered* pair of vertices that represents one possible collection of edges. An element in $\mathcal{E}_q$ is represented by $(i, j)$, for some choice of $i, j \in V$, $i \neq j$. We call $\mathcal{G}_q = (V, \mathcal{E}_q)$ a graph with nodes $V$ and edges $\mathcal{E}_q$. Here, the graph $\mathcal{G}_q$ is *undirected*. Agents $i$ and $j$ are *neighbors* of each other in graph $\mathcal{G}_q$ if $(i, j) \in \mathcal{E}_q$. The neighbors of a node $i \in V$ in a graph $\mathcal{G}_q$ are denoted by $\mathcal{N}_i^q$. Let $S$ be a set of words obtained from the alphabet $Q$ and denote an element of $S$ by $s$. We say that the graphs $\{\mathcal{G}_q\}_{q \in S}$ are *jointly connected* if $\mathcal{G}_s := (V, \mathcal{E}_s)$ is connected where $\mathcal{E}_s := \cup_{q \in s} \mathcal{E}_q$, in which 'q in s' implies $q \in \{s_1, \dots, s_{|s|}\}$.

**Assumption 3.** The set $S$ is finite, and for each $s \in S$, the collection of graphs $\{\mathcal{G}_q\}_{q \in s}$ is jointly connected.<fnref n="4" /> △

In our setup, the continuous-time oscillators described in (17) represent the nodes which communicate with another neighboring oscillator only at some discrete time instants. The edges of a graph $\mathcal{G}_q$ represent communication link between the corresponding nodes. The motivation for introducing different

graphs $\mathcal{G}_q$ is to allow for different communication patterns at different times among the oscillators.

With each $\mathcal{G}_q$, we associate a doubly stochastic matrix $W_q \in \mathbb{R}^{N \times N}$. If at a given time instant, the graph $\mathcal{G}_q$ is active, then the state $\xi_i$ is updated at that instant as a function of its neighbors' state $\xi_j$, $j \in \mathcal{N}_i^q$. In particular, we describe the update model as

$$ 
\begin{aligned}
\omega_i^+ &= \omega_i \\
\xi_i^+ &= \sum_{j \in \mathcal{N}_i^q} w_{ij}^q \xi_j
\end{aligned}
\qquad (18) $$

where $w_{ij}^q$ is the $(i, j)$-th entry of $W_q$. We assume that the graph $\mathcal{G}_q$ contains the self-loop for every node, so that $i \in \mathcal{N}_i^q$, $\forall i$, and all the diagonal entries of $W_q$ are positive.

**Example 1.** If the graph $\mathcal{G}_q$ is modeled by a symmetric adjacency matrix $A_q$, and $D_q$ is the diagonal matrix that describes the degree of each node, then $W_q$ can be given by $I_N - \varphi D_q + \varphi A_q$, where $\varphi > 0$ is chosen small enough such that $W_q$ becomes a non-negative matrix. Then, every row and column sums of $W_q$ is 1 by construction. △

**Example 2.** Consider a connected graph $\mathcal{G} = (V, \mathcal{E})$. Let $E$ denote the total number of edges in $\mathcal{E}$, and let $\mathcal{E}_q$ be a singleton which corresponds to a single edge contained in $\mathcal{E}$. We let $Q = \{1, \dots, E\}$ and consider strings $s$ from the alphabet $Q$ such that the collection of graphs $\{\mathcal{G}_q = (V, \mathcal{E}_q)\}_{q \in s}$ is jointly connected. Under connectivity assumption on $\mathcal{G} = (V, \mathcal{E})$, the set of such strings is nonempty and one possible choice is to take the strings $s$ as words of length $E$ with distinct letters from the alphabet $Q$. In this case, we can choose $W_q \in \mathbb{R}^{N \times N}$, $q \in Q$ to be a doubly stochastic matrix such that, for the index $(i, j) \in \mathcal{E}_q$,

$$ w_{ij}^q = w_{ji}^q = w_{ii}^q = w_{jj}^q = \frac{1}{2} $$

and $w_{kk}^q = 1$ for all $k \notin \{i, j\}$, and all other elements of $W_q$ are zero. △

Thus far, in this subsection, we have introduced a family of graphs $\mathcal{G}_q = (V, \mathcal{E}_q)$ and a reset rule (18) for individual agents which depends on the neighbors described by $\mathcal{G}_q$. To specify the times at which an agent resets its state and to describe the index of the active graph, we consider a timer $\tau$ and a discrete state variable $p$, whose evolution is described as,

$$ (\tau, p) \in [0, N_0] \times \mathcal{P} : \left\{ \begin{array}{l} \epsilon \dot{\tau} \in [\sigma_1, \sigma_2] \\ \dot{p} = 0 \end{array} \right. \qquad (19) $$

$$ (\tau, p) \in [1, N_0] \times \mathcal{P} : \left\{ \begin{array}{l} \tau^+ = \tau - 1 \\ p^+ \in G_p(p) \end{array} \right. \qquad (20) $$

where $N_0$, $\sigma_1$, and $\sigma_2$ are some positive constants that determine the frequency at which the graphs switch from one configuration to another. The mapping $G_p$ describes the evolution of discrete variables according to equation (16).

B. Coordinate transformation

We now collect the dynamics of individual agents from (17) and the connectivity structure from Section IV-A, to write them in the form of the hybrid system (3). This requires us

<footnote n="4">⁴We consider the case of an undirected graph for simple presentation, but the study in this section can be transparently extended to the case where $\{\mathcal{G}_q\}_{q \in s}$ is jointly strongly connected and balanced.</footnote>
to introduce a certain coordinate transformation for agents’ dynamics which we describe next.

Let $\mathbf{R} \in \mathbb{R}^{N \times (N-1)}$ be a matrix such that

$$ \mathbf{R}^{\top} \mathbf{R} = I, \quad \text{and} \quad \mathbf{1}_{N}^{\top} \mathbf{R} = 0, $$

where $\mathbf{1}_N \in \mathbb{R}^N$ is a vector with all entries equal to 1. Associated to each matrix $W_q$, let us introduce the matrix $\Lambda_q$ as follows:

$$ \Lambda_q := \mathbf{R}^{\top} W_q \mathbf{R} \in \mathbb{R}^{(N-1) \times (N-1)}. $$

Next, we define a matrix $\mathbf{Q}$ as

$$ \mathbf{Q}^{-1} := [\mathbf{1}_N \quad \mathbf{R}], \quad \text{so that} \quad \mathbf{Q} = \begin{bmatrix} \frac{1}{N} \mathbf{1}_N^{\top} \\ \mathbf{R}^{\top} \end{bmatrix}. $$

With the matrix $\mathbf{Q}$, we can write (17) and (18) in the form of (3). Let us introduce the new coordinates $(s, \zeta)$ for $\xi := \text{col}(\xi_1, \dots, \xi_N) \in \mathbb{R}^N$ as

$$ \begin{pmatrix} s \\ \zeta \end{pmatrix} := \mathbf{Q}\xi = \begin{bmatrix} \frac{1}{N} \mathbf{1}_N^{\top} \\ \mathbf{R}^{\top} \end{bmatrix} \xi \quad (21) $$

or, in other words

$$ s = \frac{1}{N} \sum_{i=1}^{N} \xi_i \in \mathbb{R} \quad \text{and} \quad \zeta = \mathbf{R}^{\top} \xi \in \mathbb{R}^{N-1}. $$

It also follows that $\xi = \mathbf{1}_N s + \mathbf{R}\zeta$. Let $\mathbf{R}_i$ be the $i$-th row of $\mathbf{R}$. Then, the system (17), (18) can be written as

$$ \begin{align} \dot{\omega}_i &= -\omega_i + s + \mathbf{R}_i \zeta, \quad i = 1, \dots, N, \notag \\ \dot{s} &= \frac{1}{N} \sum_{i=1}^{N} (1 - f_i(\omega_i))(-\omega_i + s + \mathbf{R}_i \zeta) - g_i(\omega_i) \tag{22a} \\ \dot{\zeta} &= \mathbf{R}^{\top} \phi(\omega, \mathbf{1}_N s + \mathbf{R}\zeta) \notag \\ \epsilon \dot{\tau} &\in [\sigma_1, \sigma_2] \notag \\ \dot{p} &= 0, \notag \end{align} $$

where $\omega := \text{col}(\omega_1, \dots, \omega_N)$ and

$$ \phi(\omega, \mathbf{1}_N s + \mathbf{R}\zeta) = \begin{bmatrix} \phi_1(\omega_1, s + \mathbf{R}_1\zeta) \\ \vdots \\ \phi_N(\omega_N, s + \mathbf{R}_N\zeta) \end{bmatrix} $$

with the flow set $\mathcal{C} := \mathbb{R}^{2N} \times [0, N_0] \times \mathcal{P}$, and

$$ \begin{align} \omega_i^+ &= \omega_i, \quad \text{for } i = 1, \dots, N \notag \\ s^+ &= s \notag \\ \zeta^+ &= \Lambda_p \zeta \notag \\ \tau^+ &= \tau - 1 \notag \\ \mathbf{p}^+ &\in G_p(\mathbf{p}) \tag{22b} \end{align} $$

where $\Lambda_p$ is the notation for $\Lambda_{s_r}$ considering $\mathbf{p} = (\mathbf{s}, \mathbf{r})$, with the jump set $\mathcal{D} := \mathbb{R}^{2N} \times [1, N_0] \times \mathcal{P}$. Here, if we treat $(\omega_1, \dots, \omega_N, s) = (\omega, s)$ as the slow variable and $\zeta$ as the fast variable with both the flow set and the jump set for $(\omega, s, \zeta)$ being the same as $\mathbb{R}^{2N}$, then the system is in the form of (3).

### C. Quasi-steady-state system

For stability analysis, we now consider the QSS system associated to (22). The QSS system is

$$ \begin{align} \dot{\omega}_i &= -\omega_i + s, \quad i = 1, \dots, N, \notag \\ \dot{s} &= \frac{1}{N} \sum_{i=1}^{N} (1 - f_i(\omega_i))(-\omega_i + s) - g_i(\omega_i), \tag{23} \end{align} $$

which corresponds to (6). The following assumption and lemma assure that this QSS system<fnref n="5" /> satisfies the assumption of Theorem 3.

**Assumption 4.** The second order dynamical system

$$ \begin{align} \dot{\chi}_0 &= -\chi_0 + s_0 \notag \\ \dot{s}_0 &= (1 - f(\chi_0))(-\chi_0 + s_0) - g(\chi_0) \tag{24} \end{align} $$

where

$$ f(\chi_0) := \frac{1}{N} \sum_{i=1}^{N} f_i(\chi_0), \quad g(\chi_0) := \frac{1}{N} \sum_{i=1}^{N} g_i(\chi_0) \quad (25) $$

has a locally asymptotically stable limit cycle $\mathcal{A}_0 \subset \mathbb{R}^2$. $\triangleleft$

The limit cycle $\mathcal{A}_0$ in Assumption 4 refers to an isolated periodic orbit, that is, every solution $(\chi_0, s_0)$ of system (24) with $(\chi_0(0), s_0(0)) \in \mathcal{A}_0$ is periodic, and $(\chi_0(t), s_0(t)) \in \mathcal{A}_0$, for all $t \ge 0$. Moreover, $\mathcal{A}_0$ is an isolated and compact set. Local asymptotic stability of $\mathcal{A}_0$ ensures that the nearby trajectories converge to this limit cycle in some neighborhood of $\mathcal{A}_0$. A sufficient condition on $f$ and $g$ for Assumption 4 can be found in [24, Sec. 3.8]. We use this assumption to show the existence of a stable limit cycle for system (23).

**Lemma 1.** Under Assumption 4, the QSS system (23) has a locally asymptotically stable compact set $\mathcal{A}_x \subset \mathbb{R}^{N+1}$, which is a stable limit cycle for (23), given by

$$ \mathcal{A}_x = \{(\chi_0, \dots, \chi_0, s_0) \in \mathbb{R}^{N+1} : (\chi_0, s_0) \in \mathcal{A}_0\}. \quad (26) $$

In particular, there exists $V_x$ that satisfies (9) and (10) on a neighborhood $\Omega_x$ of $\mathcal{A}_x$.

*Proof.* The basic idea behind the proof is to use a coordinate transformation and combine it with the cascade arguments [27, Proposition 4.1]. Let

$$ \chi := \frac{1}{N} \mathbf{1}_N^{\top} \omega \in \mathbb{R}, \quad \tilde{\chi} := \mathbf{R}^{\top} \omega \in \mathbb{R}^{N-1} $$

so that $\omega_i = \chi + \mathbf{R}_i \tilde{\chi}$. Then, system (23) can be rewritten as

$$ \begin{align} \dot{\tilde{\chi}} &= -\tilde{\chi}, \notag \\ \dot{\chi} &= -\chi + s, \notag \\ \dot{s} &= \frac{1}{N} \sum_{i=1}^{N} (1 - f_i(\chi + \mathbf{R}_i \tilde{\chi})) (s - \chi - \mathbf{R}_i \tilde{\chi}) - g_i(\chi + \mathbf{R}_i \tilde{\chi}) \notag \\ &= (1 - f(\chi))(-\chi + s) - g(\chi) + \frac{1}{N} \sum_{i=1}^{N} \psi_i(\tilde{\chi}, \chi, s), \tag{27} \end{align} $$

<footnote n="5">⁵The QSS system (23) is called ‘blended dynamics’ in the context of heterogeneous multi-agent systems [16].</footnote>
where $\psi_i$ vanishes when $\tilde{\chi} = 0$ because

$$ \psi_i = f_i(\chi + R_i\tilde{\chi})R_i\tilde{\chi} + (f_i(\chi + R_i\tilde{\chi}) - f_i(\chi))(\chi - s) \\ + (g_i(\chi) - g_i(\chi + R_i\tilde{\chi})). $$

Let $0_{N-1}$ denote the vector in $\mathbb{R}^{N-1}$ with all zero entries. Using Assumption 4, we next show that the set $\{0_{N-1}\} \times A_0 \subset \mathbb{R}^{N+1}$ is locally asymptotically stable for (27), which corresponds to the local asymptotic stability of the set $A_x$ for system (23). Indeed, due to Assumption 4, the set $\{0_{N-1}\} \times A_0$ is locally asymptotically stable when the dynamics of (27) are restricted to the set $\{(\tilde{\chi}, \chi, s) \in \mathbb{R}^{N+1}: \tilde{\chi} = 0\}$. From [27, Theorem 2.24], it follows that $\{0_{N-1}\} \times A_0$ is a stable compact set for system (27). To show that the set is also attractive, let $\Omega_0 \subset \mathbb{R}^2$ be the region of attraction for the set $A_0$ of system (24), and let $\Omega_{\tilde{\chi}} \subset \mathbb{R}^{N-1}$ be an arbitrary compact set containing the origin. Due to stability of $\{0_{N-1}\} \times A_0 \subset \mathbb{R}^{N+1}$, there exists a set $\Omega_x \subset \mathbb{R}^{N+1}$ such that for all initial conditions in $\Omega_x$, the solutions stay within the set $\Omega_{\tilde{\chi}} \times \Omega_0$. With all solutions bounded, and $\tilde{\chi}(t)$ converging to $\{0_{N-1}\}$ as $t \to \infty$, it follows that $(\tilde{\chi}(t), \chi(t), s(t))$ converges to the largest invariant set contained in $\{0_{N-1}\} \times \Omega_0$ which is the set $\{0_{N-1}\} \times A_0$. The existence of a continuously differentiable Lyapunov function follows from classical converse arguments. Finally, it is seen from (26) that the behavior of (23) on $A_x$ is governed by the behavior of (24) on $A_0$, which is a periodic motion. Therefore, $A_x$ is a stable limit cycle. $\square$

**Remark 3.** In Lemma 1, we did not provide an explicit construction of the Lyapunov function $V_x$ for the QSS system (23). However, in the proof we saw that the asymptotic stability of the set $A_x$ is essentially established by using a cascade argument between $(\chi, s)$ dynamics and $\tilde{\chi}$ dynamics. In some instances, one can explicitly provide a construction of Lyapunov function $V_x$ from the individual Lyapunov functions associated with $(\chi, s)$-dynamics with respect to $A_0$, and $\tilde{\chi}$ dynamics with respect to $\{0_{N-1}\}$. In our conference paper [28], we provided such a construction under the assumption that $A_0$ is locally exponentially stable for system (24). $\triangleleft$

### D. Boundary-layer system

The boundary-layer system associated to (22) is

$$ \begin{aligned} \omega_i^+ &= \omega_i, & \text{for } i = 1, \dots, N \\ s^+ &= s \\ \zeta^+ &= \Lambda_p \zeta \\ \mathbf{p}^+ &\in G_p(\mathbf{p}). \end{aligned} \quad (28) $$

We note that $\zeta$-dynamics are completely decoupled from $(\omega_i, s)$ for each $i = 1, \dots, N$. Assumption 2b, with $h \equiv 0$, therefore corresponds to checking the stability of $(\zeta, \mathbf{p})$ dynamics with respect to the set $\{0_{N-1}\} \times \mathcal{P}$. This follows from the following lemma.

**Lemma 2.** Under Assumption 3, Assumption 2b holds for the boundary-layer system (28) with

$$ V_y(\omega, s, \zeta, \mathbf{p}) = \zeta^\top \zeta $$

and with $\gamma_m = \max_{s \in S} \bar{\lambda}_s^2$, where $\bar{\lambda}_s$ denotes the largest singular value of the matrix $\bar{\Lambda}_s := \Lambda_{s|s|} \cdots \Lambda_{s_1}$.

*Proof.* Since $\{G_q\}_{q \in S}$ is jointly connected and the matrix $W_q$ is associated with the graph $G_q$, it follows that

$$ W_{s|s|} + \cdots + W_{s_1} $$

is a primitive matrix. On the other hand, by [12, Lemma 2], it holds with $c_s > 0$ that

$$ \bar{W}_s := W_{s|s|} \cdots W_{s_1} \succcurlyeq c_s (W_{s|s|} + \cdots + W_{s_1}), \quad (29) $$

where $\succcurlyeq$ implies component-wise inequality. This implies $\bar{W}_s$ is primitive (because $A \succcurlyeq B$ with a primitive matrix $B$ implies that $A$ is primitive). Moreover, we have

$$ \begin{aligned} Q\bar{W}_s Q^{-1} &= QW_{s|s|} Q^{-1} \cdots QW_{s_2} Q^{-1} QW_{s_1} Q^{-1} \\ &= \begin{bmatrix} 1 & 0 \\ 0 & \Lambda_{s|s|} \cdots \Lambda_{s_1} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & \bar{\Lambda}_s \end{bmatrix}. \end{aligned} $$

Now we note that $\bar{W}_s$ is primitive and doubly stochastic whose diagonal elements are all positive. This implies that $\lim_{k \to \infty} \bar{W}_s^k = \mathbf{1}_N c_T^T$ with a vector $c \in \mathbb{R}^N$ [3]. Since the matrix $\mathbf{1}_N c_T^T$ has rank one, the matrix $\lim_{k \to \infty} \bar{\Lambda}_s^k$ should be the zero matrix; that is, $\|\bar{\Lambda}_s\|_2 < 1$ where $\| \cdot \|_2$ represents the induced-two-norm, and thus, $\gamma_m < 1$. This in turn implies that $\zeta^\top \bar{\Lambda}_s^\top \bar{\Lambda}_s \zeta \le \gamma_m \zeta^\top \zeta$, for all $\bar{s} \in S$, and Assumption 2b holds. $\square$

### E. Stability result

Combining the statements of Lemma 1 and Lemma 2, we arrive at the following result about synchronization of trajectories modeled by (22).

**Theorem 4.** Suppose that Assumptions 3 and 4 hold and consider a compact set $I_x \subset \Omega_x$ and a compact set $I_y \subset \mathbb{R}^{N-1}$, for which $(\omega(0, 0), (1/N) \sum_{i=1}^N \xi_i(0, 0)) \in I_x$ and $R^\top \xi(0, 0) \in I_y$. Then, there are $\epsilon^* > 0$, class $KL$ functions $\beta_x$ and $\beta_y$, and class $K$ functions $\chi_x$ and $\chi_y$ such that

$$ |(\omega(t, j), s(t, j))|_{A_x} \le \beta_x(|(\omega(0, 0), s(0, 0))|_{A_x}, t) + \chi_x(\epsilon) \quad (30) $$

$$ |\zeta(t, j)| \le \beta_y(|\zeta(0, 0)|, t/\epsilon) + \chi_y(\epsilon), \quad (31) $$

for all $\epsilon \in (0, \epsilon^*]$ and $(t, j) \in \text{dom}(\omega, s, \zeta)$. In particular, for each $\rho > 0$, there is $\epsilon^*$ such that, for each $\epsilon \in (0, \epsilon^*]$, the solutions to the multi-agent system (17) with (18) satisfy, for all $i, j = 1, \dots, N$,

$$ \limsup_{t \to \infty} |(\omega_i(t, j), \xi_i(t, j)) - (\omega_j(t, j), \xi_j(t, j))| \le \rho \quad (32) $$

$$ \limsup_{t \to \infty} |(\omega_i(t, j), \xi_i(t, j))|_{A_0} \le \rho \quad (33) $$

if $\text{dom}(\omega, \xi)$ is unbounded in t-axis.

*Proof.* Proof of (30) and (31) is a consequence of Theorem 3 by taking $I_{xy} = I_x \times I_y$. To show (32), we note that $\xi = \mathbf{1}_N s + R\zeta$. For any $i$ and $j$, we have that $|\xi_i - \xi_j| \le |\xi_i - s| + |s - \xi_j| = |R_i\zeta| + |R_j\zeta|$. This in turn implies from (31) that $\limsup_{t \to \infty} |(\xi_i(t, j) - \xi_j(t, j))|$ can be made arbitrarily small by taking $\epsilon^*$ sufficiently small. Since $\limsup_{t \to \infty} |(\omega_1, \dots, \omega_N, s)|_{A_x}$ can be made arbitrarily small by taking sufficiently small $\epsilon^*$ by (30), we infer that
$\limsup_{t \to \infty} |\omega_i(t, j) - \omega_j(t, j)|$ can also be made arbitrarily small due to the definition of the set $A_x$. This proves (32). The inequality in (33) is also a direct consequence of the foregoing analysis and the definition of the set $A_x$. □

Implication of (32) is that the coupled oscillator (17) practically synchronizes under the impulsive coupling (18). Moreover, (33) implies that the limit cycle $A_0$ is practically attractive for each oscillator under the coupling. In fact, the limit cycle $A_0$ is emergent and collective one for all the oscillators in the sense that it is the limit cycle of the blended dynamics (24) with (25). This can be seen from the illustrations reported in the following example.

**Remark 4.** It can be seen from the proof of Theorem 3 that, in the current case-study, the function $C_g$ is given by $C_g(s) = C_{\Lambda s}$ with $C_{\Lambda} := \max_{q \in Q} \|\text{blockdiag}(I_{N+1}, \Lambda_{q}^{n-1})\|$. Combined with the fact that the functions $V_y$ are quadratic, the inequality (31) takes the form

$$|\zeta(t, j)| \le \sqrt{2} C_{\Lambda} \gamma_m^{\frac{\sigma_1}{2N_0 n} \frac{t}{\epsilon}} |\zeta(0, 0)| + \chi_y(\epsilon)$$

with some class $\mathcal{K}$ function $\chi_y$. Also, from the proof, we see that, $|\xi_i - \xi_j| \le |R_i \zeta| + |R_j \zeta| \le 2 \max_{1 \le i \le N} \|R_i\| \cdot |\zeta|$, which implies that the rate of synchronization for $\xi_i$ is much faster than that for $\omega_i$ when $\epsilon$ is small enough. This behavior can be seen in the bottom left and bottom right plots of Fig. 1. △

**Example 3.** We have simulated<fnref n="6" /> a network with four Van der Pol oscillators described by (17) with $f_i(\omega_i) = \mu_i(\omega_i^2 - 1)$ where $\mu_1 = -.1$, $\mu_2 = 1$, $\mu_3 = 2$, and $\mu_4 = 3$, and $g_i(\omega_i) = \omega_i$ for $i = 1, \dots, 4$. The oscillators are connected through an undirected cyclic (and hence connected) graph $G = (\mathcal{V}, \mathcal{E})$,

<figure><img src="image_2.png" /><figcaption>Fig. 1. Plot of $(\omega_i(t), \xi_i(t))$ for four agents. Numbers are time stamps in increasing order, and all trajectories circulate clock-wise. (Top left): Case of no coupling; Started from four different initial conditions (time stamps '0'), their motions are not synchronized, and the black agent (agent 1) does not have a stable limit cycle so that it converges to the origin. (Top right): Coupled with $\epsilon = 1$; Jumps are observed but synchrony is weak. (Bottom left): Coupled with $\epsilon = 0.1$. The overlapped numbers imply that they are synchronized to some extent. (Bottom right): Coupled with $\epsilon = 0.01$: approximate synchronization.</figcaption></figure>

where $\mathcal{V} = \{1, 2, 3, 4\}$ and $\mathcal{E} = \{(1, 2), (2, 3), (3, 4), (4, 1)\}$. To define the reset rule for updating the state $\xi$ of the oscillators according to (18), we split the graph into four subgraphs $G_q = (\mathcal{V}, \mathcal{E}_q)$, $q \in \mathcal{Q} = \{1, 2, 3, 4\}$, with $\mathcal{E}_1 = \{(1, 2)\}$, $\mathcal{E}_2 = \{(2, 3)\}$, $\mathcal{E}_3 = \{(3, 4)\}$, $\mathcal{E}_4 = \{(4, 1)\}$. The set $\mathcal{S}$ comprises strings of length $n \ge 4$ such that each string contains at least a single occurrence of the letters 1, 2, 3, and 4. The weight matrix $W_q$ is constructed as in Example 2. For the timer, we work with, $\sigma_1 = 0.5$, $\sigma_2 = 1.5$, $N_0 = 2$, and study different cases for $\epsilon \in \{1, 0.1, 0.01\}$. The state trajectories of all the agents are plotted in Fig. 1 for three different values of $\epsilon$. Our proposed reset rule only updates the state $\xi_i$ of the $i$-th agent whenever an edge connecting this agent with one of its neighbors is activated. Also, as we make $\epsilon$ smaller, we see that all the agents converge to the same limit cycle and the distance between their states gets smaller as well. △

## V. PROOF OF THEOREM 2

From the given compact set $I_{xy} \subset C_{xy}$, we consider the compact set $I_x := \{x \in C_x : \exists y \text{ s.t. } (x, y) \in I_{xy}\}$, and the continuously differentiable function $V_x$ obtained from Assumption 1a. We let

$$\mu := \max_{x \in I_x} V_x(x)$$

and define

$$K_x := \{x \in C_x : V_x(x) \le \mu + \delta_{\mu}\} \quad (34)$$

with a small $\delta_{\mu} > 0$.

With $K_x$ in (34), Assumption 2a yields $V_y$, and using $V_y$, define

$$\nu := \max_{(x,y) \in I_{xy}, \mathbf{p} \in \mathcal{P}} V_y(x, y, \mathbf{p}).$$

Then, it is obvious that

$$\begin{align*} I_{xy} &\subset \{(x, y) \in C_{xy} : x \in I_x, V_y(x, y, \mathbf{p}) \le \nu, \forall \mathbf{p} \in \mathcal{P}\} \\ &\subset \{(x, y) \in C_{xy} : x \in K_x, \exists \mathbf{p} \in \mathcal{P} \text{ s.t. } V_y(x, y, \mathbf{p}) \le \nu + \delta_{\nu}\} \\ &:= K_{xy} \end{align*}$$

with a small $\delta_{\nu} > 0$.

We now use the compact set $K_{xy}$ to get bounds on the vector fields. In this regard, we note that, due to continuity of $f_x$, $f_y$ and $\nabla_x V_x$ on the compact set $K_{xy} \times [0, N_0] \times \mathcal{P}$, there are $M_x$, $M_y$, and $M_V$ such that

$$|f_x(x, y, \epsilon)| \le M_x, \quad |f_y(x, y, \epsilon)| \le M_y, \quad \forall \epsilon \in [0, 1],$$

and

$$|\nabla_x V_x(x)| \le M_V.$$

Let $\tilde{f}_x(x, y, \epsilon) := f_x(x, y, \epsilon) - f_x(x, h(x), 0)$. Then, there exists a class $\mathcal{K}_{\infty}$ function $\lambda_f$ such that<fnref n="7" />

$$|\tilde{f}_x(x, y, \epsilon)| \le \lambda_f(|y - h(x)|) + \lambda_f(\epsilon), \quad \forall \epsilon \in [0, 1], \quad (35)$$

$$|f(x, y)| \le F(|x|, |y|) \le F(|x|, |x|) + F(|y|, |y|) \le \lambda(|x|) + \lambda(|y|).$$

<footnote n="6">⁶Source code for the simulation in Example 3 is available at https://github.com/aneeltanwani/coupled_osc_switching.git.</footnote>

<footnote n="7">⁷For a compact set $C$ containing $(0, 0)$ and a continuous function, $C \ni (x, y) \mapsto f(x, y)$, satisfying $f(0, 0) = 0$, consider the map $F(a, b) := \max_{|x| \le a, |y| \le b} |f(x, y)|$. Then, $\exists$ class $\mathcal{K}_{\infty}$ function $\lambda$ s.t. $\forall (x, y) \in C$,</footnote>
on the compact set $K_{xy}$. On the other hand, by (BL-C1), $V_y$ is continuously differentiable and $V_y(x, y, \mathbf{p}) = 0$ when $y = h(x)$, that is, $V_y$ attains a minimum at $(x, h(x), \mathbf{p})$ for each $x \in K_x$ and $\mathbf{p} \in \mathcal{P}$. As a result, $\nabla_x V_y$ and $\nabla_y V_y$ are continuous and take the value zero at $(x, h(x), \mathbf{p})$, for each $x \in K_x$ and $\mathbf{p} \in \mathcal{P}$. We can therefore upper bound $|\nabla_x V_y|$ and $|\nabla_y V_y|$ by continuous functions which are non-decreasing with respect to $|y - h(x)|$. In other words, there exists a class $\mathcal{K}$ function $\lambda_{yy}$ such that

$$ 
\begin{aligned}
|\nabla_y V_y(x, y, \mathbf{p})| &\leq \lambda_{yy}(|y - h(x)|) \\
|\nabla_x V_y(x, y, \mathbf{p})| &\leq \lambda_{yy}(|y - h(x)|)
\end{aligned}
\quad (36) $$

on the compact set $K_{xy} \times [0, N_0] \times \mathcal{P}$.

### A. Construction of the Lyapunov function

The next step in the proof is to construct a Lyapunov function for the system (3) using the function $V_x$ from Assumption 1a and $V_y$ stipulated in Assumption 2a, and it will be shown that this function is decreasing on a compact set.

Towards this end, let us consider the following Lyapunov function candidate:

$$ W(x, y, \tau, \mathbf{p}) := \frac{e^{c_x \epsilon \tau} V_x(x)}{\delta_\mu + \mu - V_x(x)} + \frac{e^{-c_y \tau} V_y(x, y, \mathbf{p})}{\delta_\nu + \nu - V_y(x, y, \mathbf{p})}, \quad (37) $$

where $c_x$ and $c_y$ are constants such that

$$ 0 < c_x \le 1, \quad \text{and} \quad 0 < c_y \le 1, \quad (38) $$

and the explicit bounds on these constants are provided later in the proof in (45) and (49). We recall that the parameter $\epsilon$ belongs to the interval $[0, 1]$. Now, let

$$ \Omega := \left\{ (x, y, \tau, \mathbf{p}) \in C_{xy} \times [0, N_0] \times \mathcal{P} : W(x, y, \tau, \mathbf{p}) \le e^{N_0 \frac{\mu}{\delta_\mu}} + \frac{\nu}{\delta_\nu} := w^* \right\}. $$

**Claim 1:** It holds that

$$ I_{xy} \times [0, N_0] \times \mathcal{P} \subset \Omega \quad \text{and} \quad (39a) $$

$$ \Omega \subset K_{xy} \times [0, N_0] \times \mathcal{P}. \quad (39b) $$

**Proof of Claim 1:** Pick $(x, y, \tau, \mathbf{p}) \in I_{xy} \times [0, N_0] \times \mathcal{P}$. It is clear that $V_x(x) \le \mu$ and $V_y(x, y, \mathbf{p}) \le \nu$ by construction. Thus, $W(x, y, \tau, \mathbf{p}) \le e^{N_0} \mu / \delta_\mu + \nu / \delta_\nu = w^*$ by (38), and hence (39a) holds.

To prove (39b), pick $(x, y, \tau, \mathbf{p}) \in \Omega$. Then, we have

$$ \frac{V_x(x)}{\delta_\mu + \mu - V_x(x)} \le e^{-c_x \epsilon \tau} w^* \le w^*, \quad (40a) $$

which in turn implies that

$$ V_x(x) \le \frac{w^*}{1 + w^*} (\delta_\mu + \mu). \quad (40b) $$

Thus, $x \in K_x$. We also have

$$ \frac{V_y(x, y, \mathbf{p})}{\delta_\nu + \nu - V_y(x, y, \mathbf{p})} \le e^{c_y \tau} w^* \le e^{N_0} w^*, \quad (41a) $$

which implies that

$$ V_y(x, y, \mathbf{p}) \le \frac{e^{N_0} w^*}{1 + e^{N_0} w^*} (\delta_\nu + \nu). \quad (41b) $$

Therefore, $(x, y) \in K_{xy}$ by the definitions of $K_{xy}$.

It is seen from the above claim that the set of initial conditions for $(x, y, \tau, \mathbf{p})$ given by $I_{xy} \times [0, N_0] \times \mathcal{P}$ is contained in $\Omega$, and with our choice of $K_{xy}$, all level sets of $W$ for a given $w^*$ reside within $K_{xy} \times [0, N_0] \times \mathcal{P}$ over which Assumption 1a and Assumption 2a are valid. In what follows, we denote the derivative of $W$ along the flow dynamics of (3) by $\dot{W}$, and the change in the value of $W$ at jump instants by $\Delta W$, that is, $\Delta W := W(x^+, y^+, \tau^+, \mathbf{p}^+) - W(x, y, \tau, \mathbf{p})$.

From (40) and (41), we also arrive at the following inequalities which hold on the set $\Omega$:

$$ \underline{M}_\mu := \frac{1}{\delta_\mu + \mu} \le \frac{1}{\delta_\mu + \mu - V_x(x)} \le \frac{1 + w^*}{\delta_\mu + \mu} := \bar{M}_\mu, \quad (42) $$

$$ \underline{M}_\nu := \frac{1}{\delta_\nu + \nu} \le \frac{1}{\delta_\nu + \nu - V_y(x, y, \mathbf{p})} \le \frac{1 + e^{N_0} w^*}{\delta_\nu + \nu} =: \bar{M}_\nu. \quad (43) $$

Next, pick $\rho$ such that $0 < \rho < w^*$, and define

$$ m := \min \left\{ \frac{M_\mu}{2} \alpha_x \left( \bar{\alpha}_x^{-1} \left( \frac{\rho}{2e\bar{M}_\mu} \right) \right), \alpha_y \left( \bar{\alpha}_y^{-1} \left( \frac{\rho}{2\bar{M}_\nu} \right) \right) \right\}. \quad (44) $$

Take $c_y$ such that

$$ 0 < c_y \le \min \left\{ \ln \left( \frac{1}{2\gamma_c} \right), 1 \right\}. \quad (45) $$

**Claim 2:** There are sufficiently small positive $c_x \le 1/N_0$ and $\epsilon^*$ such that, for every $\epsilon \in (0, \epsilon^*]$, time derivative of $W$, defined in (37), along the flow dynamics (3a) satisfies

$$ \dot{W} \le -\frac{M_\mu}{2} \alpha_x(|x|_{A_x}) - \underline{\alpha}_y(|y - h(x)|) + m \quad (46) $$

on the set $\Omega$.

**Proof of Claim 2:** On the set $\Omega$, we have along the flow that, by (9), (10), (42), (43), (BL-C1), and Claim 1,

$$ 
\begin{aligned}
\dot{W} = & \frac{c_x \epsilon \dot{\tau} e^{c_x \epsilon \tau} V_x}{\delta_\mu + \mu - V_x} - \frac{c_y \dot{\tau} e^{-c_y \tau} V_y}{\delta_\nu + \nu - V_y} \\
& + \frac{e^{c_x \epsilon \tau} (\delta_\mu + \mu)}{(\delta_\mu + \mu - V_x)^2} \nabla_x V_x \cdot (f_x(x, h(x), 0) + \tilde{f}_x(x, y, \epsilon)) \\
& + \frac{e^{-c_y \tau} (\delta_\nu + \nu)}{(\delta_\nu + \nu - V_y)^2} (\nabla_y V_y \cdot f_y(x, y, \epsilon) + \nabla_x V_y \cdot f_x(x, y, \epsilon)) \\
\le & c_x \sigma_2 e^{c_x \epsilon N_0} \frac{1 + w^*}{\delta_\mu + \mu} \bar{\alpha}_x(|x|_{A_x}) - \frac{c_y \sigma_1}{\epsilon} e^{-c_y N_0} \frac{\alpha_y(|y - h(x)|)}{\delta_\nu + \nu} \\
& - e^0 \frac{\delta_\mu + \mu}{(\delta_\mu + \mu)^2} \alpha_x(|x|_{A_x}) \\
& + e^{c_x \epsilon N_0} \frac{(\delta_\mu + \mu)(1 + w^*)^2}{(\delta_\mu + \mu)^2} M_V(\lambda_f(|y - h(x)|) + \lambda_f(\epsilon)) \\
& + e^0 \frac{(1 + e^{N_0} w^*)^2}{(\delta_\nu + \nu)} (|\nabla_y V_y(x, y, \mathbf{p})| M_y + |\nabla_x V_y(x, y, \mathbf{p})| M_x).
\end{aligned}
$$

With the aforementioned bounds and $M_c := c_y \sigma_1 e^{-c_y N_0}$, under the assumption $c_x \le 1/N_0$, the above inequality continues as

$$ 
\begin{aligned}
\dot{W} \le & c_x \sigma_2 e \bar{M}_\mu \bar{\alpha}_x (|x|_{A_x}) - (M_c \bar{M}_\nu / \epsilon) \underline{\alpha}_y (|y - h(x)|) \\
& - \bar{M}_\mu \alpha_x (|x|_{A_x}) + e (\bar{M}_\mu^2 / \bar{M}_\mu) M_V (\lambda_f (|y - h(x)|) + \lambda_f (\epsilon)) \\
& + (\bar{M}_\nu^2 / \bar{M}_\nu) (M_y + M_x) \lambda_{yy} (|y - h(x)|).
\end{aligned}
$$
For convenience, take

$$M_1 := e \frac{\overline{M}_{\mu}^2}{\overline{M}_{\mu}} M_V, \quad \text{and} \quad M_2 := \frac{\overline{M}_{\nu}^2}{\overline{M}_{\nu}} (M_y + M_x).$$

Now, pick a small number $o_x > 0$. Then, by Lemma 4 in Appendix A, there exists $v_x > 0$ such that it holds on the compact set $\Omega$ that

$$\bar{\alpha}_x(|x|_{\mathcal{A}_x}) \le v_x \alpha_x(|x|_{\mathcal{A}_x}) + o_x. \quad (47)$$

Now, take

$$o_z = \frac{m}{3(M_1 + M_2)}.$$

Then, again from Lemma 4 in Appendix A, there exists $v_z > 0$ such that, on $\Omega$,

$$\max\{\lambda_f(|y-h(x)|), \lambda_{yy}(|y-h(x)|)\} \le v_z \underline{\alpha}_y(|y-h(x)|) + o_z. \quad (48)$$

Putting altogether, we obtain that

$$\dot{W} \le (c_x \sigma_2 e \overline{M}_\mu v_x - \overline{M}_\mu) \alpha_x (|x|_{\mathcal{A}_x}) + c_x \sigma_2 e \overline{M}_\mu o_x + M_1 \lambda_f (\epsilon) - \left( \frac{M_c \overline{M}_\nu}{\epsilon} - (M_1 + M_2) v_z \right) \underline{\alpha}_y (|y-h(x)|) + (M_1 + M_2) o_z.$$

Therefore, if one picks $c_x$ and $\epsilon^*$ such that

$$c_x \le \min \left\{ \frac{1}{N_0}, \frac{\overline{M}_\mu}{2\sigma_2 e v_x \overline{M}_\mu}, \frac{m}{3\overline{M}_\mu \sigma_2 e o_x} \right\}, \quad (49)$$

$$\epsilon^* \le \min \left\{ 1, \frac{M_c \overline{M}_\nu}{1 + (M_1 + M_2)v_z}, \lambda_f^{-1} \left( \frac{m}{3M_1} \right) \right\}, \quad (50)$$

then the claim follows.

//

**Claim 3:** For $\epsilon \in (0, \epsilon^*]$, the inequalities $\dot{W} < 0$ and $\Delta W < 0$, hold over the set $\Omega \setminus \Omega_\rho$, where

$$\Omega_\rho = \{(x, y, \tau, \mathbf{p}) \in C_{xy} \times [0, N_0] \times \mathcal{P} : W(x, y, \tau, \mathbf{p}) \le \rho\}.$$

*Proof of Claim 3:* We first note that, on the set $\Omega \setminus \Omega_\rho$,

$$\rho < W \le e \overline{M}_\mu \bar{\alpha}_x (|x|_{\mathcal{A}_x}) + \overline{M}_\nu \bar{\alpha}_y (|y - h(x)|).$$

Then, using the fact that if $(\rho < a + b)$ then $(\rho/2 < a)$ or $(\rho/2 < b)$, it follows that, if $W(x, y, \tau, \mathbf{p}) > \rho$, then

$$|x|_{\mathcal{A}_x} > \bar{\alpha}_x^{-1} \left( \frac{\rho}{2e \overline{M}_\mu} \right) \quad \text{or} \quad |y - h(x)| > \bar{\alpha}_y^{-1} \left( \frac{\rho}{2 \overline{M}_\nu} \right).$$

Then, by (44) and (46), it follows that $\dot{W} < 0$.

To show that $\Delta W < 0$, we observe that, on the set $\Omega$,

$$W(x^+, y^+, \tau^+, \mathbf{p}^+) = e^{c_x \epsilon (\tau-1)} \frac{V_x(x)}{\delta_\mu + \mu - V_x(x)} + e^{-c_y (\tau-1)} \frac{V_y(x, g_y(x, y, \mathbf{p}), \mathbf{q})}{\delta_\nu + \nu - V_y(x, g_y(x, y, \mathbf{p}), \mathbf{q})}$$

where $q \in G_p(\mathbf{p})$. Since $V_y(x, g_y(x, y, \mathbf{p}), \mathbf{q}) \le \gamma_c V_y(x, y, \mathbf{p})$, and thus, $1/(\delta_\nu + \nu - V_y(x, g_y(x, y, \mathbf{p}), \mathbf{q})) \le 1/(\delta_\nu + \nu - V_y(x, y, \mathbf{p}))$, we have, with $\tilde{\gamma}_c := \max\{e^{-c_x \epsilon}, e^{c_y \gamma_c}\}$ which satisfies $0 < \tilde{\gamma}_c < 1$ by (45),

$$W(x^+, y^+, \tau^+, \mathbf{p}^+) \le e^{-c_x \epsilon} \left( e^{c_x \epsilon \tau} \frac{V_x(x)}{\delta_\mu + \mu - V_x(x)} \right) + e^{c_y \gamma_c} \left( e^{-c_y \tau} \frac{V_y(x, y, \mathbf{p})}{\delta_\nu + \nu - V_y(x, y, \mathbf{p})} \right) \le \tilde{\gamma}_c W(x, y, \tau, \mathbf{p})$$

from which the claim follows.

//

## B. Convergence rate analysis

We now derive the bounds on trajectories given in (7). Let $z := y - h(x)$, and the dynamics in (3) is appended by $z$-dynamics, where $\dot{z} = f_z(x, z, \epsilon) = f_y(x, z + h(x), \epsilon) - \frac{\partial h}{\partial x}(x) f_x(x, z + h(x), \epsilon)$, and $z^+ = \hat{g}_y(x, z, \mathbf{p}) := g_y(x, z + h(x), \mathbf{p}) - h(x)$. Since $\mathbf{p}$ is piecewise constant, we use the notation $\mathbf{p}_j := \mathbf{p}(t, j)$. Let $\hat{\Omega} := \{(x, z) : \exists \tau, \mathbf{p} \text{ s.t. } (x, z + h(x), \tau, \mathbf{p}) \in \hat{\Omega}\}$. For every $(x, z, \epsilon) \in \hat{\Omega} \times [0, 1]$, one can find $M_x$ and $M_z$ such that

$$|f_x(x, z + h(x), \epsilon)| \le M_x, \quad |f_z(x, z, \epsilon)| \le M_z. \quad (51)$$

Define $X := \text{col}(x, z)$ and $M_{xz} := \max\{M_x, M_z\}$. Since $g_y$ is continuous, the mapping $X \mapsto G(X, \mathbf{p}) := (x, \hat{g}_y(x, z, \mathbf{p}))$ is also continuous and, for all $\mathbf{p} \in \mathcal{P}$, it holds that $G(0, \mathbf{p}) = 0$ because $\hat{g}_y(x, 0, \mathbf{p}) = 0$. Consequently, there is a class $\mathcal{K}$ function $\lambda_G$ such that for all $X_1, X_2 \in \hat{\Omega}$, we have⁸

$$|G(X_2, \mathbf{p}) - G(X_1, \mathbf{p})| \le \lambda_G(|X_2 - X_1|).$$

Similarly, we let $L_V^x$ be the Lipschitz constant for $V_x$ on $\hat{\Omega}$ and $L_V^z$ be the Lipschitz constant for $V_z(x, z, \mathbf{p}) = V_y(x, z + h(x), \mathbf{p})$ on $\hat{\Omega}$ for all $\mathbf{p} \in \mathcal{P}$.

In what follows, we will use the symbol $t_j$ to denote the continuous time at which the $j$-th jump occurs, that is,

$$t_j := \min\{t \in \mathbb{R} \text{ such that } (t, j) \in \text{dom } X\}. \quad (52)$$

In other words, after the $j$-th jump, the system flows from time $(t_j, j)$ to $(t_{j+1}, j)$. Using the bounds in (51), combined with the fact that,

$$t_{j+1} - t_j \le \frac{N_0 \epsilon}{\sigma_1} \le \frac{N_0}{\sigma_1},$$

we get, for each $j \in \mathbb{N}$ and $t_j \le t \le t_{j+1}$,

$$|X(t, j) - X(t_j, j)| \le C_0 \epsilon \quad (53)$$

where $C_0 := \sqrt{2} M_{xz} (N_0 / \sigma_1)$. This leads to

$$|G(X(t_{j+1}, j), \mathbf{p}_j) - G(X(t_j, j), \mathbf{p}_j)| \le \lambda_G(C_0 \epsilon).$$

By letting $\hat{Y}_j := X(t_j, j) - G(X(t_{j-1}, j-1), \mathbf{p}_{j-1})$, we get,

$$\begin{aligned} |\hat{Y}_j| &= |X(t_j, j) - G(X(t_{j-1}, j-1), \mathbf{p}_{j-1})| \\ &= |G(X(t_j, j-1), \mathbf{p}_{j-1}) - G(X(t_{j-1}, j-1), \mathbf{p}_{j-1})| \\ &\le \lambda_G(C_0 \epsilon). \end{aligned} \quad (54)$$

We use this bound in the following inequalities:

$$\begin{aligned} V_z(X(t_j, j), \mathbf{p}_j) &= V_z(G(X(t_{j-1}, j-1), \mathbf{p}_{j-1}) + \hat{Y}_j, \mathbf{p}_j) \\ &\le V_z(G(X(t_{j-1}, j-1), \mathbf{p}_{j-1}), \mathbf{p}_j) + L_V^z |\hat{Y}_j| \\ &\le \gamma_c V_z(X(t_{j-1}, j-1), \mathbf{p}_{j-1}) + L_V^z \lambda_G(C_0 \epsilon) \\ &\vdots \end{aligned}$$

⁸For a continuous function $\hat{G}$ over a compact set $K$, we observe that

$$|\hat{G}(X_2) - \hat{G}(X_1)| \le \max_{|Y| \le |X_2 - X_1|, X_1 + Y \in K} |\hat{G}(X_1 + Y) - \hat{G}(X_1)| \\ =: \tilde{G}(|X_2 - X_1|)$$

where $\tilde{G} : \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$ is continuous and $\tilde{G}(0) = 0$. Therefore, there is a class $\mathcal{K}$ function $\lambda_G$ such that $\tilde{G}(s) \le \lambda_G(s)$.
$$ \leq \gamma_c^j V_z(X(0,0), \mathbf{p}_0) + \sum_{\ell=0}^{j-1} \gamma_c^\ell L_V^z \lambda_G(C_0\epsilon). $$

Therefore, for all $(t_j, j)$, it holds that,

$$ |z(t_j, j)| \leq \underline{\alpha}_y^{-1} \left( \gamma_c^j V_z(X(0,0), \mathbf{p}_0) + \frac{L_V^z \lambda_G(C_0 \epsilon)}{1 - \gamma_c} \right) . \quad (55) $$

Next, we note that, for $(t, j)$ satisfying $t_j \leq t \leq t_{j+1}$, we get

$$ \begin{aligned} |z(t,j)| &\leq |z(t,j) - z(t_j, j)| + |z(t_j, j)| \\ &\leq C_0\epsilon + |z(t_j, j)|. \end{aligned} $$

Substituting (55) in the last inequality, we immediately get,

$$ |z(t,j)| \leq \underline{\alpha}_y^{-1} \left( \gamma_c^{\frac{\sigma_1}{N_0} \frac{t}{\epsilon}} \bar{\alpha}_y(|z(0,0)|) + \frac{L_V^z \lambda_G(C_0 \epsilon)}{1 - \gamma_c} \right) + C_0 \epsilon, $$

which yields (13).

For deriving the inequality (12), we note that

$$ \begin{aligned} \dot{V}_x &= \nabla_x V_x \cdot f_x(x, h(x), 0) + \nabla_x V_x \cdot \tilde{f}_x(x, z + h(x), \epsilon) \\ &\leq -\alpha_x(V_x(x)) + L_V^x(\lambda_f(|z|) + \lambda_f(\epsilon)). \end{aligned} \quad (56) $$

The inequality (56) provides a Lyapunov characterization for ISS with $z$ and $\epsilon$ as inputs, but we need to analyze it carefully to get the estimate on $x$-trajectory independently of the initial conditions on $y$-dynamics. This is a direct consequence of Lemma 3 given in Section V-C.

### C. Derivation of $KL$ function

We use the following lemma to arrive at the inequality (7a) with the function $\beta_x$ given in (12).

**Lemma 3.** Consider system (3) and suppose that the function $V_x$ satisfies

$$ \dot{V}_x \leq -\alpha_x(V_x(x)) + \eta_z(|z|) + \eta_\epsilon(\epsilon) \quad (57) $$

for some $\alpha_x \in \mathcal{K}$, $\eta_z \in \mathcal{K}$, $\eta_\epsilon \in \mathcal{K}$, and it holds that

$$ |z(t,j)| \leq \beta_z(|z(0,0)|, t/\epsilon) + \chi_z(\epsilon) \quad (58) $$

for some $\beta_z \in \mathcal{KL}$ and $\chi_z \in \mathcal{K}$. Then there exist $\beta \in \mathcal{KL}$ and $\chi \in \mathcal{K}$ such that, for all $(t, j) \in \text{dom}(x, z)$,

$$ V_x(x(t,j)) \leq \beta(V_x(x(0,0)), t) + \chi(\epsilon). \quad (59) $$

**Proof.** We first observe that the substitution of (58) in (57) yields

$$ \dot{V}_x \leq -\alpha_x(V_x(x)) + \eta_z(2\beta_z(|z(0,0)|, t/\epsilon)) + \eta(\epsilon), \quad (60) $$

where $\eta(\epsilon) := \eta_z(2\chi_z(\epsilon)) + \eta_\epsilon(\epsilon)$. Using [2, Lemma IV.2 and Cor. IV.3], it can be shown that with the class $KL$ function $\beta_x^0$ in (11), it holds that<fnref n="9" />

$$ \begin{aligned} V_x(x(t,j)) \leq & \beta_x^0(V_x(x(0,0)), t) \\ & + \int_0^t \left( \eta_z \left( 2\beta_z \left( |z(0,0)|, \frac{s}{\epsilon} \right) \right) + \eta(\epsilon) \right) ds \end{aligned} $$

for $(t, j) \in \text{dom}(x, z)$. Recalling that, $\underline{\alpha}_y(|z(0,0)|) \leq V_z(x(0,0), z(0,0), \mathbf{p}) \leq \nu$ for all $\mathbf{p} \in \mathcal{P}$, we therefore have, for $0 \leq t \leq 1$,

$$ \begin{aligned} V_x(x(t,j)) \leq & \beta_x^0(V_x(x(0,0)), t) \\ & + \int_0^1 \eta_z \left( 2\beta_z \left( \underline{\alpha}_y^{-1}(\nu), \frac{s}{\epsilon} \right) \right) ds + \eta(\epsilon). \end{aligned} $$

Note that

$$ \begin{aligned} \int_0^1 \eta_z \left( 2\beta_z \left( \underline{\alpha}_y^{-1}(\nu), \frac{s}{\epsilon} \right) \right) ds &= \epsilon \int_0^1 \eta_z \left( \beta_z \left( \underline{\alpha}_y^{-1}(\nu), \tau \right) \right) d\tau \\ &\leq \chi_0(\epsilon), \end{aligned} $$

where $\chi_0$ is a class $K$ function defined on $[0, 1]$, whose existence follows from the fact that, for any class $KL$ function $\hat{\beta}$, we have $\lim_{\epsilon \to 0} \epsilon \int_0^{1/\epsilon} \hat{\beta}(\nu, \tau)d\tau = 0$. Hence, for $(t, j) \in \text{dom}(x, z)$ with $0 \leq t \leq 1$, we have

$$ V_x(x(t,j)) \leq \beta_x^0(V_x(x(0,0)), t) + \chi_0(\epsilon) + \eta(\epsilon). \quad (61) $$

On the other hand, define a class $K$ function

$$ \hat{\eta}(\epsilon) := \begin{cases} \eta_z \left( 2\beta_z \left( \underline{\alpha}_y^{-1}(\nu), \frac{1}{\epsilon} \right) \right), & \epsilon > 0, \\ 0, & \epsilon = 0. \end{cases} $$

Now, for $t \geq 1$, it follows from (60) that

$$ \dot{V}_x \leq -\alpha_x(V_x(x)) + \hat{\eta}(\epsilon) + \eta(\epsilon). $$

Therefore, there exist a class $K$ function $\chi_1$ such that, for $(t, j) \geq (1, j') \in \text{dom}(x, z)$,

$$ \begin{aligned} V_x(x(t,j)) &\leq \beta_x^0(V_x(x(1,j')), t-1) + \chi_1(\epsilon) \\ &\leq \beta_x^0(2\beta_x^0(V_x(x(0,0)), 1), t-1) \\ &\quad + \beta_x^0(2\chi_0(\epsilon) + 2\eta(\epsilon), 0) + \chi_1(\epsilon). \end{aligned} \quad (62) $$

We now combine the two inequalities (61) and (62) by defining

$$ \beta(r, s) := \max \left\{ \beta_x^0(r, s), \beta_x^0(2\beta_x^0(r, 1), \max\{0, s-1\}) \right\} $$

and $\chi(\epsilon) := \max\{\chi_0(\epsilon) + \eta(\epsilon), \beta_x^0(2\chi_0(\epsilon) + 2\eta(\epsilon), 0) + \chi_1(\epsilon)\}$, which yields

$$ V_x(x(t,j)) \leq \beta(V_x(x(0,0)), t) + \chi(\epsilon), \quad (t, j) \in \text{dom}(x, z). $$

This in turn implies (59). □

## VI. PROOF OF THEOREM 3

For the proof, it is convenient to introduce $z := y - h(x)$ and consider the closed set

$$ C_X := \{(x, z) : (x, z + h(x)) \in C_{xy}\}. $$

With the state $(x, z, \tau, \mathbf{p}) \in C_X \times [0, N_0] \times \mathcal{P}$, we can rewrite the flow dynamics of (3) as

$$ \begin{aligned} \dot{x} &= f_x(x, h(x), 0) + [f_x(x, z + h(x), \epsilon) - f_x(x, h(x), 0)] \\ &= f_x(x, h(x), 0) + \tilde{f}_x(x, z + h(x), \epsilon), \\ \dot{z} &= f_y(x, z + h(x), \epsilon) - \frac{\partial h}{\partial x}(x) f_x(x, z + h(x), \epsilon) \\ &=: f_z(x, z, \epsilon), \end{aligned} $$

$$ \epsilon \dot{\tau} \in [\sigma_1, \sigma_2], $$

$$ \dot{\mathbf{p}} = 0, $$

(63a)

<footnote n="9">⁹In [2, Lemma IV.2 and Cor. IV.3], the authors work with a positive definite function $\alpha_x$ which is not necessarily increasing and hence the resulting class $KL$ function is not the same as (11). However, in our case, $\alpha_x$ is a class $K$ function and we can modify the proof of [2, Lemma IV.2 and Cor. IV.3] to define $\beta_x^0$ as proposed in [18, Lemma 4.4] which coincides with (11).</footnote>
and the jump dynamics of (3), for $(x, z, \tau, \mathbf{p}) \in C_X \times [1, N_0] \times \mathcal{P}$, as

$$
\begin{align*}
x^+ &= x, \\
z^+ &= g_y(x, z + h(x), \mathbf{p}) - h(x) =: \tilde{g}(x, z, \mathbf{p}), \\
\tau^+ &= \tau - 1, \\
\mathbf{p}^+ &\in G_{\mathbf{p}}(\mathbf{p}).
\end{align*}
\tag{63b}
$$

For simplicity, define $X := \text{col}(x, z)$ so that parts of (63) can be written concisely as

$$
\begin{align*}
\dot{X} &= F(X, \epsilon) := \begin{bmatrix} f_x(x, z + h(x), \epsilon) \\ f_z(x, z, \epsilon) \end{bmatrix} \\
X^+ &= G(X, \mathbf{p}) := \begin{bmatrix} x \\ \tilde{g}(x, z, \mathbf{p}) \end{bmatrix}.
\end{align*}
$$

Let $P \subset \mathcal{P}$ be the set introduced in Assumption 2b. Without loss of generality, we assume that $p(0, 0) \in \mathcal{P}$. For a solution of (63), we can extract a sequence from the values of p-trajectory which corresponds to the change in values of p and also when it enters the set $P$. Letting $p(0, 0) = p_1^0$, we label this sequence as $p_1^0, p_2^0, \dots, p_{r_0}^0, p_1^1, p_2^1, \dots, p_1^{k+1}, \dots$. In this sequence, $p_1^k \in P$ for each $k \in \mathbb{N}$ and after $r_k \le n$ transitions, we arrive at $p_1^{k+1} \in P$ so that (BL-M2) in Assumption 2b holds. We let $J_k = \sum_{j=0}^k r_j$. Consequently, we will use the notation $\tilde{g}_r^k(x, z) := \tilde{g}(x, z, p_r^k)$ and $G_r^k(X) := G(X, p_r^k)$ with $r \le r_k \le n$. For multiple jumps indexed by the discrete variable $\mathbf{p}$, we will use the notation $G_{r,1}^k(X) := G_r^k \circ G_{r-1}^k \circ \dots \circ G_1^k$.

Now we look for a compact set over which the trajectories
evolve for a while from the initial conditions. This allows
us to find quantitative measures used for ultimate bound and
convergence rate of $z$, with which we look at the asymptotic
behavior of the overall coupled system.

Recalling our switching model, it follows that, for any $j$,

$$
t_{j+1} - t_j \le \frac{N_0 \epsilon}{\sigma_1} \le \frac{N_0}{\sigma_1} =: T_0. \quad (64)
$$

Now, let $\phi_{\epsilon}(t, X_0)$ be the solution $X$ at time $t$ to the flow
$\dot{X} = F(X, \epsilon)$ with $X(0) = X_0$. Define

$$
\Phi_{\mathbf{p}}(A) := \{G(\phi_{\epsilon}(t, X_0), \mathbf{p}) \in C_X : X_0 \in A, \epsilon \in [0, 1], t \in [0, T_0]\},
$$

which is a set of reachable states from a set $A$ within the time
$T_0$ when $\epsilon$ belongs to $[0, 1]$, followed by one jump through
$G(\cdot, \mathbf{p})$. Note that $\Phi_{\mathbf{p}}(A)$ is compact if $A$ is compact. Let $\mathcal{P}^*$
be the set of all sequences of elements of $\mathcal{P}$ of length at most
$n$ for which (BL-M2) holds, and denote a generic element of
$\mathcal{P}^*$ by $\mathbf{p}^k = \mathbf{p}_1^k \cdot \mathbf{p}_2^k \cdots \mathbf{p}_{r_k}^k$. Define, for $\mathbf{p}^k \in \mathcal{P}^*$

$$
\Phi_{\mathbf{p}^k}(A) := \Phi_{\mathbf{p}_1^k} \circ \dots \circ \Phi_{\mathbf{p}_k^k}(A),
$$

which is compact when $A$ is compact. Finally, let

$$
\Phi_{\mathcal{P}}(A) := \cup_{\mathbf{p}^k \in \mathcal{P}^*} \Phi_{\mathbf{p}^k}(A),
$$

which is a reachable set from $A$ under the jumps by any
sequence in $\mathcal{P}^*$ and the flows up to the time $T_0$ with
$\epsilon \in [0, 1]$. Since $\mathcal{P}^*$ is finite, the set $\Phi_{\mathcal{P}}(A)$ is compact
whenever $A$ is compact.

Now, with $K_x$ in (34), Assumption 2b yields $V_y$, and using
$V_y$, define

$$
\nu := \max_{(x,y) \in I_{xy}, \mathbf{p} \in \mathcal{P}} V_y(x, y, \mathbf{p}). \quad (65)
$$

Then, it is obvious that

$$
\begin{align*}
I_{xy} &\subset \{(x, y) \in C_{xy} : x \in I_x, V_y(x, y, \mathbf{p}) \le \nu, \forall \mathbf{p} \in \mathcal{P}\} \\
&\subset \{(x, y) \in C_{xy} : x \in K_x, \exists \mathbf{p} \in \mathcal{P} \text{ s.t. } V_y(x, y, \mathbf{p}) \le \nu\} \\
&=: K_{xy0}.
\end{align*}
$$

Let $K_X := \{(x, z) : (x, z + h(x)) \in K_{xy0}\}$. For convenience, let $V_z(x, z, \mathbf{p}) := V_y(x, z + h(x), \mathbf{p})$, $\mathbf{p} \in \mathcal{P}$. Then,

$$
K_X = \cup_{\mathbf{p} \in \mathcal{P}} \{(x, z) \in C_X : x \in K_x, V_z(x, z, \mathbf{p}) \le \nu\}. \quad (66)
$$

Finally, define

$$
\Omega_X := \Phi_P(K_X).
$$

And we keep in mind that, by the construction of $\Omega_X$,

$$
X(t, j) \in \Omega_X, \quad \forall (t_{J_{k-1}}, J_{k-1}) \le (t, j) \le (t_{J_k}, J_k), \quad (67)
$$

whenever $X(t_{J_{k-1}}, J_{k-1}) \in K_X$ for any $k \ge 1$.<fnref n="10" />

On the compact set $\Omega_X \times [0, 1]$ for $(x, z, \epsilon)$, one can find $M_x$ and $M_z$ such that

$$
|f_x(x, z + h(x), \epsilon)| \le M_x, \quad |f_z(x, z, \epsilon)| \le M_z. \quad (68)
$$

Define $M_{xz} := \max\{M_x, M_z\}$. Due to continuity of $\tilde{f}$, there is a class $K_\infty$ function $\lambda_f$ such that

$$
|\tilde{f}_x(x, z + h(x), \epsilon)| \le \lambda_f(|z|) + \lambda_f(\epsilon), \quad (69)
$$

for all $(x, z, \epsilon) \in \Omega_X \times [0, 1]$. By construction, the mapping $X \mapsto G(X, \mathbf{p})$ is continuous and $G(0, \mathbf{p}) = 0$, for all $\mathbf{p} \in \mathcal{P}$. Over the compact set $\Omega_X$, we can find a class $\mathcal{K}$ function $\lambda_G$ such that $|G(X_2, \mathbf{p}) - G(X_1, \mathbf{p})| \le \lambda_G(|X_2 - X_1|)$, for all $X_1, X_2 \in \Omega_X$ and all $\mathbf{p} \in \mathcal{P}$. We will use the notation $\lambda_G^1 = \lambda_G$ and $\lambda_G^\ell := \lambda_G \circ \lambda_G^{\ell-1}$, for each $\ell \ge 2$. Also, we let $L_V^x$ be Lipschitz constant for $V_x$ on $\Omega_X$ and $L_V^z$ be Lipschitz constant for $V_z(\cdot, \cdot, \mathbf{p})$ on $\Omega_X$ for all $\mathbf{p} \in \mathcal{P}$.

Suppose that $\epsilon$ is sufficiently small. In particular, $\epsilon$ is small enough to satisfy the inequalities (71) given below, for which we define:

$$
\nu^* := \lambda_f^{-1} \left( \frac{\alpha_x(\mu + \delta_\mu)}{2L_V^x} \right),
$$

$$
C_0 := \sqrt{2M_{xz}}(N_0/\sigma_1), \quad C_n(\epsilon) := \sum_{\ell=1}^{n} \lambda_G^{\ell}(C_0\epsilon),
$$

$$
C_a(\epsilon) := C_0\epsilon + C_n(\epsilon), \quad C_g(s) := \max\{s, \lambda_G^{n-1}(s)\},
$$

and choose a positive integer $N$ such that

$$
\gamma_m^N \nu \le \frac{1}{2} \alpha_y (C_g^{-1} (\nu^*/2)). \quad (70)
$$

Then, the inequalities for $\epsilon$ are:

$$
0 < \epsilon \le 1, \tag{71a}
$$

$$
L_V^z C_n(\epsilon) \le (1 - \gamma_m)\nu, \quad (71b)
$$

<footnote n="10">¹⁰In the proof, the time inequality such as $(t_a, a) \le (t, j) < (t_b, b)$ should be understood as $t_a \le t < t_b$, $a \le j < b$, and $(t, j) \in \text{dom } X$, even if $\text{dom } X$ is often omitted.</footnote>
$$L_V^x M_x n \frac{N_0}{\sigma_1} \epsilon \le \frac{\delta_\mu}{N}, \qquad (71c)$$

$$L_V^x \lambda_f(\epsilon) < \frac{\alpha_x(\mu + \delta_\mu)}{2}, \qquad (71d)$$

$$\frac{L_V^z C_n(\epsilon)}{1 - \gamma_m} \le \frac{1}{2} \alpha_y (C_g^{-1}(\nu^*/2)), \qquad (71e)$$

$$C_a(\epsilon) \le \frac{\nu^*}{2}. \qquad (71f)$$

With such $\epsilon$, our first claim is the following.

**Claim 1:** For any $k \ge 1$, suppose that $V_x(x(t_{J_{k-1}}, J_{k-1})) \le \mu + \delta_\mu$ and $X(t_{J_{k-1}}, J_{k-1}) \in K_X$. If $V_x(x(t_{J_k}, J_k)) \le \mu + \delta_\mu$, then $X(t_{J_k}, J_k) \in K_X$.

**Proof of Claim 1:** By the construction of $\Omega_X$, (67) holds, and so, with (64), the evolution of (63) satisfies

$$|X(t, j) - X(t_j, j)| \le C_0 \epsilon \qquad (72)$$

in the time interval $t_j \le t \le t_{j+1}$ for each $j$ such that $J_{k-1} \le j < J_k$. Now, with

$$\tilde{Y}_k := X(t_{J_k}, J_k) - G_{r_{k,1}}^{k-1}(X(t_{J_{k-1}}, J_{k-1})),$$

we claim that

$$|\tilde{Y}_k| = |X(t_{J_k}, J_k) - G_{r_{k,1}}^{k-1}(X(t_{J_{k-1}}, J_{k-1}))| \le \sum_{\ell=1}^{r_k} \lambda_G^\ell(C_0\epsilon) \le C_n(\epsilon). \qquad (73)$$

The last inequality can be shown by induction. To see the underlying argument with notational simplicity, let us consider the case when $k = 1$. Note that

$$|X(t_1, 1) - G_1^0(X(0, 0))| = |G_1^0(X(t_1, 0)) - G_1^0(X(0, 0))| \le \lambda_G(C_0\epsilon),$$

and with $r_1 \ge 2$,

$$|X(t_2, 2) - G_{2,1}^0(X(0, 0))| \le |G_2^0(X(t_2, 1)) - G_2^0(X(t_1, 1))| + |G_2^0(X(t_1, 1)) - G_{2,1}^0(X(0, 0))| = |G_2^0(X(t_2, 1)) - G_2^0(X(t_1, 1))| + |G_{2,1}^0(X(t_1, 0)) - G_{2,1}^0(X(0, 0))| \le \lambda_G(C_0\epsilon) + \lambda_G^2(C_0\epsilon).$$

Repeating the arguments, we obtain that

$$|X(t_{J_1}, J_1) - G_{r_1,1}^0(X(0, 0))| \le (\lambda_G + \lambda_G^2 + \dots + \lambda_G^{r_1})(C_0\epsilon) \qquad (74)$$

and (73) can be proved similarly. With (73) at hand, we see from Assumption 2b, (65), and (71b), that

$$\begin{align*} V_z(X(t_{J_k}, J_k), \mathfrak{p}_1^k) & = V_z(G_{r_{k,1}}^{k-1}(X(t_{J_{k-1}}, J_{k-1})) + \tilde{Y}_k, \mathfrak{p}_1^k) \\ & \le V_z(G_{r_{k,1}}^{k-1}(X(t_{J_{k-1}}, J_{k-1})), \mathfrak{p}_1^k) + L_V^z |\tilde{Y}_k| \quad (75) \\ & \le \gamma_m V_z(X(t_{J_{k-1}}, J_{k-1}), \mathfrak{p}_1^{k-1}) + L_V^z |\tilde{Y}_k| \\ & \le \gamma_m \nu + L_V^z C_n(\epsilon) \le \nu. \end{align*}$$

Therefore, considering the definition of $K_x$ in (34) and $K_X$ in (66), we conclude that $X(t_{J_k}, J_k) \in K_X$ under the assumption $V_x(x(t_{J_k}, J_k)) \le \mu + \delta_\mu$. //

**Claim 2:** With $N$ chosen in (70),

$$V_x(x(t, j)) \le \mu + \delta_\mu, \quad \forall (0, 0) \le (t, j) \le (t_{J_N}, J_N), \qquad (76)$$

and for $k = 0, \dots, N$,

$$|z(t_{J_k}, J_k)| \le \alpha_y^{-1} \left( \gamma_m^k \nu + \frac{L_{V_z} C_n(\epsilon)}{1 - \gamma_m} \right) \quad \text{and} \qquad (77)$$

$$|z(t, j)| \le C_a(\epsilon) + C_g(|z(t_{J_k}, J_k)|), \qquad (78) \\ t_{J_k} \le t \le t_{J_{k+1}}, J_k \le j < J_{k+1}.$$

**Proof of Claim 2:** We first note that, as long as $(x, z, \epsilon) \in \Omega_X \times [0, 1]$,

$$|\dot{V}_x| = \left| \frac{\partial V_x}{\partial x}(x) f_x(x, z + h(x), \epsilon) \right| \le L_V^x M_x. \qquad (79)$$

Therefore, it follows from $V_x(x(0, 0)) \le \mu$ and $X(0, 0) \in K_X$, and from (71c) and (64), that $V_x(x(t, j)) \le \mu + \delta_\mu/N$ for $(0, 0) \le (t, j) \le (t_{J_1}, J_1)$. Moreover, from Claim 1, we also have that $X(t_{J_1}, J_1) \in K_X$. Repeating this argument, we get

$$V_x(x(t, j)) \le \mu + \frac{k}{N}\delta_\mu, \quad \forall (0, 0) \le (t, j) \le (t_{J_k}, J_k),$$

which implies (76). Therefore, $X(t, j) \in \Omega_X$ for $(0, 0) \le (t, j) \le (t_{J_N}, J_N)$.

In order to prove (77), we proceed further with (75) as

$$\begin{align*} V_z(X(t_{J_k}, J_k), \mathfrak{p}_1^k) &= V_z(G_{r_{k,1}}^{k-1}(X(t_{J_{k-1}}, J_{k-1}))) + \tilde{Y}_k, \mathfrak{p}_1^k \\ &\le V_z(G_{r_{k,1}}^{k-1}(X(t_{J_{k-1}}, J_{k-1})), \mathfrak{p}_1^k) + L_V^z |\tilde{Y}_k| \\ &\le \gamma_m V_z(X(t_{J_{k-1}}, J_{k-1}), \mathfrak{p}_1^{k-1}) + L_V^z |\tilde{Y}_k| \\ &\quad + \gamma_m L_V^z |\tilde{Y}_{k-1}| + L_V^z |\tilde{Y}_k| \\ &\le \gamma_m^2 V_z(X(t_{J_{k-2}}, J_{k-2}), \mathfrak{p}_1^{k-2}) + \gamma_m L_V^z |\tilde{Y}_{k-1}| + L_V^z |\tilde{Y}_k| \\ &\vdots \\ &\le \gamma_m^{k-1} V_z(X(t_{J_1}, J_1), \mathfrak{p}_1^k) + \sum_{\ell=0}^{k-2} \gamma_m^\ell L_V^z |\tilde{Y}_{k-\ell}| \\ &\le \gamma_m^k V_z(X(0, 0), \mathfrak{p}_1^0) + \sum_{\ell=0}^{k-1} \gamma_m^\ell L_V^z |\tilde{Y}_{k-\ell}|. \end{align*}$$

This means that, since $\sum_{\ell=0}^\infty \gamma_m^\ell = 1/(1 - \gamma_m)$ and $|\tilde{Y}_k| \le C_n(\epsilon)$,

$$|z(t_{J_k}, J_k)| \le \alpha_y^{-1} (V_z(X(t_{J_k}, J_k), \mathfrak{p}_1^k)) \qquad (80)$$

$$\le \alpha_y^{-1} \left( \gamma_m^k V_z(X(0, 0), \mathfrak{p}_1^0) + \frac{L_V^z C_n(\epsilon)}{1 - \gamma_m} \right) \qquad (81)$$

$$\le \alpha_y^{-1} \left( \gamma_m^k \nu + \frac{L_{\tilde{V}}^z C_n(\epsilon)}{1 - \gamma_m} \right), \qquad (82)$$

which yields (77).

While (82) indicates the decrease of $z$ as $k$ increases, let us look at the behavior of $z$ where we may go outside the set $P$. For this, we note that for every $\mathfrak{p} \in \mathcal{P}$, it holds that $\tilde{g}(x, 0, \mathfrak{p}) = 0$ and $|\tilde{g}(x, z, \mathfrak{p})| \le \lambda_G(|z|)$ for all $(x, z) \in \Omega_X$. Thus, it follows that

$$|\tilde{g}_r^k \circ G_{r-1}^k \circ \dots \circ G_1^k(X)| \le \lambda_G^r(|z|), \quad 1 \le r \le r_k,$$
for $X = (x, z) \in K_X$ because, e.g., $|\tilde{g}_2^k(G_1^k(x, z))| \le \lambda_G(|\tilde{g}_1^k(x, z)|) \le \lambda_G^2(|z|)$. Then, for $k = 0, \dots, N$, if $X(t_{J_k}, J_k) \in K_X$, an upper bound of $|z(t, j)|$ for $t_j \le t \le t_{j+1}$, where $J_k \le j < J_{k+1}$, can be found as

$$ 
\begin{aligned}
|z(t,j)| &\le |z(t,j) - z(t_j, j)| + |z(t_j, j)| \\
&\le C_0\epsilon + |z(t_j, j) - \tilde{g}_j^k \circ G_{j-1}^k \circ \dots \circ G_1^k(X(t_{J_k}, J_k))| \\
&\phantom{\le} + \lambda_G^j(|z(t_{J_k}, J_k)|).
\end{aligned}
$$

An upper bound of $|z(t_j, j) - \tilde{g}_j^k \circ G_{j-1}^k \circ \dots \circ G_1^k(X(t_{J_k}, J_k))|$ can be obtained similar to (74) as $(\lambda_G + \lambda_G^2 + \dots + \lambda_G^j)(C_0\epsilon)$. Hence, for $t_j \le t \le t_{j+1}$ and $J_k \le j < J_{k+1}$,

$$ 
\begin{aligned}
|z(t,j)| &\le C_0\epsilon + \sum_{\ell=1}^{j} \lambda_G^\ell(C_0\epsilon) + \lambda_G^j(|z(t_{J_k}, J_k)|) \\
&\le C_a(\epsilon) + C_g(|z(t_{J_k}, J_k)|).
\end{aligned}
$$

Similar inequality holds for $t_{J_k} \le t \le t_{J_{k+1}}$ and $J_k \le j < J_{k+1}$ as long as $X(t_{J_k}, J_k) \in K_X$, which is (78). //

**Claim 3:** Thanks to (71d), (71e), and (71f), the inequality (76) holds for all $(t, j) \in \text{dom } X$, and thus, inequalities (77) and (78) holds for all $k \ge 0$.

**Proof of Claim 3:** Let us look at the derivative of $V_x$ on the set $\Omega_X$ with $\epsilon \in [0, 1]$, while recalling (69):

$$ 
\begin{aligned}
\dot{V}_x = & \nabla_x V_x \cdot f_x(x, h(x), 0) + \nabla_x V_x \cdot \tilde{f}_x(x, z + h(x), \epsilon) \\
& \le -\alpha_x(V_x(x)) + L_V^x \lambda_f(|z|) + L_V^x \lambda_f(\epsilon).
\end{aligned}
\quad (83) $$

Then it is seen that, if $\epsilon$ satisfies (71d) and $|z| \le \nu^*$, then $\dot{V}_x < 0$ on the level set $\{x : V_x(x) = \mu + \delta_\mu\}$. In fact, with (70) and (71e), the inequality (77) implies that $|z(t_{J_N}, J_N)| \le C_g^{-1}(\nu^*/2) \le \nu^*$, and from there, with (71f), the inequality (78) implies that $|z(t, j)| \le \nu^*$ for $t_{J_N} \le t \le t_{J_{N+1}}$ and $J_N \le j < J_{N+1}$. This implies that $V_x(x(t, j)) \le \mu + \delta_\mu$ for $(0, 0) \le (t, j) \le (t_{J_{N+1}}, J_{N+1})$, and this argument repeats. //

Putting altogether, we choose $\epsilon^*$ in Theorem 3 such that all the inequalities (71a)–(71f) hold with $\epsilon^*$. Then, the following arguments lead us to the estimate in (7b).

For $(t, j) \in \text{dom } X$, it follows from (64) that $j \ge (\sigma_1/N_0)(t/\epsilon)$, and therefore, for $k$ such that $J_k \le j < J_{k+1}$, we have

$$ k \ge \frac{\sigma_1}{N_0 n} \frac{t}{\epsilon}. \quad (84) $$

Then, from (81), (77), (78), and Assumption 2b, we get

$$ 
|z(t,j)| \le C_g \circ \underline{\alpha}_y^{-1} \left( \gamma_m^{\frac{\sigma_1}{N_0 n} \frac{t}{\epsilon}} \bar{\alpha}_y(|z(0,0)|) + \frac{L_V^z C_n(\epsilon)}{1 - \gamma_m} \right) + C_a(\epsilon). \quad (85) $$

Consider the class $KL$ function $\beta_y(s, \tau)$ in (15) and the class $K$ function

$$ \chi_y(\epsilon) := C_g \circ \underline{\alpha}_y^{-1} \left( 2 \frac{L_V^z C_n(\epsilon)}{1 - \gamma_m} \right) + C_a(\epsilon). $$

Then (85) yields the inequality in (7b). To complete the proof of Theorem 2, inequality (7a) is derived from (85) using Lemma 3 in Section V-C.

## VII. CONCLUSIONS

Motivated by the applications of hybrid dynamical systems where frequent jumps are necessary for stability, we propose analysis of such systems based on singular perturbation viewpoint. The dynamics are decomposed into a continuous-time quasi-steady-state system and a discrete-time boundary-layer system. Under the assumption that the subsystems resulting from this decomposition satisfy appropriate stability assumptions, we can derive estimates on the norm of the trajectories showing practical convergence. In the process, we see the trade-offs between the ultimate bound and the parameter describing the frequency of jumps. An important application of our results is seen in the study of network control with heterogeneous agents where the graph describing the information exchange between agents is switching.

Several interesting research directions may emanate from the framework adopted in this work. Firstly, it may be possible to improve the estimates on system trajectories reported in Theorem 3 since the proof involves computation of certain constants over the compact sets for the general class of nonlinear systems. For some specific system class such as linear systems, or globally Lipschitz dynamics, one may modify the proof techniques for getting tighter estimates. Also, it may be possible to identify conditions on system data under which we can establish asymptotic stability for the limiting value of the perturbation parameter. Naturally, these questions open up the possibility of investigating more applications where the proposed setup of singularly perturbed hybrid systems is useful for analysis. This includes especially the control design problem related to stabilization over networks where frequent jumps in the state of the system contribute to stability of the system.

## REFERENCES

[1] A.A. Ahmadi, R.M. Jungers, P.A. Parrilo, and M. Roozbehani. Joint spectral radius and path-complete graph Lyapunov functions. *SIAM Journal on Control and Optimization*, 52(1):687–717, 2014.

[2] D. Angeli, E.D. Sontag, and Y. Wang. A characterization of integral input-to-state stability. *IEEE Transactions on Automatic Control*, 45(6):1082–1097, 2000.

[3] F. Bullo. *Lectures on network systems*. Kindle Direct Publishing Seattle, DC, USA, 2020.

[4] M. Cao, A.S. Morse, and B.D.O. Anderson. Reaching a consensus in a dynamically changing environment: A graphical approach. *SIAM Journal on Control and Optimization*, 47(2):575–600, 2008.

[5] M. Coraggio, P. De Lellis, and M. di Bernardo. Distributed discontinuous coupling for convergence in heterogeneous networks. *IEEE Control Systems Letters*, 5(3):1037–1042, 2020.

[6] R. Geiselhart, R.H. Gielen, M. Lazar, and F.R. Wirth. An alternative converse Lyapunov theorem for discrete-time systems. *Systems & Control Letters*, 70:49–59, 2014.

[7] A. Girard and P. Mason. Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems. *IEEE Control Systems Letters*, 3(3):499–504, 2019.

[8] R. Goebel. Discrete-time switching systems as difference inclusions: Deducing converse Lyapunov results for the former from those for the latter. *IEEE Transactions on Automatic Control*, 68(6):3694–3697, 2023.

[9] R. Goebel, R. Sanfelice, and A.R. Teel. *Hybrid Dynamical Systems: Modeling, Stability, and Robustness*. Princeton University Press, 2012.

[10] J.P. Hespanha, D. Liberzon, and A.R. Teel. Lyapunov conditions for input-to-state stability of impulsive systems. *Automatica*, 44:2735–2744, 2008.

[11] J.P. Hespanha and A.S. Morse. Stability of switched systems with average dwell-time. In *Proceedings of 38th IEEE Conference on Decision and Control*, pages 2655–2660, 1999.
[12] A. Jadbabaie, J. Lin, and A.S. Morse. Coordination of groups of mobile autonomous agents using nearest neighbor rules. *IEEE Transactions on Automatic Control*, 48(6):988–1001, 2003.

[13] C. Kellett and A.R. Teel. Smooth Lyapunov functions and robustness of stability for difference inclusions. *Systems & Control Letters*, 52:395–405, 2004.

[14] H.K. Khalil. *Nonlinear Systems*. Prentice-Hall, Inc., 3rd edition, 2002.

[15] J. Kim, J. Yang, H. Shim, J.-S. Kim, and J.H. Seo. Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents. *IEEE Transactions on Automatic Control*, 61(10):3096–3102, 2016.

[16] J.G. Lee and H. Shim. A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling. *Automatica*, 117, 2020.

[17] D. Liberzon, D. Nesic, and A.R. Teel. Lyapunov-based small-gain theorems for hybrid systems. *IEEE Transactions on Automatic Control*, 59(6):1395–1410, 2014.

[18] Y. Lin, E.D. Sontag, and Y. Wang. A smooth converse Lyapunov theorem for robust stability. *SIAM Journal on Control and Optimization*, 34(1):124–160, 1996.

[19] S. Liu, A. Tanwani, and D. Liberzon. ISS and integral-ISS of switched systems with nonlinear supply functions. *Mathematics of Controls, Signals, and Systems*, 34(2):297–327, 2022.

[20] M. Mesbahi and M. Egerstedt. *Graph Theoretic Methods in Multiagent Networks*. Princeton University Press, 2010.

[21] A.S. Morse. Supervisory control of families of linear set-point controllers Part I: Exact matching. *IEEE Transactions on Automatic Control*, 41(10):1413–1431, 1996.

[22] E. Panteley and A. Loria. Synchronization and dynamic consensus of heterogeneous networked systems. *IEEE Transactions on Automatic Control*, 62(8):3758–3773, 2017.

[23] P. Pepe. Converse Lyapunov theorems for discrete-time switching systems with given switches digraphs. *IEEE Transactions on Automatic Control*, 64(6):2502–2508, 2019.

[24] L. Perko. *Differential Equations and Dynamical Systems*. Springer, Newyork, third edition, 2000.

[25] M. Philippe, R. Essick, G. Dullerud, and R.M. Jungers. Stability of discrete-time switching systems with constrained switching sequences. *Automatica*, 72:242–250, 2016.

[26] R.G. Sanfelice and A.R. Teel. On singular perturbations due to fast actuators in hybrid control systems. *Automatica*, 47:692–701, 2011.

[27] R. Sepulchre, M. Janković, and P.V. Kokotović. *Constructive Nonlinear Control*. Communications and Control Engineering Series. Springer-Verlag, Berlin, 1997.

[28] A. Tanwani and H. Shim. Lyapunov functions for singularly perturbed hybrid systems with frequent jump dynamics. In *Proceedings of 60th IEEE Conference on Decision and Control*, pages 5382–5387, 2021.

[29] A.R. Teel, L. Moreau, and D. Nesic. A unified framework for input-to-state stability in systems with two time scales. *IEEE Transactions on Automatic Control*, 48(9):1526–1544, 2003.

[30] A.R. Teel and L. Praly. Tools for semiglobal stabilization by partial state and output feedback. *SIAM Journal on Control and Optimization*, 33(5):1443–1488, 1995.

[31] W. Wang, A.R. Teel, and D. Nesic. Analysis for a class of singularly perturbed hybrid systems via averaging. *Automatica*, 48:1057–1068, 2012.

[32] X.-F. Wang, K.-Z. Liu, X.-M. Sun, and A.R. Teel. Lyapunov-based singular perturbation results in the framework of hybrid systems. *IFAC-PapersOnLine*, 53(2):2027–2032, 2020.

which is well-defined. We immediately get

$$ \alpha_1(s) \le \begin{cases} \delta, & 0 \le s < s^*, \\ v_{12}\alpha_2(s), & s^* \le s \le \bar{s} \end{cases} $$

which is the desired inequality. □

<figure><img src="image_3.png" /></figure>

# APPENDIX A

## A LEMMA FOR CLASS K FUNCTIONS

**Lemma 4.** Consider two class K functions $\alpha_1$ and $\alpha_2$. Fix $\bar{s} > 0$. For every $\delta > 0$, there exists $v_{12} > 0$ such that

$$ \alpha_1(s) \le v_{12}\alpha_2(s) + \delta, \quad \forall s \in [0, \bar{s}]. $$

*Proof.* For a given $\delta > 0$, pick a positive $s^* \le \bar{s}$ such that $\alpha_1(s^*) \le \delta$, and take

$$ v_{12} = \max_{s \in [s^*, \bar{s}]} \frac{\alpha_1(s)}{\alpha_2(s)} $$