# A discrete-time averaging theorem and its application to zeroth-order Nash
equilibrium seeking

Suad Krilašević  Sergio Grammatico  Delft Center for Systems and Control, TU
Delft, The Netherlands

###### Abstract

In this paper we present an averaging technique applicable to the design of
zeroth-order Nash equilibrium seeking algorithms. First, we propose a multi-
timescale discrete-time averaging theorem that requires only that the
equilibrium is semi-globally practically stabilized by the averaged system,
while also allowing the averaged system to depend on “fast” states.
Furthermore, sequential application of the theorem is possible, which enables
its use for multi-layer algorithm design. Second, we apply the aforementioned
averaging theorem to prove semi-global practical convergence of the zeroth-
order information variant of the discrete-time projected pseudogradient
descent algorithm, in the context of strongly monotone, constrained Nash
equilibrium problems. Third, we use the averaging theory to prove the semi-
global practical convergence of the asynchronous pseudogradient descent
algorithm to solve strongly monotone unconstrained Nash equilibrium problems.
Lastly, we apply the proposed asynchronous algorithm to the connectivity
control problem in multi-agent systems.

###### keywords:

Averaging theorem, equilibrium seeking, asynchronous algorithm

††thanks: This work was partially supported by the ERC under research project
COSMOS (802348). E-mail addresses: {s.krilasevic-1<EMAIL_ADDRESS>
and

##  1 Introduction

Given a complex dynamical system, averaging techniques are used to construct a
simpler system, called the averaged system, that is easier to analyze than the
given one. Ideally, the averaged system should satisfy certain properties so
that it is possible to infer stability properties of the original system based
on the averaged one. These techniques are used extensively in extremum seeking
results, in continuous-time systems [11], [6], [17], discrete-time systems
[4], [23], and hybrid systems [20], [19].

_Literature review:_ Discrete-time averaging techniques have received intense
attention over the years. In [1], the authors show that the original dynamics
render the equilibrium exponentially stable under the assumption of
exponential stability of the equilibrium for the averaged dynamics.
Furthermore, they prove a similar result with similar assumptions for a mixed
time-scale system where the fast dynamics converge to zero. The requirement
for exponential stability is relaxed in [27] to just semi-global practical
asymptotic stability, for single time-scale systems. Furthermore, the authors
include noise into their analysis and provide input-to-state stability
results. In [29], the authors provide upper bounds for the time-scale
separation parameter in the case of linear switched averaged systems by using
a time-delay approach and similar assumptions as in [1]. The previously
mentioned single time-scale results assume that the jump mapping is time-
varying and that this dependence gets “smoothed out” using the averaging
technique. Thus, the main source of “perturbations” in the original system is
this time dependence. Likewise, it is possible to assume that jump mapping is
a function of some stochastic perturbations and that the goal of the averaging
is to “smooth out” the dependence on these perturbations. In [13], the authors
prove that under certain technical assumptions, the discrete-time stochastic
algorithm can be approximated by its continuous counterpart, and that
equilibrium of the original dynamics is weakly exponentially stable if the
equilibrium of the continuous counterpart is exponentially stable.
The usual approach to design of extremum seeking algorithms consists of
choosing a well-behaved full-information gradient-based algorithm in the case
of optimization, or pseudogradient-based in the case of games, integrated with
a (pseudo)gradient zeroth-order information estimation scheme [11], [6], [19],
[18]. The produced estimate then replaces the real value of the
(pseudo)gradient in the algorithm. A typical estimation technique it that of
injecting sinusoidal perturbations into the inputs of a cost function, whose
output is then correlated with the same perturbations. Via averaging
techniques, it can be proven that this estimation behaves as the
(pseudo)gradient, on average. The theory of averaging and singular
perturbations for continuous and hybrid systems [21], [28] enables the
adaptation of a wide spectrum of algorithms. In [11], the authors adapt the
gradient descent algorithm for the zeroth-order information case, together
with the additional high-pass and low-pass filters to improve performance. An
extremum seeking variant of the pseudogradient descent algorithm used for
solving unconstrained games is presented in [6]. Recently, the authors in [18]
propose a fixed-time zeroth-order algorithm for solving games, based on a
similar full-information fixed-time algorithm. An accelerated first-order
algorithm has been adapted for optimization problems in [19]. Unfortunately,
the same variety of extremum seeking algorithms in not available in discrete-
time due to the limitations of the discrete-time averaging theory. In [4], the
authors prove exponential convergence to the optimum of a quadratic function
under the zeroth-order variant of the gradient descent algorithm with
filtering. The authors in [30] prove ultimate boundness in a similar setup
where the plant is assumed to be general dynamic nonlinear and the
trajectories of the averaged system ultimately bounded. A similar approach is
used in [23] to prove convergence to the Nash equilibrium in a game without
constraints. In [13], the authors prove stability of its stochastic variant.
On the other hand, zeroth-order methods that use other approaches for gradient
estimation appear to be more successful and a recent overview for methods in
optimization can be found here [14]. The authors in [16] solve an N-coalition
game without local constraints for strongly monotone games by using Gaussian
smoothing to estimate the pseudogradient, while the authors in [24] propose an
algorithm for solving cooperative multi-agent cost minimization problem with
local constraints, also with Gaussian smoothing. Both papers assume
synchronous sampling of the agents, albeit with possible information delay.
Similar approach to Gaussian smoothing is the residual feedback estimator that
uses a previous evaluation of the cost function for the second point of the
pseudogradient approximation, thus reducing the numbers of cost functions
samples that need to be taken in one iteration. Using this approach, the
authors in [8] adapt two extra-gradient algorithms and prove convergence to
the Nash equilibrium in pseudo-monotone plus games for diminishing step sizes
and query radiuses. Authors in [25] and [26] estimate the pseudogradient using
the idea of continuous action-set learning automaton and prove convergence for
strictly monotone games and merely monotone games, respectively, via
diminishing step sizes, and Tikhonov regularization.
Asynchronous zeroth-order optimization algorithms have been well studied and
an overview can be found here [12]. For example, the authors in [22] use the
residual feedback estimator in an asynchronous gradient decent scheme to prove
convergence. In the current state of the art, zeroth-order discrete-time Nash
equilibrium seeking algorithms based on averaging use pseudogradient descent
without projections, while algorithms based on other methods are more general,
yet still assume synchronous sampling.
_Contribution_ : Motivated by the above literature and open research problems,
to the best of our knowledge, we consider an averaging technique for mixed
time-scale discrete-time systems and merely semi-globally practically
convergent averaged systems, with the application to the problem of learning
Nash equilibria via zeroth-order discrete-time algorithms, in the cases of
locally constrained agents, and asynchronous sampling. Specifically, our main
technical contributions are summarized next:

  * •

We extend the current results on averaging theory by using a mixed time-scale
formulation of the original system and requiring that the averaged systems
renders the equilibrium set SGPAS, unlike [27, Thm. 2], where a single time-
scale, time-variant system is used, and differently from [1, Thm. 2.2.4], [15,
Thm. 8.2.28] where exponential stability is needed and the fast subsystem
state is assumed to converge to the origin. Furthermore, we allow certain
types of additive perturbation dynamics to interfere with the nominal
averaging dynamics, and that the averaged jump mapping is a function of the
fast states, thus enabling easier _consecutive_ application of the averaging
theorem and the design of more complex algorithms.

  * •

Enabled by our extended averaging theory, we propose two novel zeroth-order
algorithms for game equilibrium seeking in discrete time. The first algorithm
solves the equilibrium in games with local constraints, differently from [23],
[16] where agents have no constraints; while the second one solves the problem
in the case where the agents are asynchronous, i.e. the agents do not sample
at the same time, nor do they coordinate in any way, differently from [23],
[16] where the agents sample synchronously.

_Notation_ : The set of real numbers and the set of nonnegative real numbers
are denoted by $\mathbb{R}$ and $\mathbb{R}_{+}$, respectively. Given a set
$\mathcal{Z}$, $\mathcal{Z}^{n}$ denotes the Cartesian product of $n$ sets
$\mathcal{Z}$. For a matrix $A\in\mathbb{R}^{n\times m}$, $A^{\top}$ denotes
its transpose. For vectors $x,y\in\mathbb{R}^{n}$ and $M\in\mathbb{R}^{n\times
n}$ a positive semi-definite matrix and $\mathcal{A}\subset\mathbb{R}^{n}$,
$\left\langle x\ \middle|\ y\right\rangle$, $\|x\|$, $\|x\|_{M}$ and
$\|x\|_{\mathcal{A}}$ denote the Euclidean inner product, norm, weighted norm
and distance to set respectively. Given $N$ vectors $x_{1},\dots,x_{N}$,
possibly of different dimensions, $\operatorname{col}\left({x_{1},\dots
x_{N}}\right)\coloneqq\left[x_{1}^{\top},\dots,x_{N}^{\top}\right]^{\top}$.
Collective vectors are denoted in bold, i.e,
$\boldsymbol{x}\coloneqq\operatorname{col}\left({x_{1},\dots,x_{N}}\right)$
and for each $i=1,\dots,N$,
$\boldsymbol{x}_{-i}\coloneqq\operatorname{col}\left({x_{1},\dots,x_{i-1},x_{i+1},\dots,x_{N}}\right)$
as they collect vectors from multiple agents. Given $N$ matrices $A_{1}$,
$A_{2}$, …, $A_{N}$, $\operatorname{blkdiag}\left({A_{1},\ldots,A_{N}}\right)$
denotes the block diagonal matrix with $A_{i}$ on its diagonal. Given a vector
$x$, $\operatorname{diag}(x)$ represents a diagonal matrix whose diagonal
elements are equal to the elements of the vector $x$. For a function
$v:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ differentiable in
the first argument, we denote the partial gradient vector as
$\nabla_{x}v(x,y)\coloneqq\operatorname{col}\left({\frac{\partial
v(x,y)}{\partial x_{1}},\ldots,\frac{\partial v(x,y)}{\partial
x_{N}}}\right)\in\mathbb{R}^{n}$. We use
$\mathbb{S}^{1}:=\left\\{z\in\mathbb{R}^{2}:z_{1}^{2}+z_{2}^{2}=1\right\\}$ to
denote the unit circle in $\mathbb{R}^{2}$. The set-valued mapping
$\text{N}_{S}:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ denotes the
normal cone operator for the set $S\subseteq\mathbb{R}^{n}$, i.e.,
$\text{N}_{S}(x)=\varnothing$ if $x\notin S$,
$\left\\{v\in\mathbb{R}^{n}|\sup_{z\in S}v^{\top}(z-x)\leq 0\right\\}$
otherwise. $\operatorname{Id}$ is the identity operator; $I_{n}$ is the
identity matrix of dimension $n$ and $\boldsymbol{0}_{n}$ is vector column of
$n$ zeros; their index is omitted where the dimensions can be deduced from
context. The unit ball of appropriate dimensions depending on context is
denoted with $\mathbb{B}$. A continuous function
$\gamma:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is of class $\mathcal{K}$ if
it is zero at zero and strictly increasing. A continuous function
$\alpha:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is of class $\mathcal{L}$ if
is non-increasing and converges to zero as its arguments grows unbounded. A
continuous function
$\beta:\mathbb{R}_{+}\times\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is of
class $\mathcal{KL}$ if it is of class $\mathcal{K}$ in the first argument and
of class $\mathcal{L}$ in the second argument. UGES and SGPAS refer to uniform
global exponential stability and semi-global practical asymptotic stability,
respectively, as defined in [19, Def. 2.2, 2.3].

##  2 Discrete-time averaging

We consider the following discrete-time system written in hybrid system
notation [7, Eq. 1.1, 1.2]

$\displaystyle\left\\{\begin{array}[]{ll}{u}^{+}&=u+\varepsilon G(u,\mu)\\\
{\mu}^{+}&=M(u,\mu)\end{array}\right.,\begin{array}[]{rl}(u,\mu)&\in\mathcal{U}\times\Omega.\end{array}$
(4)

where $u$ and $\mu$ are the state variables,
$\mathcal{U}\subset\mathbb{R}^{m}$, $\Omega\subset\mathbb{R}^{l}$,
$G:\mathcal{U}\times\Omega\rightarrow\mathbb{R}^{m}$ and
$M:\mathcal{U}\times\Omega\rightarrow\Omega$ are the state jump functions for
the states $u$ and $\mu$ respectively, and $\varepsilon>0$ is a small
parameter. Furthermore, the mapping $G$ is parametrized by a small parameter
$\gamma>0$, i.e. $G=G_{\gamma}$, but for notational convenience, this
dependence is not written explicitly in the equations.

Next we consider now the trajectories of the state $\mu$ oscillate
indefinitely and in turn create oscillations in the trajectory of the state
$u$. An equivalent system that produces trajectories $\bar{u}$ without
oscillations should be easier to analyze. We refer to such systems as averaged
systems [1, Eq. 2.2.12], [15, Eq. 8.33], and we focus on those of the
following form:

$\displaystyle\left\\{\begin{array}[]{ll}\tilde{u}^{+}=\tilde{u}+\varepsilon
G_{\textup{avg}}(\tilde{u},\tilde{\mu})\\\
{\tilde{\mu}}^{+}=M(\tilde{u},\tilde{\mu})\end{array}\right.,\begin{array}[]{rl}(\tilde{u},\tilde{\mu})&\in\mathcal{U}\times\Omega,\end{array}$
(8)

where $G_{\textup{avg}}:\mathcal{U}\times\Omega\rightarrow\mathbb{R}^{m}$ and
is also parametrized by $\gamma>0$. Unlike [1, Thm. 2.2.4], [15, Thm. 8.2.28],
we take into consideration the case where the function $G_{\textup{avg}}$
depends on the fast state $\tilde{\mu}$, not only on $\tilde{u}$.

To postulate the required relation between the function $G$ and the mapping
$G_{\textup{avg}}$, we should introduce an auxiliary system that describes the
behaviour of system (4) when the state $u$ is kept constant, i.e.
$\varepsilon=0$, the so-called boundary layer system [28, Eq. 6]:

$\displaystyle\left\\{\begin{array}[]{l}{{u}_{\textup{bl}}^{+}}=u_{\textup{bl}}\\\
{{\mu}_{\textup{bl}}^{+}}=M(u_{\textup{bl}},\mu_{\textup{bl}})\end{array}\right.,(u_{\textup{bl}},\mu_{\textup{bl}})\in\mathcal{U}\times\Omega.$
(11)

Thus, a function $G_{\textup{avg}}$ is called an average of the mapping $G$
with the boundary layer dynamics in (2) if the following condition holds true:

######  Assumption 1.

For any compact set $K\subset U$ and any solution
$(u_{\textup{bl}},\mu_{\textup{bl}})$ of (2) where $u_{\textup{bl}}$ is
contained in the compact set $K$, it holds that:

$\displaystyle\left\|\frac{1}{N}\sum_{i=0}^{N-1}\left[G\left(u_{\textup{bl}}(i),\,\mu_{\textup{bl}}(i)\right)-G_{\textup{avg}}\left(u_{\textup{bl}}(i),\mu_{\textup{bl}}(i)\right)\right]\right\|$
$\displaystyle\leq\sigma\left(N\right),$ (12)

for some function $\sigma:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ of class
$\mathcal{L}$. ∎

In plain words, in Assumption 1 we postulate that by using more samples over
time, our approximation of the mapping $G$ becomes better. Furthermore, let us
assume local Lipschitz continuity of the mappings as in [15, Assum. 8.2.13],
[27, Eq. 1, Def. 1] and compactness :

######  Assumption 2.

The functions $G$, $M$ and $G_{\textup{avg}}$ in (4), (8) are continuous in
their arguments and locally bounded; the mapping $G_{\textup{avg}}$ is locally
Lipschitz continuous in its first argument. The set $\Omega$ is compact. ∎

The averaging method can be used in unison with other algorithms via time-
scale separation. In such cases, often the averaged system does not
exponentially or asymptotically stabilize the equilibrium as in [1, Thm.
2.22], [15, Thm. 8.2.28], due to the introduction of perturbations from other
dynamics. Here, we assume the weaker property of semi-global practical
stability of the set $\mathcal{A}\times\Omega$ under $v$-perturbed dynamics of
the averaged system, where the perturbations are given by the dynamical system

$\displaystyle v^{+}=U(\tilde{u},\tilde{\mu},v)$ (13)

with $v\in\mathbb{R}^{m}$,
$U:\mathcal{U}\times\Omega\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ being
a function parametrized by a some $\varepsilon>0$.

######  Assumption 3.

Consider the system in (8) and perturbation dynamics in (13), respectively.
For any set $K\subset\mathcal{U}$, and all trajectories
$(\tilde{u},\tilde{\mu})$ contained in $K\times\Omega$, there exists a
function of class $\mathcal{K}$, $\overline{v}$, such that
$\max_{k\in\mathrm{dom}(v)}\left\|v(k)\right\|\leq\overline{v}(\varepsilon)$.
∎

######  Assumption 4.

The set $\mathcal{A}\times\Omega$ is SGPAS as $\gamma\rightarrow 0$ for the
dynamics in (8), perturbations in (13), and the corresponding Lyapunov
function $V_{\textup{a}}$ satisfies:

$\displaystyle\underline{\alpha}_{\textup{a}}\left(\left\|z\right\|_{\mathcal{A}}\right)\leq
V_{\textup{a}}(z,\mu)\leq\overline{\alpha}_{\textup{a}}\left(\left\|z\right\|_{\mathcal{A}}\right)$
(14a) $\displaystyle
V_{\textup{a}}(z^{+},\mu^{+})-V_{\textup{a}}(z,\mu)\leq-\tilde{\alpha}_{\varepsilon}\left({\varepsilon}\right)\alpha_{\textup{a}}\left(\left\|z\right\|_{\mathcal{A}}\right)$
$\displaystyle\text{ for
}\left\|z\right\|_{\mathcal{A}}\geq\alpha_{\gamma}(\gamma),$ (14b)

where $z=\tilde{u}+v$,
$\underline{\alpha}_{\textup{a}},\overline{\alpha}_{\textup{a}},\tilde{\alpha}_{\varepsilon},\alpha_{\textup{a}},\alpha_{\gamma}$
are functions of class $\mathcal{K}$, and the function
$\frac{\varepsilon}{\tilde{\alpha}_{\varepsilon}\left(\varepsilon\right)}$ is
bounded for $\varepsilon\in(0,\overline{\varepsilon})$. ∎

Under these assumptions, we claim that the original system is semi-global
practically asymptotically stable, as formalized next:

######  Theorem 1.

Let Assumptions 1, 2, 3 and 4 hold. The set $\mathcal{A}\times\Omega$ is SGPAS
as $(\varepsilon,\gamma)\rightarrow 0$ for the discrete dynamics in (4) with
perturbations in (13). The corresponding Lyapunov function $V_{\textup{a}}$
satisfies:

$\displaystyle\underline{\alpha}_{\textup{a}}\left(\left\|\xi\right\|_{\mathcal{A}}\right)\leq
V_{\textup{a}}(\xi,\mu)\leq\overline{\alpha}_{\textup{a}}\left(\left\|\xi\right\|_{\mathcal{A}}\right)$
$\displaystyle
V_{\textup{a}}(\xi^{+},\mu^{+})-V_{\textup{a}}(\xi,\mu)\leq-\hat{\alpha}_{\varepsilon}\left({\varepsilon}\right)\alpha_{\textup{a}}\left(\left\|\xi\right\|_{\mathcal{A}}\right)$
$\displaystyle\text{ for
}\left\|\xi\right\|_{\mathcal{A}}\geq\max\\{\alpha_{\gamma}(\gamma),{\alpha}_{\varepsilon}(\varepsilon)\\},$

where $\xi\coloneqq u+v+\eta$, $\eta$ is the perturbation state with dynamics

$\displaystyle\eta^{+}=(1-\varepsilon)\eta+\varepsilon[G_{\textup{avg}}(u,\mu)-G(u,\mu)],$
$\displaystyle\max_{k\in\mathbb{N}}\left\|\eta(k)\right\|\leq\overline{\eta}(\varepsilon),$
(15)

the $v$ dynamics are given by (13), and $\hat{\alpha}_{\varepsilon}$,
${\alpha}_{\varepsilon}$, $\overline{\eta}$ are functions of class
$\mathcal{K}$. ∎

See Appendix A. $\blacksquare$

##  3 Applications of the averaging theorem

In this section, we apply our averaging theorem, Theorem 1, to derive two
novel convergence results for NEPs. First, we propose a zeroth-order algorithm
for solving strongly monotone NEPs with local constraints _in discrete time_.
Secondly, we propose an algorithm for solving strongly monotone unconstrained
NEPs where the agents sample their states _asynchronously_.

###  3.1 Zeroth-order discrete time forward-backward algorithm

Let us consider a multi-agent system with $M$ agents indexed by
$i\in\mathcal{I}\coloneqq\\{1,2,\dots M\\}$, each with cost function

$\displaystyle J_{i}(x_{i},\boldsymbol{x}_{-i}),$ (16)

where $x_{i}\in\Omega_{i}\subset\mathbb{R}^{m_{i}}$ is the decision variable,
$J_{i}:\mathbb{R}^{m_{i}}\times\mathbb{R}^{m_{-i}}\rightarrow\mathbb{R}$,
$m\coloneqq\sum_{j\in\mathcal{I}}m_{j}$, $m_{-i}\coloneqq\sum_{j\neq i}m_{j}$,
$\Omega\coloneqq\Omega_{i}\times\dots\times\Omega_{N}$. Formally, let the goal
of each agent be to reach a steady state that minimizes their own cost
function, i.e.,

$\displaystyle\forall i\in\mathcal{I}:\ $
$\displaystyle\min_{x_{i}\in\Omega_{i}}J_{i}(x_{i},\boldsymbol{x}_{-i}).$ (17)

A popular solution to this problem is the so-called Nash equilibrium:

######  Definition 1 (Nash equilibrium).

A set of decision variables
$\boldsymbol{x}^{*}\coloneqq\operatorname{col}\left({x_{i}^{*}}\right)_{i\in\mathcal{I}}$
is a Nash equilibrium if, for all $i\in\mathcal{I}$,

$\displaystyle
x_{i}^{*}\in\underset{v_{i}\in\Omega_{i}}{\operatorname{argmin}}\
J_{i}\left(v_{i},\boldsymbol{x}_{-i}^{*}\right).\qed$ (18)

A fundamental mapping in NEPs is the pseudogradient mapping
$F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$, which is defined as:

$\displaystyle
F(\boldsymbol{x}):=\operatorname{col}\left(\left(\nabla_{x_{i}}J_{i}\left(x_{i},\boldsymbol{x}_{-i}\right)\right)_{i\in\mathcal{I}}\right).$
(19)

Let us also define
$C_{\textup{F}}\coloneqq\overline{\operatorname{co}}\\{F(\Omega)\\}$, the
convex hull of the image of the pseudogradient. To ensure existence and
uniqueness of the Nash equilibrium, we assume certain regularity properties
[2, Thm. 4.3]:

######  Assumption 5.

For each $i\in\mathcal{I}$, the function $J_{i}$ in (16) is continuously
differentiable in $x_{i}$ and continuous in $\boldsymbol{x}_{-i}$; the
function $J_{i}\left(\cdot,\boldsymbol{x}_{-i}\right)$ is strictly convex for
every fixed $\boldsymbol{x}_{-i}$. ∎

Furthermore, let us assume that no agent can compute their part of the the
pseudogradient $F$ directly, but they can only measure their instantaneous
cost $h_{i}=J_{i}(x_{i},\boldsymbol{x}_{-i})$, a common assumption in
extremum-seeking problems [10], [17], [20], [23]. The full-information problem
where $F$ is known can be solved in many ways, depending on the technical
assumptions on the problem data. Here we choose to study a simple forward-
backward algorithm [3, Equ. 26.14]:

$\displaystyle\boldsymbol{x}^{+}=(1-\lambda)\boldsymbol{x}+\lambda\mathrm{proj}_{C}\left(\boldsymbol{x}-\gamma
F(\boldsymbol{x})\right),$ (20)

for which the Lyapunov function
$V(\boldsymbol{x})=\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|^{2}$
satisfies the inequality

$\displaystyle
V(\boldsymbol{x}^{+})-V(\boldsymbol{x})\leq-\lambda(1-c)(2-\lambda
c)V(\boldsymbol{x}),$ (21)

where $c\coloneqq\frac{\sqrt{1+\gamma^{2}L^{2}}}{1+\gamma\mu_{F}}$ and
$\boldsymbol{x}^{*}$ is the Nash equilibrium of the game in (16). We note that
this Lyapunov function satisfies Assumption 4.
A naïve approach to adapting the algorithm in (20) for zeroth-order
implementation would be to use a gradient estimation scheme as in [4], [23]
and plug in the estimate directly into (20). However, because of the
projection, Assumption 1 would not be satisfied. Thus, an additional time-
scale separation is hereby proposed:

$\displaystyle\left\\{\begin{array}[]{ll}\boldsymbol{x}^{+}&=(1-\alpha\beta)\boldsymbol{x}+\alpha\beta\mathrm{proj}_{C}\left(\boldsymbol{x}-\gamma\boldsymbol{\xi}\right)\\\
\boldsymbol{\xi}^{+}&=(1-\alpha)\boldsymbol{\xi}+\alpha
2A^{-1}J(\boldsymbol{x}+A\mathbb{D}\boldsymbol{\mu})\mathbb{D}\boldsymbol{\mu}\\\
\boldsymbol{\mu}^{+}&=\mathcal{R}\boldsymbol{\mu}\end{array}\right.,$ (25)

where $\boldsymbol{\xi}\in\mathbb{R}^{m}$ are filter states,
$\boldsymbol{\mu}\in\mathbb{S}^{m}$ are the oscillator states,
$\alpha,\beta>0$ are small time-scale separation parameters,
$\mathcal{R}\coloneqq\operatorname{blkdiag}\left({(\mathcal{R}_{i})_{i\in\mathcal{I}}}\right)$,
$\mathcal{R}_{i}\coloneqq\operatorname{blkdiag}\left({\left[\begin{smallmatrix}\cos(\omega_{i}^{j})&-\sin(\omega_{i}^{j})\\\
\sin(\omega_{i}^{j})&\cos(\omega_{i}^{j})\end{smallmatrix}\right]_{j\leq
m_{i}}}\right)$, $\omega_{i}^{j}>0$ for all $i$ and $j$,
$\mathbb{D}\in\mathbb{R}^{m\times 2m}$ is a matrix that selects every odd row
from the vector of size $2m$, $a_{i}>0$ are small perturbation amplitude
parameters, $A\coloneqq\operatorname{diag}\left({(a_{i})_{i\leq m}}\right)$
and
$J(\boldsymbol{x})=\operatorname{blkdiag}\left({(J_{i}(x_{i},\boldsymbol{x}_{-i})I_{m_{i}})_{i\in\mathcal{I}}}\right)$.
We claim that the dynamics in (25) render the set
$\\{\boldsymbol{x}^{*}\\}\times C_{\textup{F}}\times\mathbb{S}^{m}$
practically stable. To the best of our knowledge, it is not possible to prove
convergence of the algorithm in (25), using the current averaging theory for
discrete-time systems, since [15, Thm. 8.2.28], [1, Thm. 2.22] require
exponential stability of the origin via the averaged system, and [27, Thm. 2]
does not incorporate boundary-layer dynamics. We claim that under the strong
monotonicity assumption of the pseudogradient, and a proper choice of the
perturbation frequencies, the algorithm in (25) converges to a Nash
equilibrium.

######  Assumption 6.

The pseudogradient mapping $F$ is $\mu_{\text{f}}$-strongly monotone and
$L$-Lipschitz continuous, i.e. $\left\langle\boldsymbol{x}-\boldsymbol{y}\
\middle|\
F(\boldsymbol{x})-F(\boldsymbol{y})\right\rangle\geq\mu_{\text{f}}\left\|\boldsymbol{x}-\boldsymbol{y}\right\|$,
$\left\|F(\boldsymbol{x})-F(\boldsymbol{y})\right\|\leq
L\left\|\boldsymbol{x}-\boldsymbol{y}\right\|$, for all
$(\boldsymbol{x},\boldsymbol{y})\in\mathbb{R}^{2m}$. ∎

######  Assumption 7.

The sets $(\Omega_{i})_{i\in\mathcal{I}}$ are convex, closed and bounded. ∎

######  Assumption 8.

The rotational frequencies of each agent $i$,
$\boldsymbol{\omega}_{i}=\operatorname{col}\left({(\omega_{i}^{j})_{j\leq
m_{i}}}\right)$, are chosen so that $\omega_{i}^{j}\pm\omega_{u}^{v}\neq 2\pi
z^{\prime},z^{\prime}\in\mathbb{Z}$, for every $u\in\mathcal{I}$, for every
$j\in\\{1,\dots,m_{i}\\}$, for every $v\in\\{1,\dots,m_{u}\\}$, apart for the
case when $i=u$ and $j=v$. ∎

######  Theorem 2.

Let Assumptions 5, 6, 7 and 8 hold. The set $\\{\boldsymbol{x}^{*}\\}\times
C_{\textup{F}}\times\mathbb{S}^{m}$ is SGPAS as
$(\alpha,\overline{a},\beta)\rightarrow 0$ for the dynamics in (25). ∎

See Appendix B. $\blacksquare$

###  3.2 Asynchronous zeroth-order discrete time forward algorithm

####  3.2.1 First-order information feedback

We now consider the same NEP as in Section 3.1 with
$\Omega_{i}\coloneqq\mathbb{R}^{m_{i}}$, but where the agent are
_asynchronous_ , i.e. each agent samples their states independently of others
without a global clock for synchronization. For ease of exposition, we assume
that the initial conditions are chosen so that simultaneous sampling never
occurs. In the full-information case, such algorithm can be represented in the
following form:

$\displaystyle\operatorname{col}\left({\dot{x}_{i},\dot{\tau}_{i},\dot{\kappa_{i}},\dot{t}}\right)=\operatorname{col}\left({\boldsymbol{0},\tfrac{1}{T_{i}},0,1}\right)$
$\displaystyle\text{ if
}\operatorname{col}\left({x_{i},\tau_{i},\kappa_{i},t}\right)\in\mathbb{R}^{m_{i}}\times[0,1]\times\mathbb{N}\times\mathbb{R}$
(26a)
$\displaystyle\left\\{\begin{array}[]{ll}x_{i}^{+}&=x_{i}-\,\alpha\nabla_{x_{i}}J_{i}(x_{i},\boldsymbol{x}_{-i})\\\
\tau_{i}^{+}&=0\\\ \kappa_{i}^{+}&=\kappa_{i}+1\\\ t^{+}&=t\end{array}\right.$
(26f) $\displaystyle\text{ if
}\operatorname{col}\left({x_{i},\tau_{i},\kappa_{i},t}\right)\in\mathbb{R}^{m_{i}}\times\\{1\\}\times\mathbb{N}\times\mathbb{R},$
(26g)

which in collective form reads as

$\displaystyle\operatorname{col}\left({\dot{\boldsymbol{x}},\dot{\boldsymbol{\tau}},\dot{\boldsymbol{\kappa}},\dot{k},\dot{t}}\right)=\operatorname{col}\left({\boldsymbol{0},\boldsymbol{T}^{-1},\boldsymbol{0},0,1}\right)$
$\displaystyle\text{ if
}\operatorname{col}\left({\boldsymbol{x},\boldsymbol{\tau},\boldsymbol{\kappa},k,t}\right)\in\mathbb{R}^{m}\times\mathcal{T}\times\mathbb{N}^{N+1}\times\mathbb{R}$
(27a)
$\displaystyle\left\\{\begin{array}[]{ll}\boldsymbol{x}^{+}&=\boldsymbol{x}-\alpha
S_{x}(\boldsymbol{\tau})F(\boldsymbol{x})\\\
\boldsymbol{\tau}^{+}&=(I-S_{\tau}(\boldsymbol{\tau}))\boldsymbol{\tau}\\\
\boldsymbol{\kappa}^{+}&=\boldsymbol{\kappa}+S_{\tau}(\boldsymbol{\tau})\\\
k^{+}&=k+1\\\ t^{+}&=t\end{array}\right.$ (27g) $\displaystyle\text{ if
}\operatorname{col}\left({\boldsymbol{x},\boldsymbol{\tau},\boldsymbol{\kappa},k,t}\right)\in\mathbb{R}^{m}\times\mathcal{T}_{\textup{R}}\times\mathbb{N}^{N+1}\times\mathbb{R},$
(27h)

where $\tau_{i}$ are timer states, $t$ is the “experienced” global time,
$\boldsymbol{T}^{-1}\coloneqq\operatorname{col}\left({(\tfrac{1}{T_{i}})_{i\in\mathcal{I}}}\right)$
is the vector of inverse sampling times, $\mathcal{T}\subset[0,1]^{N}$ is a
closed invariant set in which all of the timers evolve and it excludes the
initial conditions and their neighborhood for which we have concurrent
sampling,
$\mathcal{T}_{\textup{R}}\coloneqq\left(\cup_{i\in\mathcal{I}}[0,1]^{i-1}\times\\{1\\}\times[0,1]^{N-i}\right)\cap\mathcal{T}$
is the set of timer intervals where at least one agent has triggered its
sampling, $\kappa_{i}$ are private event counters, $k$ is the global event
counter, $S_{x}:\mathcal{T}\rightarrow\mathbb{R}^{m\times m}$ and
$S_{\tau}:\mathcal{T}\rightarrow\mathbb{R}^{N\times N}$ are continuous
functions that output diagonal matrices with ones on the positions that
correspond to states and timers of agents with $\tau_{i}=1$, respectively,
while other elements are equal to zero, when evaluating at
$\boldsymbol{\tau}\in\mathcal{T}_{\textup{R}}$. We note that the functions
$S_{x},S_{\tau}$ are introduced only to write down the algorithm in the
collective form, while the agents themselves do not require them for their
dynamics and in fact just follow (26). Furthermore, the counter states
$\kappa_{i},k$ and global time $t$ are not necessary for the algorithm
convergence, yet they help with understanding the setup of the algorithm. We
choose to represent the algorithm in the hybrid dynamical system framework to
replicate the behaviour of sampled systems with different sampling periods,
and to see its effects on the functions $S_{x},S_{\tau}$. Later, we represent
and model the system as a fully discrete-time system in order to study
convergence.
First, we show that the solution $\boldsymbol{\tau}(t,j)$ is periodic and that
implies that $S_{x}(\boldsymbol{\tau}(t,j))$ and
$S_{\tau}(\boldsymbol{\tau}(t,j))$ are also periodic. We make the following
assumption:

######  Assumption 9.

There exist natural numbers $\left(p_{i}\right)_{i\in\mathcal{I}}$, such that
the proportion $T_{1}:T_{2}:\dots:T_{N}=p_{1}:p_{2}\dots:p_{N}$ holds, where
$\left(T_{i}\right)_{i\in\mathcal{I}}$ are the sampling times. ∎

######  Lemma 1.

Let Assumption 9 hold. Denote $r_{i}=\tfrac{p}{p_{i}}$ and
$r=\sum_{i\in\mathcal{I}}r_{i}$, where $p$ is the least common multiple of
$(p_{i})_{i\in\mathcal{I}}$. For any trajectory
$S_{x}(\boldsymbol{\tau}(t,j))$ and $S_{\tau}(\boldsymbol{\tau}(t,j))$, where
$\boldsymbol{\tau}(t,j)$ is a solution of the system in (27), there exists
$T>0$ such that
$S_{x}(\boldsymbol{\tau}(t,j))=S_{x}(\boldsymbol{\tau}(t+T,j+r))$ and
$S_{\tau}(\boldsymbol{\tau}(t,j))=S_{\tau}(\boldsymbol{\tau}(t+T,j+r))$ for
all $(t,j)\in\mathrm{dom}(\boldsymbol{\tau})$ such that a jump occurred at
time $t$. ∎

See Appendix F. $\blacksquare$Because the values of $S_{x}$ and $S_{\tau}$ are
used only during jumps, we define

$\displaystyle\hat{S}_{x}(k;\boldsymbol{\tau}(0,0))=S_{x}(\max_{(t\in\mathrm{dom}(\boldsymbol{\tau}(\cdot,k))}t,k)$
(28)
$\displaystyle\hat{S}_{\tau}(k;\boldsymbol{\tau}(0,0))=S_{\tau}(\max_{(t\in\mathrm{dom}(\boldsymbol{\tau}(\cdot,k))}t,k),$
(29)

where functions $\hat{S}_{x}:\mathbb{N}\rightarrow\mathbb{R}^{m\times m}$ and
$\hat{S}_{\tau}:\mathbb{N}\rightarrow\mathbb{R}^{N\times N}$ are parametrized
by the vector of initial conditions of the timers, since different initial
conditions can change the order in which the agents are sampling their
actions. Due to Lemma, 1, for every initial condition
$\boldsymbol{\tau}(0,0)=\boldsymbol{\tau}^{0}$, it follows that
$\hat{S}_{x}(k,\boldsymbol{\tau}^{0})=\hat{S}_{x}(k+r,\boldsymbol{\tau}^{0})$
for all $k\in\mathbb{N}$.
Now we consider the following discrete time systems

$\displaystyle x_{i}(k+1)$
$\displaystyle=x_{i}(k)-\,\alpha\hat{S}^{i}_{x}(k;\boldsymbol{\tau}_{0})\nabla_{x_{i}}J_{i}(x_{i}(k),\boldsymbol{x}_{-i}(k))$
(30)

which in collective form read as

$\displaystyle\boldsymbol{x}(k+1)$
$\displaystyle=\boldsymbol{x}(k)-\alpha\hat{S}_{x}\left(k;\boldsymbol{\tau}_{0}\right)F(\boldsymbol{x}(k)),$
(31)

where the function $\hat{S}^{i}:\mathbb{N}\rightarrow\mathbb{R}^{m_{i}\times
m_{i}}$ returns the rows of $\hat{S}_{x}^{i}(k;\boldsymbol{\tau}_{0})$
corresponding to agent $i$. We can show that for every solution of (31) there
exists a corresponding solution of (27) and vice versa. We claim that under
the strong monotonicity assumption, an additional regularity assumption due to
the unboundedness of the decision set, and proper choice of the parameter
$\alpha$, the dynamics in (31) converge to the solution of the game, with the
same minimal convergence rate, regardless of the initial conditions of the
timers.

######  Assumption 10.

For each $i\in\mathcal{I}$, the function $J_{i}(\cdot,\boldsymbol{x}_{-i})$ in
(16) is radially unbounded for every fixed $\boldsymbol{x}_{-i}$. ∎

######  Theorem 3.

Let Assumptions 5, 6, 9 and 10 hold. Then, for all vectors of initial
conditions $\boldsymbol{\tau}_{0}$, there exists $\alpha^{*}$, such that for
any $\alpha\in(0,\alpha^{*})$, the NE solution $\boldsymbol{x}^{*}$ is UGES
for the dynamics in (31). Furthermore, the corresponding Lyapunov function
satisfies Assumption 4. ∎

See Appendix E. $\blacksquare$Moreover, for the hybrid system representation
in (27), since the trajectories of
$(\boldsymbol{\tau},\boldsymbol{\kappa},k,t)$ are invariant to the set
$\mathcal{T}\times\mathbb{N}^{N+1}\times\mathbb{R}$, and by the structure of
the flow and jump sets in (27) that assures complete solutions with unbounded
time and jump domains, it follows that the dynamics in (27) render the set
$\\{\boldsymbol{x}^{*}\\}\times\mathcal{T}\times\mathbb{N}^{N+1}\times\mathbb{R}$
UGES, as formalized next.

######  Corollary 1.

Let the Assumptions 5, 6, 9 and 10 hold. Then, the set
$\\{\boldsymbol{x}^{*}\\}\times\mathcal{T}\times\mathbb{N}^{N+1}\times\mathbb{R}$
is UGES for the dynamics in (27). Furthermore, the corresponding Lyapunov
function satisfies Assumption 4. ∎

####  3.2.2 Zeroth-order information feedback

Now, consider that each agent only has access to the measurements of the cost
function. They can modify the algorithm in (26) by implementing a
pseudogradient estimation scheme, similar to the one in Equation (25):

$\displaystyle\operatorname{col}\left({\dot{x}_{i},\dot{\xi}_{i},\dot{\mu}_{i},\dot{\tau}_{i},\dot{\kappa_{i}},\dot{t}}\right)=\operatorname{col}\left({\boldsymbol{0},\boldsymbol{0},\boldsymbol{0},\tfrac{1}{T_{i}},0,1}\right)$
(32a) $\displaystyle\text{ if
}\operatorname{col}\left({x_{i},\xi_{i},\mu_{i},\tau_{i},\kappa_{i},t}\right)\in\mathbb{R}^{2m_{i}}\times\mathbb{S}^{m}\times[0,1]\times\mathbb{N}\times\mathbb{R},$
$\displaystyle\left\\{\begin{array}[]{ll}x_{i}^{+}=x_{i}-\,\alpha\beta\xi_{i}\\\
\xi_{i}^{+}=(1-\alpha)\xi+\alpha\tfrac{2}{a_{i}}J_{i}(x+A\mathbb{D}\mu)\mathbb{D}_{i}\mu_{i}\\\
\mu_{i}^{+}=\mathcal{R}_{i}\mu_{i}\\\ \tau_{i}^{+}=0\\\
\kappa_{i}^{+}=\kappa_{i}+1\\\ t^{+}=t\end{array}\right.$ (32h)
$\displaystyle\text{ if
}\operatorname{col}\left({x_{i},\xi_{i},\mu_{i},\tau_{i},\kappa_{i},t}\right)\in\mathbb{R}^{2m_{i}}\times\mathbb{S}^{m}\times\\{1\\}\times\mathbb{N}\times\mathbb{R},$

which in the collective form reads as:

$\displaystyle\operatorname{col}(\dot{\boldsymbol{x}},\dot{\boldsymbol{\xi}},\dot{\boldsymbol{\mu}},\dot{\boldsymbol{\tau}},\dot{\boldsymbol{\kappa}},\dot{k},\dot{t})=\operatorname{col}\left({\boldsymbol{0},\boldsymbol{0},\boldsymbol{0},\boldsymbol{T}^{-1},\boldsymbol{0},0,1}\right)$
(33a) $\displaystyle\text{ if
}\operatorname{col}\left({\boldsymbol{x},\boldsymbol{\xi},\boldsymbol{\mu},\boldsymbol{\tau},\boldsymbol{\kappa},k,t}\right)\in\mathbb{R}^{2m}\times\mathbb{S}^{m}\times\mathcal{T}\times\mathbb{N}^{N+1}\times\mathbb{R},$
$\displaystyle\left\\{\begin{array}[]{ll}\boldsymbol{x}^{+}=\boldsymbol{x}-\alpha\beta
S_{x}(\boldsymbol{\tau})\boldsymbol{\xi}\\\
\boldsymbol{\xi}^{+}=\boldsymbol{\xi}+\alpha
S_{x}(\boldsymbol{\tau})\left(2A^{-1}J(\boldsymbol{x}+A\mathbb{D}\boldsymbol{\mu})\mathbb{D}\boldsymbol{\mu}-\boldsymbol{\xi}\right)\\\
\boldsymbol{\mu}^{+}=(I-S_{\mu}(\boldsymbol{\tau}))\boldsymbol{\mu}+S_{\mu}(\boldsymbol{\tau}))\mathcal{R}\boldsymbol{\mu}\\\
\boldsymbol{\tau}^{+}=(I-S_{\tau}(\boldsymbol{\tau}))\boldsymbol{\tau}\\\
\boldsymbol{\kappa}^{+}=\boldsymbol{\kappa}+S_{\tau}(\boldsymbol{\tau})\\\
k^{+}=k+1\\\ t^{+}=t\end{array}\right.$ (33i) $\displaystyle\text{ if
}\operatorname{col}\left({\boldsymbol{x},\boldsymbol{\xi},\boldsymbol{\mu},\boldsymbol{\tau},\boldsymbol{\kappa},k,t}\right)\in\mathbb{R}^{2m}\times\mathbb{S}^{m}\times\mathcal{T}_{\textup{R}}\times\mathbb{N}^{N+1}\times\mathbb{R},$

where $S_{\mu}:\mathcal{T}\rightarrow\mathbb{R}^{2m\times 2m}$ is a continuous
functions that outputs a diagonal matrix with ones on the positions that
correspond to oscillator states of agents with $\tau_{i}=1$, while other
elements are equal to zero, when evaluating at
$\boldsymbol{\tau}\in\mathcal{T}_{\textup{R}}$, and other notation is defined
as in (25) and (27).
Under the assumption of properly chosen perturbation frequencies, we claim
semi-global practical stability of the set of solutions.

######  Assumption 11.

The rotational frequencies of every agent $i$,
$\boldsymbol{\omega}_{i}=\operatorname{col}\left({(\omega_{i}^{j})_{j\leq
m_{i}}}\right)$, are chosen so that
$\omega_{i}^{j}r_{i}\pm\omega_{u}^{v}r_{j}\neq 2\pi
z^{\prime},z^{\prime}\in\mathbb{Z}$,
$r_{i}=\tfrac{p}{p_{i}},r_{j}=\tfrac{p}{p_{j}}$, for every $u\in\mathcal{I}$,
for every $j\in\\{1,\dots,m_{i}\\}$, for every $v\in\\{1,\dots,m_{u}\\}$,
apart for the case when $i=u$ and $j=v$. ∎

######  Theorem 4.

Let the Assumptions 5, 6, 9, 11 hold. The set
$\\{x^{*}\\}\times\mathbb{R}^{m}\times\mathbb{S}^{m}\times\mathcal{T}\times\mathbb{N}^{N+1}\times\mathbb{R}$
is SGPAS as $(\alpha,\overline{a},\beta)\rightarrow 0$ for the dynamics in
(33). ∎

The result is proven by following the same steps as the proof of Theorem 2 and
by using system in (27) with additional filtering state $\boldsymbol{\xi}$
like in (93) as the second averaged system. $\blacksquare$

##  4 Illustrative example

The connectivity control problem has been considered in [23] as a Nash
equilibrium problem. In many practical scenarios, multi-agent systems, besides
their primary objective, are designed to uphold certain connectivity as their
secondary objective. In what follows, we consider a similar problem in which
each agent is tasked with finding a source of an unknown signal while
maintaining certain connectivity. Unlike [23], we only consider the case
without vehicle dynamics.
Consider a system consisting of multiple agents indexed by
$i\in\mathcal{I}\coloneqq\\{1,\dots N\\}$. Each agent is tasked with locating
a source of a unique unknown signal. The strength of all signals abides by the
inverse-square law, i.e. proportional to $1/r^{2}$. Therefore, the inverse of
the signal strength can be used as a cost function. Additionally, the agents
must not drift apart from each other too much, as they should provide quick
assistance to each other in case of critical failure. This is enforced by
incorporating the signal strength of the fellows agents in the cost functions.
Thus, we design the cost

$\displaystyle\forall
i\in\mathcal{I}:J_{i}(\boldsymbol{x})=\|x_{i}-x_{i}^{s}\|^{2}+c\sum_{j\in\mathcal{I}_{-i}}\|x_{i}-x_{j}\|^{2},$
(34)

where $\mathcal{I}_{-i}\coloneqq\mathcal{I}\setminus\\{i\\}$, $c,b>0$ and
$x_{i}^{s}$ represents the position of the source assigned to agent $i$. Goal
of each agent is to minimize their cost function, and the solution to this
problem is a Nash equilibrium. Furthermore, agents are mutually independent so
their sampling time are not synchronised. To solve this problem, we use the
asynchronous pseudogradient descent algorithm in (33).
For our numerical simulations, we choose the parameters: $x^{s}_{1}=(-4,-8)$,
$x^{s}_{2}=(-12,-3)$, $x^{s}_{3}=(1,7)$, $x^{s}_{4}=(16,8)$, $c=0.04$,
$\gamma=0.1$, $\alpha=0.1$, $\beta=0.003$, $a_{i}=0$ for all $i$,
$T=(0.01,0.015,0.02,0.01)$, $\boldsymbol{\tau}(0,0)=(0,0.002,0.004,0.006)$,
the perturbation frequencies $\omega_{i}^{j}$ were chosen as different natural
numbers with added random numbers of maximal amplitude of 0.5.
The numerical results are illustrated on Figures 1 and 2. We note that the
trajectories converge to a small neighborhood of the Nash equilibrium. This
can be partially attributed to the robustness properties of the pseudogradient
descent with strongly monotone operators.

Figure 1: State trajectories in the $x_{1}-x_{2}$ plane. Circle symbols
represent locations of the sources, while the $\times$ symbols represent
locations of the NE. Perturbations signals are added to the states. Figure 2:
Time response of the states. The dashed lines correspond to the states of the
Nash equilibrium.

##  5 Conclusion

Averaging theory can be adapted for use in discrete systems with multiple
timescales. Furthermore, strongly monotone Nash equilibrium problem with
constrained action sets, or with asynchronous action sampling, can be solved
via zeroth-order discrete-time algorithms that leverage novel averaging theory
results.

## References

  * [1] E-W Bai, L-C Fu, and Sosale Shankara Sastry.  Averaging analysis for discrete time and sampled data adaptive systems.  IEEE Transactions on Circuits and Systems, 35(2):137–148, 1988\.
  * [2] Tamer Başar and Geert Jan Olsder.  Dynamic noncooperative game theory.  SIAM, 1998.
  * [3] Heinz H Bauschke, Patrick L Combettes, et al.  Convex analysis and monotone operator theory in Hilbert spaces, volume 408.  Springer, 2 edition, 2011.
  * [4] Joon-Young Choi, Miroslav Krstic, Kartik B Ariyur, and Jin Soo Lee.  Extremum seeking control for discrete-time systems.  IEEE Transactions on automatic control, 47(2):318–323, 2002.
  * [5] Yat Tin Chow, Tianyu Wu, and Wotao Yin.  Cyclic coordinate-update algorithms for fixed-point problems: Analysis and applications.  SIAM Journal on Scientific Computing, 39(4):A1280–A1300, 2017.
  * [6] Paul Frihauf, Miroslav Krstic, and Tamer Basar.  Nash equilibrium seeking in noncooperative games.  IEEE Transactions on Automatic Control, 57(5):1192–1207, 2011.
  * [7] Rafal Goebel, Ricardo G Sanfelice, and Andrew R Teel.  Hybrid dynamical systems.  Princeton University Press, 2012.
  * [8] Yuanhanqing Huang and Jianghai Hu.  Zeroth-order learning in continuous games via residual pseudogradient estimates, 2023.
  * [9] Hassan K Khalil.  Nonlinear systems.  Prentice Hall, 2002.
  * [10] Suad Krilašević and Sergio Grammatico.  Learning generalized Nash equilibria in multi-agent dynamical systems via extremum seeking control.  Automatica, 133:109846, 2021.
  * [11] Miroslav Krstić and Hsin-Hsiung Wang.  Stability of extremum seeking feedback for general nonlinear dynamic systems.  Automatica, 36(4):595–601, 2000.
  * [12] Xiangru Lian, Huan Zhang, Cho-Jui Hsieh, Yijun Huang, and Ji Liu.  A comprehensive linear speedup analysis for asynchronous stochastic parallel optimization from zeroth-order to first-order.  Advances in Neural Information Processing Systems, 29, 2016.
  * [13] Shu-Jun Liu and Miroslav Krstic.  Stochastic averaging in discrete time and its applications to extremum seeking.  IEEE Transactions on Automatic control, 61(1):90–102, 2015.
  * [14] Sijia Liu, Pin-Yu Chen, Bhavya Kailkhura, Gaoyuan Zhang, Alfred O Hero III, and Pramod K Varshney.  A primer on zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications.  IEEE Signal Processing Magazine, 37(5):43–54, 2020.
  * [15] Iven Mareels and Jan Willem Polderman.  Adaptive systems.  In Adaptive Systems, pages 1–26. Springer, 1996.
  * [16] Yipeng Pang and Guoqiang Hu.  Nash equilibrium seeking in n-coalition games via a gradient-free method.  Automatica, 136:110013, 2022.
  * [17] Jorge I Poveda and Miroslav Krstić.  Fixed-time gradient-based extremum seeking.  In 2020 IEEE American Control Conference (ACC), pages 2838–2843, 2020.
  * [18] Jorge I Poveda and Miroslav Krstić.  Nonsmooth extremum seeking control with user-prescribed fixed-time convergence.  IEEE Transactions on Automatic Control, 66(12):6156–6163, 2021\.
  * [19] Jorge I Poveda and Na Li.  Robust hybrid zero-order optimization algorithms with acceleration via averaging in time.  Automatica, 123:109361, 2021.
  * [20] Jorge I Poveda and Andrew R Teel.  A framework for a class of hybrid extremum seeking controllers with dynamic inclusions.  Automatica, 76:113–126, 2017.
  * [21] Ricardo G Sanfelice and Andrew R Teel.  On singular perturbations due to fast actuators in hybrid control systems.  Automatica, 47(4):692–701, 2011.
  * [22] Yi Shen, Yan Zhang, Scott Nivison, Zachary I Bell, and Michael M Zavlanos.  Asynchronous zeroth-order distributed optimization with residual feedback.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 3349–3354. IEEE, 2021.
  * [23] Miloš S Stankovic, Karl H Johansson, and Dušan M Stipanovic.  Distributed seeking of Nash equilibria with applications to mobile sensor networks.  IEEE Transactions on Automatic Control, 57(4):904–919, 2011.
  * [24] Yujie Tang, Zhaolin Ren, and Na Li.  Zeroth-order feedback optimization for cooperative multi-agent systems.  Automatica, 148:110741, 2023.
  * [25] Tatiana Tatarenko and Maryam Kamgarpour.  Learning generalized Nash equilibria in a class of convex games.  IEEE Transactions on Automatic Control, 64(4):1426–1439, 2018.
  * [26] Tatiana Tatarenko and Maryam Kamgarpour.  Bandit online learning of nash equilibria in monotone games.  arXiv preprint arXiv:2009.04258, 2020.
  * [27] Wei Wang and Dragan Nešić.  Input-to-state stability analysis via averaging for parameterized discrete-time systems.  In Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, pages 1399–1404. IEEE, 2009.
  * [28] Wei Wang, Andrew R Teel, and Dragan Nešić.  Analysis for a class of singularly perturbed hybrid systems via averaging.  Automatica, 48(6):1057–1068, 2012.
  * [29] Xuefei Yang, Jin Zhang, and Emilia Fridman.  Periodic averaging of discrete-time systems: A time-delay approach.  IEEE Transactions on Automatic Control, 2022.
  * [30] Hassan Zargarzadeh, Sarangapani Jagannathan, and James A Drallmeier.  Extremum-seeking for nonlinear discrete-time systems with application to hcci engines.  In 2014 American Control Conference, pages 861–866. IEEE, 2014\.

##  Appendix A Proof of Theorem 1

_Sketch of the proof:_ First, we show that under a change of coordinates, the
system in (4) can be represented as an inflated version of the averaged system
in (8). Then we show that the inflation can be arbitrarily small for small
enough $\varepsilon$. Finally, we use the stability properties of the averaged
system and the bounded inflation property to prove SGPAS.

By introducing an additional state, we construct the augmented system:

$\displaystyle\left\\{\begin{array}[]{l}{u}^{+}=u+\varepsilon G(u,\mu)\\\
{\mu}^{+}=M(u,\mu)\\\
\eta^{+}=(1-\varepsilon)\eta+\varepsilon\left[G_{\textup{avg}}(u,\mu)-G(u,\mu)\right]\\\
\end{array}\right.$ (38)
$\displaystyle(u,\mu,\eta)\in\mathcal{U}\times\Omega\times\mathbb{R}^{m}.$

With a change of coordinates $u=\tilde{u}-\eta$, $\mu=\tilde{\mu}$ the system
is transformed to:

$\displaystyle\left\\{\begin{array}[]{l}\tilde{u}^{+}=\tilde{u}+\varepsilon
G_{\textup{avg}}(\tilde{u}-\eta,\tilde{\mu})-\varepsilon\eta\\\
\tilde{\mu}^{+}=M(\tilde{u}-\eta,\tilde{\mu})\\\
\eta^{+}=(1-\varepsilon)\eta+\varepsilon\left[G_{\textup{avg}}(\tilde{u}-\eta,\tilde{\mu})-G(\tilde{u}-\eta,\tilde{\mu})\right]\end{array}\right.$
(42)
$\displaystyle(\tilde{u}-\eta,\tilde{\mu},\eta)\in\mathcal{U}\times\Omega\times\mathbb{R}^{m}$

We note that the $\tilde{u}$ dynamics in (42) are perturbed dynamics of the
averaged system in (8) and $\left\|\eta\right\|$ is the upper bound on the
perturbation amplitude. To prove our desired stability, we characterize the
bound of this amplitude:

######  Lemma 2.

For every $a>0$ and compact set $K\in\mathcal{U}$, there exists
$\varepsilon^{*}$ such that $\left\|\eta(t,j)\right\|<a$ holds for any
$\varepsilon\in(0,\varepsilon^{*}]$ and any trajectory of system in (42) where
$\tilde{u}$ is contained in the set $K$. ∎

For the purposes of the proof, we construct the concatenated trajectory
($u^{L},\,\mu^{L}$), which is created by taking solutions of length $L$ of the
system in (2) and concatenating them together.
We derive a similar bound to (12) for the concatenated trajectory type using
Assumption 1:

$\displaystyle\left\|\sum_{i=1}^{N}\left[G\left(u^{L}_{i},\mu_{i}^{L}\right)-G_{\textup{avg}}\left(u^{L}_{i},\mu_{i}^{L}\right)\right]\right\|$
$\displaystyle\leq\sum_{j=1}^{n}\left\|\sum_{i=1+L(j-1)}^{1+jL}\left[G\left(u_{i}^{L},\mu_{i}^{L}\right)-G_{\textup{avg}}\left(u^{L}_{i},\mu_{i}^{L}\right)\right]\right\|$
$\displaystyle+\left\|\sum_{i=nL+1}^{N}\left[G\left(u_{i}^{L},\mu_{i}^{L}\right)-G_{\textup{avg}}\left(u^{L}_{i},\mu_{i}^{L}\right)\right]\right\|$
$\displaystyle\leq nL\sigma(L)+(N-nL)\sigma(N-nL)$ $\displaystyle\leq
N\sigma(L)+L\sigma(0)$ (43)

Note that the bounds in (12) and (43) use (concatenated) boundary layer
trajectories instead of the $u$ trajectory in (4). In order to use the bound
in (43), we rewrite the $\eta$ dynamics in (38) as

$\displaystyle\eta^{+}=(1-\varepsilon)\eta+v_{1}+v_{2},$ (44)

where

$\displaystyle v_{1}$
$\displaystyle=\varepsilon\left[{G\left(u^{L},\mu^{L}\right)-G_{\textup{avg}}\left(u^{L}\right)}\right],$
(45) $\displaystyle v_{2}$
$\displaystyle=\varepsilon\left[{G\left(u,\mu\right)-G_{\textup{avg}}\left(u\right)-G\left(u^{L},\mu^{L}\right)+G_{\textup{avg}}\left(u^{L}\right)}\right].$
(46)

We use the superposition principle to determine the maximum value of $\eta$ by
analysing the inputs $v_{1}$ and $v_{2}$ separately. Let us start with
$v_{1}$. We append the subscripts to the notation of states to denote the time
index. The discrete dynamics are given by:

$\displaystyle\eta_{k+1}$
$\displaystyle=(1-\varepsilon)\eta_{k}+\varepsilon\left[G\left(u_{k}^{L},\mu_{k}^{L}\right)-G_{\textup{avg}}\left(u_{k}^{L}\right)\right].$
(47)

We define two additional variables:

$\displaystyle\phi_{k+1}$
$\displaystyle=\varepsilon\sum_{i=1}^{k}\left[G\left(u_{k}^{L},\mu_{k}^{L}\right)-G_{\textup{avg}}\left(u_{k}^{L},\mu_{k}^{L}\right)\right]$
(48) $\displaystyle\theta_{k}$ $\displaystyle=\eta_{k}-\phi_{k}$ (49)

From (38) and (50) it holds

$\displaystyle\theta_{k+1}$
$\displaystyle=\eta_{k+1}-\phi_{k+1}=(1-\varepsilon)\eta_{k}-\phi_{k}$
$\displaystyle=(1-\varepsilon)(\eta_{k}-\phi_{k})-(\varepsilon-1)\phi_{k}-\phi_{k}$
$\displaystyle=(1-\varepsilon)\theta_{k}-\varepsilon\phi_{k}$
$\displaystyle=(1-\varepsilon)\left[(1-\varepsilon)\theta_{k-1}-\varepsilon\phi_{k-1}\right]-\varepsilon\phi_{k}$
$\displaystyle=(1-\varepsilon)^{2}\theta_{k-1}-\varepsilon\left[(1-\varepsilon)\phi_{k-1}+\phi_{k}\right]$
$\displaystyle=\dots$
$\displaystyle=-\varepsilon\left[\sum_{i=0}^{k}(1-\varepsilon)^{i}\phi_{k-i}\right]$
$\displaystyle=\varepsilon\left[\sum_{i=0}^{k}(1-\varepsilon)^{i}\left(\phi_{k}-\phi_{k-i}\right)\right]-\varepsilon\left[\sum_{i=0}^{k}(1-\varepsilon)^{i}\phi_{k}\right]$
$\displaystyle=\varepsilon\left[\sum_{i=0}^{k}(1-\varepsilon)^{i}\left(\phi_{k}-\phi_{k-i}\right)\right]-\phi_{k}\left[1-(1-\varepsilon)^{k+1}\right]$
(50)

From (47), (48) and (50) we have

$\displaystyle\eta_{k+1}$ $\displaystyle=\theta_{k+1}+\phi_{k+1}$
$\displaystyle={\varepsilon\left[\sum_{i=0}^{k}(1-\varepsilon)^{i}\left(\phi_{k}-\phi_{k-i}\right)\right]}+\left(\phi_{k+1}-\phi_{k}\right)$
$\displaystyle+(1-\varepsilon)^{k+1}\phi_{k}.$ (51)

We use (43) in (51) to derive:

$\displaystyle\left\|\eta_{k+1}\right\|$
$\displaystyle\leq\underbrace{\varepsilon^{2}\left[\sum_{i=0}^{k}(1-\varepsilon)^{i}\left(i\,\sigma(L)+L\,\sigma(0)\right)\right]}_{S_{1}}+\varepsilon\sigma(L)$
$\displaystyle+\varepsilon
L\,\sigma(0)+\underbrace{\varepsilon(1-\varepsilon)^{k+1}\left(k\,\sigma(L)+L\,\sigma(0)\right)}_{S_{2}}$
(52)

To compute $S_{1}$, we start by find the sum

$\displaystyle\sum_{i=1}^{\infty}i(1-\varepsilon)^{i}$
$\displaystyle=\sum_{k=1}^{\infty}\sum_{i=k}^{\infty}(1-\varepsilon)^{i}=\sum_{k=1}^{\infty}\frac{(1-\varepsilon)^{k}}{\varepsilon}=\frac{1-\varepsilon}{\varepsilon^{2}}.$
(53)

Thus, we bound $S_{1}$ as follows:

$\displaystyle S_{1}\leq(1-\varepsilon)\sigma(L)+\varepsilon
L\sigma(0)\leq\sigma(L)+\varepsilon L\sigma(0).$ (54)

For $S_{2}$, we define the function $z(x)=x\,(1-\varepsilon)^{x}$. It is an
easy exercises to check that the maximum of the function is given by
$z\left(\frac{-1}{\log(1-\varepsilon)}\right)=\frac{-e}{\log(1-\varepsilon)}$.
Therefore, for the bound of $S_{2}$ we have

$\displaystyle S_{2}$
$\displaystyle\leq\varepsilon(1-\varepsilon)\left[\frac{-e}{\log\left(1-\varepsilon\right)}\sigma(L)+L\sigma(0)\right].$
(55)

As $\lim_{\varepsilon\rightarrow
0^{+}}\tfrac{\varepsilon}{\log\left(1-\varepsilon\right)}=1$, for small enough
$\varepsilon$, it follows:

$\displaystyle\leq
e(1-\varepsilon)\sigma(L)+\varepsilon(1-\varepsilon)L\,\sigma(0)$
$\displaystyle\leq e\sigma(L)+\varepsilon\,L\sigma(0).$ (56)

Finally, we have

$\displaystyle\left\|\eta_{k+1}\right\|$
$\displaystyle\leq(1+\varepsilon+e)\sigma(L)+3\varepsilon\,L\sigma(0),$ (57)

which holds for all $k$. The norm can be made arbitrarily small by the right
choice of parameters $L$ and $\varepsilon$.

Now, we move on to the input $v_{2}$. We define the inflated boundary layer
system:

$\displaystyle\left.\begin{array}[]{l}{{u}_{\textup{bl}}^{\delta+}}\in
u_{\textup{bl}}^{\delta}+\delta\mathbb{B}\\\
{{\mu}_{\textup{bl}}^{\delta+}}=M(u_{\textup{bl}}^{\delta},\mu_{\textup{bl}}^{\delta})\end{array}\right\\},(u_{\textup{bl}}^{\delta},\mu_{\textup{bl}}^{\delta})\in\mathcal{U}\times\Omega.$
(60)

We claim the following:

######  Lemma 3.

For any period $L$, positive real number $a$ and compact set
$K\in\mathcal{U}$, there exist a $\delta^{*}$ such that for every
$\delta\in(0,\delta^{*}]$ and for any trajectory
$(u_{\textup{bl}}^{\delta},\mu_{\textup{bl}}^{\delta})$ of the system in (60)
that is contained in $K\times\Omega$, there exist a concatenated trajectory
$(u^{L},\mu^{L})$ such that

$\displaystyle\|G\left(u_{\textup{bl}}^{\delta}(k),\mu_{\textup{bl}}^{\delta}(k)\right)-G_{\textup{avg}}\left(u_{\textup{bl}}^{\delta}(k)\right)-$
$\displaystyle
G\left(u^{L}(k),\mu^{L}(k)\right)+G_{\textup{avg}}\left(u^{L}(k)\right)\|\leq
a\qed$ (61)

$L,a,K$ are given. Let $\xi\coloneqq\operatorname{col}\left({u,\mu}\right)$.
Based on the continuity property of functions $G,G_{\textup{avg}}$, there
exists $\rho>0$ such that

$\displaystyle\left\|\xi_{1}-\xi_{2}\right\|\leq\rho\Rightarrow\left\|G(\xi_{1})-G(\xi_{2})\right\|\leq\tfrac{a}{2},$
$\displaystyle\left\|G_{\textup{avg}}(\xi_{1})-G_{\textup{avg}}(\xi_{2})\right\|\leq\tfrac{a}{2}.$
(62)

Next, we use [28, Lemma 2] for closeness of solutions of the inflated systems
with parameters $(0,L+1,\rho)$ and set $K\times\Omega$ to determine
$\delta^{*}$. That means that for every trajectory of the system in (60) where
$\xi_{\textup{bl}}^{\delta}(k)\in K\times\Omega$ for all
$k\in\mathrm{dom}(\xi_{\textup{bl}}^{\delta})$, there exists a trajectory
$\xi_{\textup{bl}}$ of the boundary layer system in (2), such that for each
$k\in\mathrm{dom}(\xi_{\textup{bl}}^{\delta})$ with $k\leq L+1$, we have
$\left\|\xi_{\textup{bl}}^{\delta}(k)-\xi_{\textup{bl}}(k)\right\|\leq\rho$.
As the inflated boundary layer system is time invariant, any sample shifted
trajectory is also a trajectory of the original system. Thus, for trajectories
starting in
$\xi_{\textup{bl}}^{\delta}(0),\xi_{\textup{bl}}^{\delta}(L),\dots,\xi_{\textup{bl}}^{\delta}(nL)$
with $n\in\mathbb{N},nL\in\mathrm{dom}(\xi_{\textup{bl}}^{\delta})$ there
exist trajectories (not necessary the same one) $\xi_{\textup{bl}}$ such that
the previous inequality holds for each segment of length $L$. We concatenate
these boundary layer trajectories into $\xi^{L}$ and write

$\displaystyle\left\|\xi_{\textup{bl}}^{\delta}(k)-\xi^{L}(k)\right\|\leq\rho\text{,
for }k\in\mathrm{dom}(\xi_{\textup{bl}}^{\delta}).$ (63)

From (62) and (63) we conclude (61). The trajectories $\tilde{u}(k),{\mu}(k)$
of the transformed system in (42), where $\tilde{u}(k)\in K$ for all
$k\in\mathrm{dom}\left(\tilde{u}\right)$ and $\left\|\eta(0)\right\|\leq a$,
are also trajectories of the inflated boundary layer system in (60) with

$\displaystyle\delta=\varepsilon\max_{u\in K,\mu\in\Omega,\eta\in
K_{\eta}}\\{G_{\textup{avg}}\left(u-\eta,\mu\right)-\eta\\},$ (64)

where $K_{\eta}$ is the set in which $\eta$ is contained during the trajectory
of the system. Let us prove that $K_{\eta}\subset a\mathbb{B}$ by first
showing that i.e. $\left\|\eta(1)\right\|\leq a$. First we find
$\varepsilon_{1}$ and $L$ such that
$(1+\varepsilon_{1}+e)\sigma(L)+3\varepsilon_{1}\,L\sigma(0)\leq\tfrac{a}{2}$.
Then we use the same $L$, positive number $\tfrac{a}{2}$ and set $K$ with
Lemma 3 to find $\delta^{*}$. For
$\varepsilon_{2}=\frac{\delta^{*}}{\max_{u\in K,\eta\in
a\mathbb{B}}\\{G_{\textup{avg}}\left(u-\eta\right)-\eta\\}}$, we guarantee
that for one step, the solution of (42) is also a solution of the inflated
boundary layer system in (60). Thus, for
$\varepsilon^{*}=\min\\{\varepsilon_{1},\varepsilon_{2}\\}$, we have that
variables in (45) and (46) are bounded as
$\left\|v_{1}\right\|\leq\tfrac{a}{2}$ and
$\left\|v_{2}\right\|\leq\tfrac{a}{2}$, and it follows from (44) that

$\displaystyle\left\|\eta(1)\right\|\leq\left\|(1-\varepsilon)\eta(0)\right\|+\varepsilon
a\leq a.$

The next sample will also be a solution of the $\delta$-inflated boundary
layer system and all of the previous bounds hold. Hence, the procedure can be
repeated with the same $\delta^{*}$ for all
$k\in\mathrm{dom}\left(\overline{u}\right)$, and it holds
$\left\|\varepsilon\eta(k)\right\|\leq a$. $\blacksquare$Now, we return to the
proof of Theorem 1. Let the set of initial conditions $K$ be given. From the
stability of the set $\mathcal{A}\times\Omega$ in Assumption 4 and the
dynamics in (42), we have:

$\displaystyle
V_{\textup{a}}(u^{+}+\eta^{+}+v^{+},\mu^{+})-V_{\textup{a}}(u+\eta+v,\mu)$
$\displaystyle=V_{\textup{a}}(\tilde{u}^{+}+v^{+},\tilde{\mu})-V_{\textup{a}}(\tilde{u}+v,\tilde{\mu})$
$\displaystyle\leq V_{\textup{a}}(\tilde{u}+\varepsilon
G_{\textup{avg}}(\tilde{u},\tilde{\mu})+v^{+},\tilde{\mu}^{+})-V_{\textup{a}}(\tilde{u}+v,\mu)$
$\displaystyle-V_{\textup{a}}(\tilde{u}+\varepsilon
G_{\textup{avg}}(\tilde{u},\tilde{\mu})+v^{+},\tilde{\mu}^{+})$
$\displaystyle+V_{\textup{a}}(\tilde{u}+\varepsilon
G_{\textup{avg}}(\tilde{u}-\eta,\tilde{\mu})-\varepsilon\eta+v^{+})$
$\displaystyle\leq-\tilde{\alpha}_{\varepsilon}\left({\varepsilon}\right)\alpha_{\textup{a}}\left(\left\|\tilde{u}+v\right\|_{\mathcal{A}}\right)+\varepsilon
L_{V_{\textup{a}}}(1+L_{\text{G}})\left\|\eta\right\|,$
$\displaystyle\quad\text{ for
}\left\|{\tilde{u}+v}\right\|_{\mathcal{A}}\geq\alpha_{\gamma}\left(\gamma\right)$
$\displaystyle\leq-\hat{\alpha}_{\varepsilon}\left({\varepsilon}\right)\alpha_{\textup{a}}\left(\left\|\tilde{u}+v\right\|_{\mathcal{A}}\right)$
$\displaystyle\quad\text{ for
}\left\|{\tilde{u}+v}\right\|_{\mathcal{A}}\geq\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\},$
(65)

where $L_{\text{G}}$ and $L_{V_{\textup{a}}}$ are Lipschitz constants of the
mapping $G_{\text{avg}}$ and function $V_{\text{a}}$ respectively,
${\alpha}_{\varepsilon}\left(\varepsilon\right)\geq\alpha_{\textup{a}}^{-1}\left[\frac{\varepsilon}{k\tilde{\alpha}_{\varepsilon}\left({\varepsilon}\right)}L_{V_{\textup{a}}}(1+L_{\text{G}})\left\|\eta\right\|\right]$,
$\hat{\alpha}_{\varepsilon}\left({\varepsilon}\right)\coloneqq(1-k)\tilde{\alpha}_{\varepsilon}\left({\varepsilon}\right)$
and $k\in(0,1)$. The function $\alpha_{\varepsilon}$ is a function of class
$\mathcal{K}$ on interval $(0,\bar{\varepsilon})$, as due to Assumption 4,
$\tfrac{\varepsilon}{\tilde{\alpha}_{\varepsilon}(\varepsilon)}$ is bounded on
that interval and $\eta$ can become arbitrarily small for proper choice of
$\varepsilon$, per Lemma 2. Finally, we plug in the states of the original
system to get

$\displaystyle V_{\textup{a}}(u^{+}+\eta^{+}+v^{+})-V_{\textup{a}}(u+\eta+v)$
$\displaystyle\leq-\hat{\alpha}_{\varepsilon}\left({\varepsilon}\right)\alpha_{\textup{a}}\left(\left\|u+\eta+v\right\|_{\mathcal{A}}\right)$
$\displaystyle\quad\text{ for
}\left\|{u+\eta+v}\right\|_{\mathcal{A}}\geq\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}.$
(66)

Let $\xi=u+\eta+v$. From the previous equation it follows

$\displaystyle V_{\textup{a}}(\xi(k))\leq
V_{\textup{a}}(\xi(0))-\sum_{i=0}^{k-1}\hat{\alpha}_{\varepsilon}\left({\varepsilon}\right)\alpha_{\textup{a}}\left(\left\|\xi(i)\right\|_{\mathcal{A}}\right)$
$\displaystyle\quad\text{ for
}\left\|{\xi(k)}\right\|_{\mathcal{A}}\geq\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}.$
(67)

Now, we move onto proving semi-global practical stability, Let $\Delta>\delta$
be any strictly positive real numbers. We choose parameters $\varepsilon$ and
$\gamma$ such that
$\overline{\eta}(\varepsilon)+\overline{v}(\gamma)\leq\tfrac{\delta}{4}$ and
$\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}\leq\tfrac{\delta}{4}$.
The conditional inequality in (67) is satisfied when
$\left\|{u(k)}\right\|_{\mathcal{A}}\geq\tfrac{\delta}{2}$.
_Semi-global stability_
For ease of notation, we drop the explicit dependence on $\varepsilon$ and
$\gamma$ in $\overline{\eta}(\varepsilon)$ and $\overline{v}(\gamma)$. We have
to show that for any $R>\delta$, there exists $r>0$, so that
$\left\|u(0)\right\|_{\mathcal{A}}\leq r$ implies that
$\left\|u(k)\right\|_{\mathcal{A}}\leq R$ for all
$k\in\mathrm{dom}\left(u\right)$. From (14a) and (67), it follows that

$\displaystyle\underline{\alpha}(\left\|\xi(k)\right\|_{\mathcal{A}})\leq
V_{\textup{a}}(\xi(k))\leq
V_{\textup{a}}(\xi(0))\leq\overline{\alpha}(\left\|\xi(0)\right\|_{\mathcal{A}})$
$\displaystyle\left\|\xi(k)\right\|_{\mathcal{A}}\leq\underline{\alpha}^{-1}\left(\overline{\alpha}(\left\|\xi(0)\right\|_{\mathcal{A}})\right)$
$\displaystyle\left\|u(k)\right\|_{\mathcal{A}}-\overline{\eta}-\overline{v}\leq\underline{\alpha}^{-1}\left(\overline{\alpha}(\left\|u(0)\right\|_{\mathcal{A}}+\overline{\eta}+\overline{v})\right)$
$\displaystyle\left\|u(k)\right\|_{\mathcal{A}}\leq\underline{\alpha}^{-1}\left(\overline{\alpha}(\left\|u(0)\right\|_{\mathcal{A}}+\overline{\eta}+\overline{v})\right)+\overline{\eta}+\overline{v}$
$\displaystyle\quad\text{ for
}\left\|{u}\right\|_{\mathcal{A}}\geq\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}+\overline{\eta}+\overline{v}$

From last equation it follows that
$R=\underline{\alpha}^{-1}\left(\overline{\alpha}(r+\overline{\eta}+\overline{v})\right)+\overline{\eta}+\overline{v}$.
Thus, it holds
$r=\overline{\alpha}^{-1}\left(\underline{\alpha}(R-\overline{\eta}-\overline{v})\right)-\overline{\eta}-\overline{v}$.
Considering that the infimum value of $R$ is $\delta$ and that
$\overline{\eta}+\overline{v}\leq\tfrac{\delta}{4}$, to ensure $r$ is
positive, we assume $\varepsilon$ and $\gamma$ are chosen so that
$\overline{\eta}+\overline{v}\leq\tfrac{1}{2}\overline{\alpha}^{-1}\left(\underline{\alpha}\left(\tfrac{3}{4}\delta\right)\right)$
which implies that
$r\leq\tfrac{1}{2}\overline{\alpha}^{-1}\left(\underline{\alpha}\left(\tfrac{3}{4}\delta\right)\right)$.
Furthermore, do assure that the Lyapunov difference is defined for those
radiuses, we impose an additional inequality on the tuning parameters:
$\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}+\overline{\eta}+\overline{v}<\tfrac{1}{2}\overline{\alpha}^{-1}\left(\underline{\alpha}\left(\tfrac{3}{4}\delta\right)\right)$.

_Practical attractivity_
We have to show that for any $R,r$ that satisfy $\Delta>R>r>\delta>0$, there
exists $T$, such that $\left\|u(0)\right\|_{\mathcal{A}}\leq R$ implies that
$\left\|u(k)\right\|_{\mathcal{A}}\leq r$ for all
$k\in\mathrm{dom}\left(u\right)$ and $k\geq T$. First, we use the bound we
derived in the proof of stability to define
$r^{\prime}\coloneqq\overline{\alpha}^{-1}\left(\underline{\alpha}(r-\overline{\eta}-\overline{v})\right)-\overline{\eta}-\overline{v}$,
from which we can conclude that $\left\|u(0)\right\|_{\mathcal{A}}\leq
r^{\prime}$ implies that $\left\|u(k)\right\|_{\mathcal{A}}\leq r$ for all
$k\in\mathrm{dom}\left(u\right)$. Let us define

$\displaystyle
T\coloneqq\left\lceil\frac{\overline{\alpha}\left(R+\overline{\eta}+\overline{v}\right)-\underline{\alpha}\left(r^{\prime}-\overline{\eta}-\overline{v}\right)}{\hat{\alpha}_{\varepsilon}\left(\varepsilon\right)\alpha\left(r-\overline{\eta}-\overline{v}\right)}\right\rceil+1.$
(68)

To prove via contradiction, we assume that $\left\|u(k)\right\|>r^{\prime}$
for all $k\leq T$. Now, by using the upper and lower bound of the Lyapunov
function on Equation in (67), it follows

$\displaystyle\underline{\alpha}\left(\left\|u(k)\right\|_{\mathcal{A}}-\overline{\eta}-\overline{v}\right)\leq\overline{\alpha}\left(R+\overline{\eta}+\overline{v}\right)$
$\displaystyle\quad\quad-k\hat{\alpha}_{\varepsilon}\left(\varepsilon\right)\alpha(r^{\prime}-\overline{\eta}-\overline{v})$
$\displaystyle\left\|u(k)\right\|_{\mathcal{A}}\leq\underline{\alpha}^{-1}\left(\overline{\alpha}\left(R+\overline{\eta}+\overline{v}\right)-k\hat{\alpha}_{\varepsilon}\left(\varepsilon\right)\alpha(r^{\prime}-\overline{\eta}-\overline{v})\right)$
$\displaystyle\quad\quad+\overline{\eta}+\overline{v}$ (69)

Let us choose $k=T-1$. When we plug in the chosen value of $k$ into inequality
(69), it follows that:

$\displaystyle\left\|u(k)\right\|_{\mathcal{A}}\leq r^{\prime},$ (70)

which leads us to a contradiction. Thus, in the first $T$ steps, $u(k)$
trajectory will enter at least once the set $A+r^{\prime}\mathbb{B}$. From the
stability properties, we know that once the trajectory enters aforementioned
set, it will never leave the set $\mathcal{A}+r\mathbb{B}$, which proves
practical attractivity.

Hence, to have semi-global practical asymptotic stability we have to choose
our parameters $\varepsilon,\gamma$ so that they satisfy inequalities

$\displaystyle\overline{\eta}(\varepsilon)+\overline{v}(\gamma)\leq\tfrac{\delta}{4}$
(71)
$\displaystyle\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}\leq\tfrac{\delta}{4}$
(72)
$\displaystyle\max\\{\alpha_{\gamma}\left(\gamma\right),{\alpha}_{\varepsilon}\left(\varepsilon\right)\\}+\overline{\eta}+\overline{v}<\tfrac{1}{2}\overline{\alpha}^{-1}\left(\underline{\alpha}\left(\tfrac{3}{4}\delta\right)\right)$
(73)

That concludes the proof of and semi-global practical asymptotic stability.
$\blacksquare$

##  Appendix B Proof of Theorem 2

First we show how to derive the boundary layer and averaged systems. Then we
show that we can apply Theorem 1 to prove stability.
The parameter $\alpha$ can be used as a time-scale separation parameter of the
first layer in algorithm (25). We derive the first boundary layer system
($\alpha=0$):

$\displaystyle\left\\{\begin{array}[]{ll}{\boldsymbol{x}_{\textup{bl}}^{1}}^{+}&=\boldsymbol{x}_{\textup{bl}}^{1}\\\
{\boldsymbol{\xi}_{\textup{bl}}^{1}}^{+}&=\boldsymbol{\xi}_{\textup{bl}}^{1}\\\
{\boldsymbol{\mu}_{\textup{bl}}^{1}}^{+}&=\mathcal{R}\boldsymbol{\mu}_{\textup{bl}}^{1}\end{array}\right.,$
(77)

and the first averaged system

$\displaystyle\left\\{\begin{array}[]{ll}\hat{\boldsymbol{x}}^{+}&=(1-\alpha\beta)\hat{\boldsymbol{x}}+\alpha\beta\mathrm{proj}_{C}\left(\hat{\boldsymbol{x}}-\gamma\boldsymbol{\xi}\right)\\\
\hat{\boldsymbol{\xi}}^{+}&=(1-\alpha)\hat{\boldsymbol{\xi}}+\alpha(F(\hat{\boldsymbol{x}})+\mathcal{O}(\overline{a}))\\\
\hat{\boldsymbol{\mu}}^{+}&=\mathcal{R}\hat{\boldsymbol{\mu}}\end{array}\right.,$
(81)

which is an $\mathcal{O}(\overline{a})$ inflation of the nominal averaged
system

$\displaystyle\left\\{\begin{array}[]{ll}\hat{\boldsymbol{x}}^{+}&=(1-\alpha\beta)\hat{\boldsymbol{x}}+\alpha\beta\mathrm{proj}_{C}\left(\hat{\boldsymbol{x}}-\gamma\hat{\boldsymbol{\xi}}\right)\\\
\hat{\boldsymbol{\xi}}^{+}&=(1-\alpha)\hat{\boldsymbol{\xi}}+\alpha
F(\hat{x})\\\
\hat{\boldsymbol{\mu}}^{+}&=\mathcal{R}\hat{\boldsymbol{\mu}}\end{array}\right..$
(85)

Furthermore, we use $\alpha\beta$ for the parameter of the second-time layer
separation to determine the second boundary layer system

$\displaystyle\left\\{\begin{array}[]{ll}{\boldsymbol{x}_{\textup{bl}}^{2}}^{+}&=\boldsymbol{x}_{\textup{bl}}^{2}\\\
{\boldsymbol{\xi}_{\textup{bl}}^{2}}^{+}&=(1-\alpha)\boldsymbol{\xi}_{\textup{bl}}^{2}+\alpha
F(\boldsymbol{x}_{\textup{bl}}^{2})\\\
{\boldsymbol{\mu}_{\textup{bl}}^{2}}^{+}&=\mathcal{R}{\boldsymbol{\mu}_{\textup{bl}}^{2}}\end{array}\right.,$
(89)

and the second averaged system

$\displaystyle\left\\{\begin{array}[]{ll}\tilde{\boldsymbol{x}}^{+}&=(1-\alpha\beta)\tilde{\boldsymbol{x}}+\alpha\beta\mathrm{proj}_{C}\left(\tilde{\boldsymbol{x}}-\gamma
F(\tilde{\boldsymbol{x}})\right)\\\
\tilde{\boldsymbol{\xi}}^{+}&=(1-\alpha)\tilde{\boldsymbol{\xi}}+\alpha
F(\tilde{\boldsymbol{x}})\\\
\tilde{\boldsymbol{\mu}}^{+}&=\mathcal{R}\tilde{\boldsymbol{\mu}}\end{array}\right.,$
(93)

which is the algorithm in (20) with additional bounded dynamics that renders
the set $\\{\boldsymbol{x}^{*}\\}\times C_{\textup{F}}\times\mathbb{S}^{m}$
UGAS.
In order to satisfy Assumption 1 for both averaged systems, we establish the
following result:

######  Lemma 4.

For any solution of the first boundary layer system
$(\boldsymbol{x}_{\textup{bl}}^{1},\boldsymbol{\xi}_{\textup{bl}}^{1},\boldsymbol{\mu}_{\textup{bl}}^{1})$
and compact set $C$ such that $\boldsymbol{x}_{\textup{bl}}^{1}\in C$ for all
$k\in\mathrm{dom}\left(\boldsymbol{x}_{\textup{bl}}^{1}\right)$, it holds
that:

$\displaystyle\Biggl{\|}\frac{1}{N}\sum_{i=1}^{N}\Bigl{[}{2}{A^{-1}}J\left(\boldsymbol{x}_{\textup{bl}}^{1}(i)+A\mathbb{D}\boldsymbol{\mu}_{\textup{bl}}^{1}(i)\right)\mathbb{D}\boldsymbol{\mu}_{\textup{bl}}^{1}(i)$
$\displaystyle-F(\boldsymbol{x}_{\textup{bl}}^{1}(i))-\mathcal{O}(\bar{a})\Bigr{]}\Biggr{\|}\leq\sigma_{1}\left(N\right),$
(94)

where $\sigma_{1}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ is a function of
class $\mathcal{L}$. ∎

See Appendix C. $\blacksquare$

######  Lemma 5.

For any solution of the second boundary layer system
$(\boldsymbol{x}_{\textup{bl}}^{2},\boldsymbol{\xi}_{\textup{bl}}^{2})$ and
compact set $C$ such that
$\operatorname{col}\left({\boldsymbol{x}_{\textup{bl}}^{1},\boldsymbol{\xi}_{\textup{bl}}^{1}}\right)\in
C$ for all $k\in\mathrm{dom}\left(\boldsymbol{x}_{\textup{bl}}^{1}\right)$, it
holds that:

$\displaystyle\left\|\frac{\gamma}{N}\sum_{i=1}^{N}\left[\boldsymbol{\xi}_{\textup{bl}}^{2}(i)-F(\boldsymbol{x}_{\textup{bl}}^{2}(i))\right]\right\|\leq\sigma_{2}\left(N\right),$
(95)

where $\sigma_{2}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ is a function of
class $\mathcal{L}$. ∎

See Appendix D. $\blacksquare$To prove stability, we start from the second
layer and “move” upwards. As the second averaged system satisfies Assumptions
1 due to Lemma 4, Assumption 2 due the nonexpansivnes of the projection
mapping [3, Prop. 12.28, 29.1], Assumption 3 due to Lemma 2 and Assumption 4
due to (21), we have that due to Theorem 1, the nominal averaged system in
(85) renders the set $\\{x^{*}\\}\times C_{\textup{F}}\times\mathbb{S}^{m}$
SGPAS as $\alpha\beta\rightarrow 0$, with the Lyapunov difference given by

$\displaystyle
V(\tilde{\boldsymbol{x}}^{+}+\eta_{1}^{+})-V(\tilde{\boldsymbol{x}}+\eta_{1})\leq$
$\displaystyle\quad-\tfrac{1}{2}\alpha\beta(1-c)(2-\alpha\beta
c)\left\|\tilde{\boldsymbol{x}}+\eta_{1}-\boldsymbol{x}^{*}\right\|_{\mathcal{A}}$
$\displaystyle\text{for
}\left\|\tilde{\boldsymbol{x}}+\eta_{1}\right\|_{\mathcal{A}}\geq\alpha_{\varepsilon}(\alpha\beta)\geq\sqrt{\tfrac{\alpha\beta\,L\left\|\eta_{1}\right\|}{2\alpha\beta(1-c)(2-\alpha\beta
c)}}$ (96)

and the perturbation dynamics

$\displaystyle\eta_{1}^{+}=(1-\alpha\beta)\eta_{1}+\alpha\beta(\mathrm{proj}_{C}\left(\tilde{\boldsymbol{x}}-\gamma
F(\tilde{\boldsymbol{x}})\right)-\tilde{\boldsymbol{\xi}}).$ (97)

We note that we had to take $\alpha\beta$ as the time-scale separation
parameter. If we had chosen only $\beta$, as might be the intuition, the
function $\alpha_{\varepsilon}$ of class $\mathcal{K}$ that appears in the
inequality in (96), would have an implicit dependence on the parameter
$\alpha$. In fact, decreasing $\alpha$ would increase the value of the
function $\alpha_{\varepsilon}$, as it would hold

$\displaystyle\alpha_{\varepsilon}(\beta)\geq\sqrt{\tfrac{\beta\,L\left\|\eta_{1}\right\|}{2\alpha\beta(1-c)(2-\alpha\beta
c)}},$ (98)

which would invalidate all of the following stability analysis. Thus, it is
important to capture all of the parameters that affect the speed of
convergence. Nevertheless, if we assume that the parameter $\alpha$ is
contained in the set $(0,\overline{\alpha})$, it is possible to construct a
function of class $\mathcal{K}$, such that it holds
$\left\|\tilde{x}+\eta_{1}\right\|_{\mathcal{A}}\geq\alpha_{\beta}(\beta)\geq\alpha_{\varepsilon}(\alpha\beta)$.
Hence, the averaged system in (85) renders the set $\\{x^{*}\\}\times
C_{\textup{F}}\times\mathbb{S}^{m}$ SGPAS as $\beta\rightarrow 0$.
The first averaged system is an $\mathcal{O}(\overline{a})$ inflation of the
nominal averaged system, and it can be shown that the inflation introduces a
small perturbation in the Lyapunov difference inequality which can be made
arbitrarily small by choosing $\overline{a}$ small enough. For the sake of the
proof, we set $\overline{a}=\alpha_{a}(\beta)$, where $\alpha_{a}$ is a
function of class $\mathcal{K}$. Thus, it also satisfies Assumption 4,
Assumption 1 due to (97), Assumption 2 because of [3, Prop. 12.28, 29.1], and
Assumption 3 as a result of (97). Hence, the system in (25) renders the set
$\\{x^{*}\\}\times C_{\textup{F}}\times\mathbb{S}^{m}$ SGPAS as
$(\alpha,\overline{a},\beta)\rightarrow 0$. $\blacksquare$

##  Appendix C Proof of Lemma 4

Without the loss of generality, let a solution of the first boundary layer
system is given by
$\boldsymbol{x}_{\textup{bl}}^{1}(k)=\boldsymbol{x}_{\textup{bl}}^{1}=const.$,
$\boldsymbol{\xi}_{\textup{bl}}^{1}=const.$,
$\mu_{\textup{bl}}^{1,i}(k)=\operatorname{col}\left({(\sin(\omega_{i}^{j}k),\cos(\omega_{i}^{j}k))_{j\leq
m_{i}}}\right)$,
$\mu_{\textup{bl}}^{1}(k)=\operatorname{col}\left({(\mu_{\textup{bl}}^{1,i}(k)))_{i\in\mathcal{I}}}\right)$.
First, with the following Lemma, we characterize the properties of average
discrete-time sinusoidal signals:

######  Lemma 6.

For any $\phi,\phi_{i},\phi_{l}\in\mathbb{R}$ such that $\phi\neq 2\pi t$,
$\phi_{i}\pm\phi_{l}\neq 2\pi p,t,p\in\mathbb{Z}$, it holds that:

$\displaystyle\frac{1}{N}\left\|\sum_{k=0}^{N-1}\sin(\phi
k)\right\|\leq\frac{c_{1}}{N},$ (99)
$\displaystyle\frac{1}{N}\left\|\sum_{k=0}^{N-1}\cos(\phi
k)\right\|\leq\frac{c_{2}}{N},$ (100)
$\displaystyle\frac{1}{N}\left\|\sum_{k=0}^{N-1}\sin(\phi_{i}k)\sin(\phi_{l}k)\right\|\leq\frac{c_{3}}{N},$
(101) $\displaystyle\frac{1}{N}\left\|\sum_{k=0}^{N-1}(\sin^{2}(\phi
k)-\tfrac{1}{2})\right\|\leq\frac{c_{4}}{N},$ (102)

for some $c_{1},c_{2},c_{3},c_{4}>0$. ∎

We note that for $\phi\neq 2\pi p,p\in\mathbb{Z}$, it follows:

$\displaystyle\left\|\sum_{k=0}^{N-1}e^{j\phi
k}\right\|^{2}=\left\|\sum_{k=0}^{N-1}\cos(\phi
k)\right\|^{2}+\left\|\sum_{k=0}^{N-1}\sin(\phi k)\right\|^{2}$
$\displaystyle=\left\|\frac{e^{j\phi
N}-1}{e^{j\phi}-1}\right\|^{2}\leq\frac{4}{\left\|e^{j\phi}-1\right\|^{2}}.$
(103)

Therefore, we have

$\displaystyle\left\|\sum_{k=0}^{N-1}\cos(\phi
k)\right\|\leq\frac{2}{\left\|e^{j\phi}-1\right\|}\coloneqq c_{1},$ (104)
$\displaystyle\left\|\sum_{k=0}^{N-1}\sin(\phi
k)\right\|\leq\frac{2}{\left\|e^{j\phi}-1\right\|}\coloneqq c_{2}.$ (105)

Equations (99), (100) follow from the previous equations.
Let $\phi=\phi_{i}\pm\phi_{l}$. From (104), we have

$\displaystyle\left\|\sum_{k=0}^{N-1}\left(\sin(\phi_{i}k)\sin(\phi_{l}k)\mp\cos(\phi_{i}k)\cos(\phi_{l}k)\right)\right\|$
$\displaystyle\leq\frac{2}{\min_{\phi\in\phi_{i}\pm\phi_{l}}\left\|e^{j\phi}-1\right\|}\coloneqq
c_{3}$ (106)

For any scalars $a,b\in\mathbb{R},c\in\mathbb{R}_{+}$, that satisfy equation
$\left\|a\pm b\right\|\leq c$, it holds

$\displaystyle-c\leq a+b$ $\displaystyle\leq c\text{, and}$
$\displaystyle-c\leq a-b$ $\displaystyle\leq c.$ (107)

By summing the last two inequalities, we have

$\displaystyle\left\|a\right\|\leq c.$ (108)

Thus from (106), (107) and (108), we conclude

$\displaystyle\left\|\sum_{k=0}^{N-1}{\sin(\phi_{i}k)\sin(\phi_{l}k)}\right\|\leq
c_{3}.$ (109)

Again, (101) follows trivially. Finally, using the identity
$\sin^{2}(x)=\frac{1-\cos(2x)}{2}$, we rewrite Equation (102) as

$\displaystyle\frac{1}{2N}\left\|\sum_{k=0}^{N-1}\cos(2\phi
k)\right\|\leq\frac{c_{4}}{N}.$ (110)

By switching $2\phi$ instead of $\phi$ in (104), analogously it is possible to
prove (102). $\blacksquare$Via the Taylor expansion of the an addend in (94),
we have

$\displaystyle{2}{A^{-1}}J\left(\boldsymbol{x}_{\textup{bl}}^{1}+A\mathbb{D}\mu_{\textup{bl}}^{1}(i)\right)\mathbb{D}\mu_{\textup{bl}}^{1}(i)-F(\boldsymbol{x}_{\textup{bl}}^{1})-\mathcal{O}(\bar{a})$
$\displaystyle={{2}A^{-1}}J(\boldsymbol{x}_{\textup{bl}}^{1})\mathbb{D}\mu_{\textup{bl}}^{1}(i)+F(\boldsymbol{x}_{\textup{bl}}^{1})^{\top}\mathbb{D}\mu_{\textup{bl}}^{1}(i)\mathbb{D}\mu_{\textup{bl}}^{1}(i)$
$\displaystyle-F(\boldsymbol{x}_{\textup{bl}}^{1})$ (111)

Due to the inequality
$\sqrt{\sum_{i=0}^{m}x_{i}^{2}}\leq\sum_{i=0}^{m}|x_{i}|$, we can bound the
expression in (94) via the bounds for each row:

$\displaystyle\Biggl{\|}\frac{1}{N}\sum_{i=1}^{N}\Bigl{[}{2}{A^{-1}}J\left(\boldsymbol{x}_{\textup{bl}}^{1}+A\mathbb{D}\mu_{\textup{bl}}^{1}(i)\right)\mathbb{D}\mu_{\textup{bl}}^{1}(i)$
$\displaystyle-F(\boldsymbol{x}_{\textup{bl}}^{1})-\mathcal{O}(\bar{a})\Bigr{]}\Biggr{\|}\leq$
$\displaystyle\sum_{j=1}^{m}\left\|\sum_{k=0}^{N-1}\left[\tfrac{2}{a_{j}}J_{j}(\boldsymbol{x}_{\textup{bl}}^{1})\sin(\phi_{j}k)\right.\right.$
$\displaystyle+\nabla_{x_{j}}J_{j}(\boldsymbol{x}_{\textup{bl}}^{1})(2\sin^{2}(\phi_{j}k)-1)$
$\displaystyle+2\sum_{l\neq
j}^{n}\nabla_{x_{l}}J_{j}(\boldsymbol{x}_{\textup{bl}}^{1})\sin(\phi_{l}k)\sin(\phi_{j}k)\Bigr{]}\Biggr{\|}\leq$
$\displaystyle\frac{2m}{\underline{a}}\left\|J(\boldsymbol{x}_{\textup{bl}}^{1})\right\|\frac{c_{1}}{N}+\left\|\nabla
J(\boldsymbol{x}_{\textup{bl}}^{1})\right\|_{\infty}\frac{2\overline{c}\,m^{2}}{N},$
(112)

where $\underline{a}\coloneqq\min_{i}a_{i}$ and
$\overline{c}\coloneqq\max\\{c_{3},c_{4}\\}$. Thus for the compact set $C$, we
define $\sigma_{1}(N)\coloneqq\frac{2m}{\underline{a}}\max_{x\in
C}\left\|J(x)\right\|\frac{c_{1}}{N}+\max_{x\in C}\left\|\nabla
J(x)\right\|_{\infty}\frac{2\overline{c}\,m^{2}}{N}$ which belongs to the
class of $\mathcal{L}$ functions. $\blacksquare$

##  Appendix D Proof of Lemma 5

A solutions of the second boundary layer system are given by
$\boldsymbol{x}_{\textup{bl}}^{2}(k)=\boldsymbol{x}_{\textup{bl}}^{2}=const.$,
and
$\boldsymbol{\xi}_{\textup{bl}}^{2}(k)=(1-\alpha)^{k}\left(\boldsymbol{\xi}_{\textup{bl}}^{2}(0)-F(\boldsymbol{x}_{\textup{bl}}^{2})\right)+F(\boldsymbol{x}_{\textup{bl}}^{2})$.
Thus, the norm in (95) can be rewritten as

$\displaystyle\left\|\frac{\gamma}{N}\sum_{i=0}^{N-1}\left[\boldsymbol{\xi}_{\textup{bl}}^{2}(i)-F(\boldsymbol{x}_{\textup{bl}}^{2}(i))\right]\right\|\leq$
$\displaystyle\frac{\gamma}{N}\left\|\sum_{i=0}^{N-1}\left[(1-\alpha)^{i}\left(\boldsymbol{\xi}_{\textup{bl}}^{2}(0)-F(\boldsymbol{x}_{\textup{bl}}^{2}\right))\right]\right\|\leq$
$\displaystyle\frac{\gamma}{N\alpha}\left\|{{\boldsymbol{\xi}_{\textup{bl}}^{2}(0)-F(\boldsymbol{x}_{\textup{bl}}^{2})}}\right\|.$
(113)

Thus, for the compact set $C$, we define
$\sigma_{2}(N)\coloneqq\frac{\gamma}{N\alpha}\max_{\operatorname{col}\left({\boldsymbol{\xi},\boldsymbol{x}}\right)\in
C}\left\|{{\boldsymbol{\xi}-F(\boldsymbol{x})}}\right\|$, which belongs to the
class of $\mathcal{L}$ functions. $\blacksquare$

##  Appendix E Proof of Theorem 3

For notational simplicity, we denote
$S(k)\coloneqq\hat{S}_{x}\left(k;\boldsymbol{\tau}_{0}\right)$. Thus, the
algorithm reads as

$\displaystyle\left\\{\begin{array}[]{ll}\boldsymbol{x}^{+}&=\boldsymbol{x}-\alpha{S}\left(k\right)F(\boldsymbol{x})\\\
k^{+}&=k+1.\end{array}\right.$ (116)

One epoch is defined as $r$ iterations of the algorithm in (116), where $r$ is
the period of the function $S$ from Lemma 1. From the proof of the Lemma, it
follows that every agent _individually_ jumps $r_{i}$ times in one epoch. Let
$F_{j}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ be the mapping that returns
the rows of the pseudogradient that correspond to the agents that sample at
$k=j+rn$, $n\in\mathbb{N}$, $0\leq j\leq r-1$, i.e.
$F_{j}(\boldsymbol{x})\coloneqq S(j)F(\boldsymbol{x})$. We define the full
update operator, the asynchronous update operator, and error operator
respectively, as

$\displaystyle T\coloneqq I-\alpha\Gamma F,$ (117) $\displaystyle
E\coloneqq(I-\alpha F_{1})\dots(I-\alpha F_{r}),$ (118) $\displaystyle
R\coloneqq\tfrac{1}{\alpha}(T-E).$ (119)

where
$\Gamma\coloneqq\operatorname{diag}\left({(r_{i}I_{m_{i}})_{i\in\mathcal{I}}}\right)$.
We note that the $E$ operator represents one epoch of the algorithm in (116),
i.e. $\boldsymbol{x}(k+r)=E(x(k))$.
The proof of convergence is analogous to the proof in [5], and here we just
provide the outlines. The error operator can be bounded as

$\displaystyle\left\|R(x)\right\|^{2}\leq\frac{\alpha^{2}L^{2}r^{2}\overline{r}^{2}}{2}(1+\alpha
L)^{2m}\left\|S(x)\right\|^{2},$ (120)

where $\overline{r}=\max_{i\in\mathcal{I}}r_{i}$. For the Lyapunov function
candidate, we propose

$\displaystyle
V(\boldsymbol{x})=\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}.$
(121)

It can be proven that

$\displaystyle\left\|E(\boldsymbol{x})-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}\leq\left(1-\frac{\alpha\mu_{\textup{F}}^{2}\underline{r}^{2}}{2}\right)\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2},$
(122)

where $\underline{r}=\min_{i\in\mathcal{I}}r_{i}$, if $\alpha$ is chosen such
that

$\displaystyle\tfrac{1}{2}+\tfrac{\overline{r}\eta}{\mu_{\textup{F}}^{2}}-(1-\alpha\eta)\left(1-\alpha-\tfrac{\left(\alpha
Lr\overline{r}(1+\alpha L)^{r}\right)^{2}}{2\eta}\right)\leq 0.$

The inequality is satisfied for $\eta$ and $\alpha$ small enough, and
$\alpha\ll\eta$. We note that the inequality does not depend on the initial
conditions of the timers $\boldsymbol{\tau}_{0}$. Equation (122) holds for
_epochs_ , not necessarily the individual samples. Due to Lipschitz continuity
of the pseudogradient, it follows that

$\displaystyle\left\|\boldsymbol{x}^{+}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}\leq(1+\alpha
L\tfrac{\overline{r}}{\underline{r}})\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}.$

Thus, for some $k=ir+j$, where $i,j\in\mathbb{N}$ we have

$\displaystyle\left\|\boldsymbol{x}(k)-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}\leq$
$\displaystyle\left(1+\alpha
L\tfrac{\overline{r}}{\underline{r}}\right)^{2r}\left(1-\frac{\alpha\mu_{\textup{F}}^{2}\underline{r}^{2}}{2}\right)^{i}\left\|\boldsymbol{x}(0)-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}.$
(123)

It holds $i=\lfloor\tfrac{k}{r}\rfloor\geq\tfrac{k}{r}-1$. Hence, the previous
inequality becomes

$\displaystyle\left\|\boldsymbol{x}(k)-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}\leq$
$\displaystyle\left(1+\alpha
L\tfrac{\overline{r}}{\underline{r}}\right)^{2r}\left(1-\frac{\alpha\mu_{\textup{F}}^{2}\underline{r}^{2}}{2}\right)^{\tfrac{k}{r}-1}\left\|\boldsymbol{x}(0)-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}.$

The last inequality is the KL exponential stability bound for all initial
conditions of timer states $\boldsymbol{\tau}_{0}$. Thus, the dynamics in
(116) render $\boldsymbol{x}^{*}$ UGES.

Additionally, we need to establish the Lyapunov difference convergence speed.
To do this, we construct a Lyapunov function using a similar procedure as in
the proof of [9, Thm. 4.14], which we omit due to space limitations. Let

$\displaystyle a$ $\displaystyle\coloneqq\frac{\left(1+\alpha
L\tfrac{\overline{r}}{\underline{r}}\right)^{2r}}{1-\frac{\alpha\mu_{\textup{F}}^{2}\underline{r}^{2}}{2}}$
$\displaystyle b$
$\displaystyle\coloneqq\left(1-\frac{\alpha\mu_{\textup{F}}^{2}\underline{r}^{2}}{2}\right)^{\tfrac{1}{r}}.$
(124)

Then, the Lyapunov function satisfies the following properties

$\displaystyle\frac{1-\left(1-\alpha
L\tfrac{\overline{r}}{\underline{r}}\right)^{2\delta\,r}}{2\alpha
L\tfrac{\overline{r}}{\underline{r}}(1-\alpha
L\tfrac{\overline{r}}{\underline{r}})}\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}\leq
V(\boldsymbol{x})$
$\displaystyle\leq\frac{a(b^{\delta}\,r-1)}{b-1}\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2}$
$\displaystyle
V(\boldsymbol{x}^{+})-V(\boldsymbol{x})\leq-(1-ab^{\delta\,r})\left\|\boldsymbol{x}-\boldsymbol{x}^{*}\right\|_{\Gamma^{-1}}^{2},$

where $\delta$ is large integer. Using the Taylor expansion, for small values
of $\alpha$, it holds

$\displaystyle(1-ab^{\delta\,r})\approx\alpha\left(\frac{\mu^{2}\underline{r}^{2}(\delta-1)}{2}-2rL\tfrac{\overline{r}}{\underline{r}}\right).$
(125)

Thus, for $\delta$ large enough, we can guarantee that the Lyapunov difference
is negative. Furthermore, we have that $\tfrac{\alpha}{(1-ab^{\delta\,r})}$ is
bounded on some interval $(0,\alpha^{*})$ and the Lyapunov function satisfies
the conditions from Assumption (4). $\blacksquare$

##  Appendix F Proof of Lemma 1

First, let us denote least common sampling time as $T=T_{i}\tfrac{p}{p_{i}}$.
_Claim 1._ The number of jumps in any time interval $[\hat{t},\hat{t}+T)$ is
constant. ∎
Let us denote as the number of jumps in this interval as $q$. We “slide” the
interval by some $\Delta T$, i.e. $[\hat{t}+\Delta T,\hat{t}+T+\Delta T)$, so
that we exclude one event in $[\hat{t},\hat{t}+\Delta T)$. As $T=T_{i}r_{i}$
for all $i\in\mathcal{I}$, it follows that there must be an event in the
interval $[\hat{t}+T,\hat{t}+T+\Delta T)$, thus the total number of jumps in
the interval $[\hat{t}+\Delta T,\hat{t}+T+\Delta T)$ remains the same. We can
repeat this procedure for any $\Delta\hat{T}$ by sliding the interval for
every jump by $\Delta T_{j}$ until $\Delta\hat{T}=\sum_{j}\Delta T_{j}$. The
number of jumps in the interval is equal to $q$.
_Claim 2._ For every $(\hat{t},\hat{j})\in\mathrm{dom}(\boldsymbol{\tau})$
where a jump occurred, it holds
$\boldsymbol{\tau}(\hat{t},\hat{j})=\boldsymbol{\tau}(\hat{t}+T,\hat{j}+r)$. ∎
We observe that if a jump is initiated by an agent at $t=\hat{t}$, these same
agent will also initiate a jump at $t=\hat{t}+T$. Furthermore, if an agents
did not jump at $\hat{t}$, it will also not jump at $\hat{t}+T$. Thus, in the
moment $t=\hat{t}+T$, the same agent will jump, i.e. $\tau_{i}=T_{i}$ and it
follows that
$\boldsymbol{\tau}(\hat{t},\hat{j})=\boldsymbol{\tau}(\hat{t}+T,\hat{j}+r)$.
As the functions $S_{x}$ and $S_{\tau}$ are single-valued, the claim of the
Lemma holds. $\blacksquare$

##  Appendix G Proof of Theorem 4

As the proof is analogous to the proof of Theorem 2, we provide just the
required system definitions and averaging Lemmas. The equivalent discrete-time
system of (33) is given by

$\displaystyle\left\\{\begin{array}[]{ll}\boldsymbol{x}^{+}&=\boldsymbol{x}-\alpha\beta
S(k)\boldsymbol{\xi}\\\ \boldsymbol{\xi}^{+}&=\boldsymbol{\xi}+\alpha
S(k)\left(2A^{-1}J(\boldsymbol{x}+A\mathbb{D}\boldsymbol{\mu})\mathbb{D}\boldsymbol{\mu}-\boldsymbol{\xi}\right)\\\
\boldsymbol{\mu}^{+}&=(I-S_{\mu}(k))\boldsymbol{\mu}+S_{\mu}(k))\mathcal{R}\boldsymbol{\mu}\\\
\boldsymbol{\kappa}^{+}&=\boldsymbol{\kappa}+S_{\tau}(\boldsymbol{\tau})\\\
k^{+}&=k+1.\end{array}\right.$ (131)

The first boundary-layer system is defined as

$\displaystyle\left\\{\begin{array}[]{ll}{\boldsymbol{x}^{1}_{\textup{bl}}}^{+}&=\boldsymbol{x}^{1}_{\textup{bl}}\\\
{\boldsymbol{\xi}^{1}_{\textup{bl}}}^{+}&=\boldsymbol{\xi}^{1}_{\textup{bl}}\\\
{\boldsymbol{\mu}^{1}_{\textup{bl}}}^{+}&=(I-S_{\mu}(k))\boldsymbol{\mu}^{1}_{\textup{bl}}+S_{\mu}(k))\mathcal{R}\boldsymbol{\mu}^{1}_{\textup{bl}}\\\
{\boldsymbol{\kappa}}^{+}&=\boldsymbol{\kappa}+S_{\tau}(k)\\\
k^{+}&=k+1,\end{array}\right.$ (137)

while the first averaged system is given as

$\displaystyle\left\\{\begin{array}[]{ll}\tilde{\boldsymbol{x}}^{+}&=\tilde{\boldsymbol{x}}-\alpha\beta
S(k)\tilde{\boldsymbol{\xi}}\\\
\tilde{\boldsymbol{\xi}}^{+}&=\tilde{\boldsymbol{\xi}}+\alpha
S(k)\left(F(\tilde{\boldsymbol{x}})-\tilde{\boldsymbol{\xi}}\right)\\\
\tilde{\boldsymbol{\mu}}^{+}&=(I-S_{\mu}(k))\tilde{\boldsymbol{\mu}}+S_{\mu}(k))\mathcal{R}\tilde{\boldsymbol{\mu}}\\\
{\boldsymbol{\kappa}}^{+}&=\boldsymbol{\kappa}+S_{\tau}(k)\\\
k^{+}&=k+1.\end{array}\right.$ (143)

The second boundary-layer system follows the dynamics

$\displaystyle\left\\{\begin{array}[]{ll}{\boldsymbol{x}^{2}_{\textup{bl}}}^{+}&=\boldsymbol{x}^{2}_{\textup{bl}}\\\
{\boldsymbol{\xi}^{2}_{\textup{bl}}}^{+}&=\boldsymbol{\xi}^{2}_{\textup{bl}}-\alpha
S(k)(F(\boldsymbol{x}^{2}_{\textup{bl}})-\boldsymbol{\xi}^{2}_{\textup{bl}})\\\
{\boldsymbol{\mu}^{2}_{\textup{bl}}}^{+}&=(I-S_{\mu}(k))\boldsymbol{\mu}^{2}_{\textup{bl}}+S_{\mu}(k))\mathcal{R}\boldsymbol{\mu}^{2}_{\textup{bl}}\\\
{\boldsymbol{\kappa}}^{+}&=\boldsymbol{\kappa}+S_{\tau}(k)\\\
k^{+}&=k+1,\end{array}\right.$ (149)

whereas the second averaged system is defined as

$\displaystyle\left\\{\begin{array}[]{ll}\tilde{\boldsymbol{x}}^{+}&=\tilde{\boldsymbol{x}}-\alpha\beta
S(k)F(\tilde{\boldsymbol{x}})\\\
\tilde{\boldsymbol{\xi}}^{+}&=\tilde{\boldsymbol{\xi}}+\alpha
S(k)\left(F(\tilde{\boldsymbol{x}})-\tilde{\boldsymbol{\xi}}\right)\\\
\tilde{\boldsymbol{\mu}}^{+}&=(I-S_{\mu}(k))\tilde{\boldsymbol{\mu}}+S_{\mu}(k))\mathcal{R}\tilde{\boldsymbol{\mu}}\\\
{\boldsymbol{\kappa}}^{+}&=\boldsymbol{\kappa}+S_{\tau}(k)\\\
k^{+}&=k+1.\end{array}\right.$ (155)

To prove Assumption 1, the following two Lemmas are needed:

######  Lemma 7.

For any solution of the first boundary layer system
$(\boldsymbol{x}_{\textup{bl}}^{1},\boldsymbol{\xi}_{\textup{bl}}^{1},\boldsymbol{\mu}_{\textup{bl}}^{1})$
and compact set $C$, such that $\boldsymbol{x}_{\textup{bl}}^{1}\in C$ for all
$k\in\mathrm{dom}\left(\boldsymbol{x}_{\textup{bl}}^{1}\right)$, it holds
that:

$\displaystyle\Biggl{\|}\frac{1}{N}\sum_{i=1}^{N}\Bigl{[}{2}{A^{-1}}J\left(\boldsymbol{x}_{\textup{bl}}^{1}(i)+A\mathbb{D}\boldsymbol{\mu}_{\textup{bl}}^{1}(i)\right)\mathbb{D}\boldsymbol{\mu}_{\textup{bl}}^{1}(i)$
$\displaystyle-F(\boldsymbol{x}_{\textup{bl}}^{1}(i))-\mathcal{O}(\bar{a})\Bigr{]}S(i)\Biggr{\|}\leq\sigma_{1}\left(N\right),$
(156)

where $\sigma_{1}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ is a function of
class $\mathcal{L}$. ∎

First, we note that the difference in the previous inequality for rows
corresponding to agent $i$ is equal to zero whenever the agent is not jumping.
This motivates us to study the “isolated” system of agent $i$, instead of the
group dynamics in (137). Consider the dynamics

$\displaystyle\mu^{i}(\kappa_{i}+1)=\mathcal{R}_{i}\mu^{i}(\kappa_{i}).$ (157)

Without the loss of generality, let a solution of the previous system be given
by $\mu^{i}(\kappa_{i})=\operatorname{col}((\cos(\omega_{i}^{j}\kappa_{i}),$
$\sin(\omega_{i}^{j}\kappa_{i}))_{j\leq m_{i}})$. The solution of
$\mu_{\textup{bl}}^{1,i}$ is similar to the solution of $\mu^{i}$, as it also
has the same samples, but they “persist” for more iterations, i.e. until the
agent $i$ jumps again. If we define the set valued mapping
$k(v,i)\coloneqq\\{u\mid\kappa_{i}(u)=v\\}$ that relates the global jump
counter to the internal counter of agent $i$, it holds that
$\mu_{\textup{bl}}^{1,i}(k(v,i))=\mu^{i}(v)$. Furthermore, at the sample $v$
given by the internal counter for agent $i$, or given by
$\overline{k}(v,i)=\max k(v,i)$ by the global counter, the internal counter of
some other agent $j$ is given by

$\displaystyle\kappa_{j}(v,i)\coloneqq\left\lfloor\Delta_{j}^{i}+\tfrac{T_{i}}{T_{j}}v\right\rfloor.$
(158)

where
$\Delta_{j}^{i}\coloneqq\kappa_{j}^{0}+\tfrac{\tau_{i}-\kappa_{i}^{0}T_{i}}{T_{j}}$.
Hence it holds
$\mu_{\textup{bl}}^{1,j}(\overline{k}(v,i))=\mu^{j}(\kappa_{j}(v,i))$. Lastly,
we see that diagonal elements of $S(k)$ corresponding to agent $i$ are
different from zero for $k\in\\{u\mid u=\bar{k}(v,i),v\in\mathbb{N}\\}$. Let
us denote the norm in (156) as $\Phi$ and write $N=r\,l+o,l\in\mathbb{N},0\leq
o<r$. The previous iterator relations and properties of $S(k)$ allows us bound
the inequality as

$\displaystyle
N\Phi\leq\sum_{i=1}^{M}\Biggl{\|}\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}\,l-1}\Bigl{[}-F(\boldsymbol{x}_{\textup{bl}})-\mathcal{O}(\bar{a})+$
$\displaystyle\tfrac{2}{a_{i}}J_{i}\left(\boldsymbol{x}_{\textup{bl}}^{1}+A\mathbb{D}\boldsymbol{\mu}_{\textup{bl}}^{1}(\overline{k}(v,i)\right)\mathbb{D}_{i}{\mu}_{\textup{bl}}^{1,i}(\overline{k}(v,i))\Bigr{]}\Biggr{\|}+\Phi_{\textup{r}}$
$\displaystyle\leq\sum_{i=1}^{M}\Biggl{\|}\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}\,l-1}\Bigl{[}\tfrac{2}{a_{i}}J_{i}\left(\boldsymbol{x}_{\textup{bl}}^{1}\right)\mathbb{D}_{i}\mu_{i}(v)$
$\displaystyle+2\nabla_{x_{i}}J_{i}(\boldsymbol{x}_{\textup{bl}}^{1})^{\top}\mathbb{D}_{i}\mu_{i}(v)\mathbb{D}_{i}\mu_{i}(v)-\nabla_{x_{i}}J_{i}(\boldsymbol{x}_{\textup{bl}}^{1})$
$\displaystyle+2\sum_{j\neq
i}^{M}\nabla_{x_{j}}J_{i}(\boldsymbol{x}_{\textup{bl}}^{1})^{\top}\mathbb{D}_{j}\mu_{i}(\kappa_{j}(v,i))\mathbb{D}_{i}\mu_{i}(v)\Bigr{]}\Biggr{\|}+\Phi_{\textup{r}}$
$\displaystyle\leq\sum_{i=1}^{M}\sum_{u=1}^{m_{i}}\Biggl{\|}\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}\,l-1}\Bigl{[}\tfrac{2}{a_{i}}J_{i}\left(\boldsymbol{x}_{\textup{bl}}^{1}\right)\sin\left(\omega_{i}^{u}v\right)$
$\displaystyle\nabla_{x_{i}^{u}}J_{i}(\boldsymbol{x}_{\textup{bl}}^{1})\left[2\sin\left(\omega_{i}^{u}v\right)\sum_{j=1}^{m_{i}}\sin(\omega_{i}^{j}v)-1\right]$
$\displaystyle+2\sum_{j\neq
i}^{M}\sum_{s=1}^{m_{j}}\nabla_{x_{j}^{s}}J_{i}(\boldsymbol{x}_{\textup{bl}}^{1})\sin(\omega_{j}^{s}\kappa_{j}(v,i))\sum_{j=1}^{m_{i}}\sin(\omega_{i}^{j}v)\Bigr{]}\Biggr{\|}$
$\displaystyle+\Phi_{\textup{r}},$ (159)

where

$\displaystyle\Phi_{\textup{r}}\coloneqq\sum_{i=1}^{M}\,r_{i}\,m\,M\,\max_{\boldsymbol{x}\in
C}\left\|\tfrac{1}{a_{i}}J\left(\boldsymbol{x}+A\mathbb{B}\right)\right\|+\left\|F(\boldsymbol{x})+\mathcal{O}(\overline{a})\right\|.$
(160)

Using Lemma (5) and Assumption 11, we can derive the upper bounds of all of
the sums in the norm, apart for the last one, which contains addends of form
$\sin(\omega_{i}^{j}v)\sin(\omega_{j}^{s}\kappa_{j}(v,i))$. Using the same
procedure as in proof of Lemma 4 and Equation (158), we find the equivalent
exponential representation. For some $\omega_{1},\omega_{2}$, consider the sum

$\displaystyle\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}\,l-1}e^{j\left(\omega_{1}v+\omega_{2}\left\lfloor\Delta_{j}^{i}+\tfrac{p_{i}}{p_{j}}v\right\rfloor\right)}$
$\displaystyle=\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}-1}e^{j\left(\omega_{1}v+\omega_{2}\left\lfloor\Delta_{j}^{i}+\tfrac{p_{i}}{p_{j}}v\right\rfloor\right)}\times\sum_{u=0}^{l-1}e^{j\left(\omega_{1}r_{i}+\omega_{2}r_{j}\right)u}$
$\displaystyle=\frac{1-e^{j\left(\omega_{1}r_{i}+\omega_{2}r_{j}\right)l}}{1-e^{j\left(\omega_{1}r_{i}+\omega_{2}r_{j}\right)}}\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}-1}e^{j\left(\omega_{1}v+\omega_{2}\left\lfloor\Delta_{j}^{i}+\tfrac{p_{i}}{p_{j}}v\right\rfloor\right)},$

where the second equality follows from Assumption 9, the properties of the
least common multiple $p=p_{i}\,r_{i}=p_{j}\,r_{j}$, and last equality holds
for $\omega_{1}r_{i}+\omega_{2}r_{j}\neq 2\pi z,z\in\mathbb{Z}$. Thus, we have

$\displaystyle\left\|\sum_{v=\kappa_{i}(0)}^{\kappa_{i}(0)+r_{i}\,l-1}e^{j\left(\omega_{1}v+\omega_{2}\left\lfloor\Delta_{j}^{i}+\tfrac{p_{i}}{p_{j}}v\right\rfloor\right)}\right\|^{2}\leq$
$\displaystyle\left\|\frac{r_{i}}{1-e^{j\left(\omega_{1}r_{i}+\omega_{2}r_{j}\right)}}\right\|^{2}\leq
C_{e},$

where $C_{e}$ is the supremum with respect to all possible combinations of
$r_{i}$ and $\omega_{i}^{j}$. The rest of the procedure follows the same steps
as after Equation (106) in the proof of Lemma 5. The bound that holds
regardless of the initial conditions of the timers $\boldsymbol{\tau}(0)$.
Thus, the Lemma holds. $\blacksquare$

######  Lemma 8.

For any solution of the second boundary layer system
$(\boldsymbol{x}_{\textup{bl}}^{2},\boldsymbol{\xi}_{\textup{bl}}^{2})$, it
holds that:

$\displaystyle\left\|\frac{\gamma}{N}\sum_{i=1}^{N}S(i)\left[\boldsymbol{\xi}_{\textup{bl}}^{2}(i)-F(\boldsymbol{x}_{\textup{bl}}^{2}(i))\right]\right\|\leq\sigma_{2}\left(N\right),$
(161)

where $\sigma_{2}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ is a function of
class $\mathcal{L}$. ∎

The proof is analogous to the proof of Lemmas 5 and 7 and is omitted due to
space constraints. $\blacksquare$