# Bounded Extremum Seeking for Static Quadratic Maps using Nonlinear Transformation and Lyapunov Method

Frederic Mazenc, Michael Malisoff, Emilia Fridman

► To cite this version:

Frederic Mazenc, Michael Malisoff, Emilia Fridman. Bounded Extremum Seeking for Static Quadratic Maps using Nonlinear Transformation and Lyapunov Method. IEEE Transactions on Automatic Control, 2025, 70 (5), pp. 10.1109/TAC.2024.3504356. hal-05239785

HAL Id: hal-05239785

https://inria.hal.science/hal-05239785v1

Submitted on 4 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Bounded Extremum Seeking for Static Quadratic Maps using Nonlinear Transformation and Lyapunov Method

Frederic Mazenc, Michael Malisoff, Senior Member, IEEE, and Emilia Fridman, Fellow, IEEE

**Abstract**—We present a new practical stability analysis for a bounded gradient based extremum seeking problem for two variable static quadratic maps that contain a time-varying additive measurement uncertainty. Instead of using earlier averaging-based approaches, we introduce a new state transformation, a time-varying quadratic Lyapunov function, and a comparison principle to obtain essentially less conservative bounds on the dither frequency and on the ultimate bound of the estimation error compared with earlier results. Our numerical example illustrates the efficiency of the method.

**Index Terms**—Extremum seeking, uncertainty, time-varying

## I. INTRODUCTION

Extremum seeking is a central current research topic in control theory, because of its ability to provide model free, online, real time optimization methods to find extrema when objective functions that need to be maximized or minimized contain significant uncertainties [1]. While basic extremum seeking was used by Leblanc in 1922 [2], the first mathematical stability analysis for extremum seeking appears to be in averaging and singular perturbation approaches of M. Krstic and his collaborators, e.g., in [3]; see [4] for a history of extremum seeking. This motivated widespread use of extremum seeking, including in aerospace models and source seeking [5], [6]. More recent theoretical studies of extremum seeking include [5], [7]–[19]. For instance, [5], [14], [15] included bounded extremum seeking where unknown maps arise in arguments of a sine or cosine, providing bounds for update rates.

When using extremum seeking, one often encounters measurement uncertainties [1]. In the previous works on extremum seeking, one finds an input-to-state stability (or ISS) analysis. For instance, additive measurement uncertainty as we consider in this paper was considered in [19], [20]. A related problem arises when an extremum seeking algorithm is used to control a given system, but where the given system contains uncertainties, such as additive uncertainties on the coefficient matrices

in linear systems. Such approaches have been pursued using averaging [20]–[22] or Lyapunov function methods, by placing positive lower bounds on the dither periods in extremum seeking algorithms. See also the works [1], [6], [23]–[26] on extremum seeking under delays.

Here, we pursue a different objective under additive measurement uncertainty on the unknown objective function. We incorporate the effects of the uncertainty in overshoot terms in the upper bound on the norm of the estimation error. We introduce a state transformation, which allows us to use a new time-varying quadratic Lyapunov function. This function is reminiscent of Lyapunov functions that are generated by ‘strictification’ [27] (which transforms quadratic functions into strict Lyapunov functions, and so provides a constructive variant of the Matrosov approach [28]). We also apply a comparison lemma to a differential inequality having a square root of the state on its right side. We illustrate how our new methods can achieve much larger bounds on the parameter $\epsilon$ in the first dither period $\omega_1 = 2\pi/\epsilon$ (and so a usefully smaller $\omega_1$) while reducing the ultimate bounds on the estimation error, compared with averaging-based works such as [20].

Instead of averaging, a key ingredient that makes our analysis work is our new state transformation and Lyapunov function that enable us to cancel undesirable overshoot terms, to obtain an ultimate bound that is $O(\sqrt{\epsilon})$ in the special case where no measurement uncertainty is present. This note therefore provides a higher dimensional analog of the conference version [1] of this paper, which was confined to bounded single variable extremum seeking, meaning, the unknown objective function to be maximized or minimized was a function of one variable, instead of two variables as in this note. Moreover, whereas [1] placed the measurement uncertainty inside the unknown objective function as additive uncertainty on the input, here we cover important cases having measurable locally essentially bounded added measurement uncertainties on the objective function. This requires a new state transformation and a new time-varying Lyapunov function construction that were beyond the scope of all earlier works.

## II. EXTREMUM SEEKING PROBLEM AND THEOREM

Following [1] but allowing the objective function to instead be a function of two variables, we consider a two input real valued function $Q(\theta(t))$ that has the quadratic form

$$Q(\theta(t)) = Q^* + \frac{1}{2} \begin{pmatrix} \theta_1(t) - \theta_1^* \\ \theta_2(t) - \theta_2^* \end{pmatrix}^\top H \begin{pmatrix} \theta_1(t) - \theta_1^* \\ \theta_2(t) - \theta_2^* \end{pmatrix} + \delta(t) \quad (1)$$

Submitted March 31, 2024. Supported by US National Science Foundation Grant 2308282 and Israel Science Foundation Grant 446/24. A preliminary version [1] for one variable static maps was accepted for presentation at the 63rd IEEE Conference on Decision and Control; see Section I for a comparison of the conference version with this version.

F. Mazenc is with INRIA EPI DISCO, L2S-CNRS-CentraleSupélec, 3 rue Joliot Curie, 91192, Gif-sur-Yvette, France (e-mail: frederic.mazenc@l2s.centralesupelec.fr).

M. Malisoff is with Department of Mathematics, 303 Lockett Hall, 175 Field House Drive, Louisiana State University, Baton Rouge, LA 70803, USA (e-mail: malisoff@lsu.edu).

E. Fridman is with the School of Electrical Engineering, Tel Aviv University, Tel Aviv 69978, Israel (e-mail: emilia@tauex.tau.ac.il).
where $H = [h_{ij}] \in \mathbb{R}^{2 \times 2}$ is an unknown $2 \times 2$ matrix, the real constants $Q^*$, $\theta_1^*$, and $\theta_2^*$ are unknown, and the unknown measurable locally essentially bounded function $\delta$ represents measurement uncertainty. A key assumption throughout this note is that $\theta = [\theta_1, \theta_2]^\top$ is valued in $\mathbb{R}^2$, which is needed for the existence of a known ratio $\ell$ of the dither periods; see (6). Although $H$, $\theta_1^*$, $\theta_2^*$, and $\delta$ are unknown, we assume the following, where $|\cdot|_\infty$ is the essential supremum norm, and $M_{2+}$ is the set of all $2 \times 2$ real matrices whose eigenvalues both have positive real parts and $|\cdot|$ is the usual Euclidean norm, but analogs where $H$ is negative definite can be proven by replacing $Q^*$ by $-Q^*$ in what follows:

**Assumption 1.** There is a known compact set $\mathcal{H} \subseteq M_{2+}$ such that $H \in \mathcal{H}$. Also, there are known real values $\underline{\theta}_i$ and $\bar{\theta}_i$ for $i=1,2$ such that $\theta_i^* \in [\underline{\theta}_i, \bar{\theta}_i]$ for $i=1,2$. Also, there is a known constant $\bar{\delta} > 0$ such that $|\delta|_\infty \le \bar{\delta}$. $\square$

We then set

$$ \underline{\theta} = [\underline{\theta}_1, \underline{\theta}_2]^\top, \theta^* = [\theta_1^*, \theta_2^*]^\top, \text{ and } \bar{\theta} = [\bar{\theta}_1, \bar{\theta}_2]^\top, $$

so $\underline{\theta} \le \theta^* \le \bar{\theta}$ in the componentwise sense. By using the compactness of the set $\mathcal{H}$ in Assumption 1, we can construct positive constants $\underline{p}, \bar{p}$, and $q_0$ such that there exist a positive definite $2 \times 2$ matrix $P$ that satisfies

$$ -\frac{1}{2} P H - \frac{1}{2} H^\top P \le -q_0 P \quad (2) $$

and $pI \le P$ and $|P| \le \bar{p}$, where inequalities $A \le B$ for square real matrices of the same size mean that $B-A$ is nonnegative definite, $I$ is the identity matrix, and $|\cdot|$ is the matrix operator 2-norm; see Remark 1 for a way to find the required constants $\underline{p}, \bar{p}$, and $q_0$, using uniqueness of solution properties of Lyapunov equations and continuity properties of eigenvalues from [29], which cover the cases that we study in this note where the $H$'s (and therefore also the $P$'s satisfying (2) for some constant $q_0 > 0$) are unknown. We fix $\underline{p}, \bar{p}$, and $q_0$ satisfying the preceding requirements in what follows.

The objective of extremum seeking is to find a real-time estimate $\theta(t)$ of the extremum point $\theta^*$ (meaning, $\lim_{t \to +\infty} \theta(t) = \theta^*$) based on the disturbed measurements of the map

$$ y(t) = Q^* + \frac{1}{2}(\theta(t) - \theta^*)^\top H(\theta(t) - \theta^*) + \delta(t). \quad (3) $$

In terms of the estimation error vector

$$ \tilde{\theta}(t) = \begin{bmatrix} \tilde{\theta}_1 \\ \tilde{\theta}_2 \end{bmatrix}, \text{ where } \tilde{\theta}_i(t) = \theta_i(t) - \theta_i^* \text{ for } i=1,2, \quad (4) $$

and following [14], we therefore consider the two state gradient based bounded extremum seeking dynamics

$$ \dot{\tilde{\theta}}_i(t) = \sqrt{\alpha \omega_i} \cos(\omega_i t + k y(t)) \quad \text{for } i=1,2, \quad (5) $$

for constants $\alpha > 0$ and $k > 0$, where the known positive constants $\omega_i$ for $i=1,2$, $\ell > 1$, and $\epsilon > 0$ are such that

$$ \omega_2 = \ell \omega_1 \text{ and } \omega_1 = \frac{2\pi}{\epsilon}. \quad (6) $$

Throughout this work, all equalities and inequalities that include $\delta$ should be regarded to be holding for almost all times $t \ge 0$, in the Lebesgue measure sense. We will introduce conditions involving $\bar{\delta}$ and $\epsilon$ that allow us to obtain a suitable

ultimate bound on $\tilde{\theta}(t)$, i.e., a constant $B_U > 0$ such that $\lim \sup_{t \to +\infty} |\tilde{\theta}(t)| \le B_U$. To express our conditions and ultimate bound, it is convenient to introduce the following constants, in addition to the constants $\alpha, k, q_0, \underline{p}, \bar{p}, \epsilon$, and $\ell$ that we introduced above. We use the constants

$$ \bar{h} = k\alpha \sup_{H \in \mathcal{H}} |H|, \quad k^\# = k \left(1 + \bar{h} \left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}\right), \quad (7) $$

$$ \bar{M} = \frac{1}{2\pi} \left( \frac{1}{4} + \frac{\ell^{3/2}}{\ell^2 - 1} \right), \quad \bar{N} = \frac{1}{\pi} \left( \frac{1}{4} + \frac{\sqrt{\ell}}{2(\ell+1)} \right), \quad (8) $$

$$ \nu_1 = \sqrt{2}(1 + \sqrt{\ell})\bar{h}\sqrt{\frac{\ell+\ell}{\ell}\bar{p}}(1 + 2\overline{M}\epsilon\bar{h})\frac{\sqrt{\epsilon}}{\sqrt{\pi}} + 2\bar{p}(1 + 2\overline{M}\epsilon\bar{h})\sqrt{1+\ell}\sqrt{\frac{2\pi}{\epsilon}}k^\#\bar{\delta}, $$

$$ \nu_2 = 2\bar{p}\overline{N}\bar{h}^2\epsilon\sqrt{\ell+1}\sqrt{1+\frac{1}{\ell}} + 4\bar{p}\overline{M}\bar{h}^2(1+\sqrt{\ell})\epsilon + 2\bar{p}k\bar{h}\bar{\delta}(1+2\overline{M}\epsilon\bar{h})\sqrt{1+\ell}\sqrt{1+\frac{1}{\ell}}, \quad (9) $$

$$ \nu_3 = 2\bar{p}\overline{N}\bar{h}^2\sqrt{2\pi\epsilon(\ell+1)}, \quad \lambda_s = \frac{1}{2c} \left[ \frac{k\alpha q_0}{3} - \sqrt{\left(\frac{k\alpha q_0}{3}\right)^2 - \frac{3c\nu_1}{\sqrt{2p}}} \right], \quad c = \frac{27\nu_3}{32\sqrt{2p}\sqrt{3}}, \quad (10) $$

$$ \lambda_l = \frac{1}{2c} \left[ \frac{k\alpha q_0}{3} + \sqrt{\left(\frac{k\alpha q_0}{3}\right)^2 - \frac{3c\nu_1}{\sqrt{2p}}} \right], \quad \text{and} \quad \sigma_0 = |\bar{\theta} - \theta|, \quad (11) $$

where the fact that $\lambda_s$ and $\lambda_l$ are positive real numbers will follow from the following assumptions on our constants, which can be satisfied for any given values of the other constants that we defined above when $\epsilon > 0$ and $\bar{\delta} > 0$ are small enough:

**Assumption 2.** The four conditions

$$ \nu_2 \le \frac{1}{4}pk\alpha q_0, \quad (12) $$

$$ \epsilon \le \frac{p}{18\bar{p}\overline{M}\bar{h}}, \quad (13) $$

$$ \frac{3c\nu_1}{\sqrt{2p}} < \frac{(k\alpha q_0)^2}{9}, \quad (14) $$

and

$$ \frac{3}{2} \sqrt{\frac{p}{2}} \left( \frac{\sigma_0}{\sqrt{\alpha}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \right) < \lambda_l \quad (15) $$

are satisfied. $\square$

See Section V for examples showing how Assumption 2 allows larger $\epsilon$'s and smaller ultimate bounds on the estimation errors, as compared with previous results. The constants $\nu_i$ from (7)-(11) will play an essential role in defining coefficients arising in a decay estimate for a quadratic Lyapunov function in the proof of our theorem, because this decay estimate (in conjunction with a comparison lemma argument) will be key to proving our ultimate bound on the norm of the error variable. Also, $\lambda_s$, $\lambda_l$, and $c$ from (10)-(11) are essential for defining the rate of exponential convergence of the norm of the estimation error $\tilde{\theta}$ to the ultimate bound on the estimation error. Hence, our main result goes beyond finding an ultimate bound, by also computing rates of convergence of $|\tilde{\theta}(t)|$ towards the ultimate bound. To express the convergence rate and simplify the analysis, we also use the positive constant

$$ d_* = \lambda_l - \max \left\{ \lambda_s, \sqrt{\frac{9}{8}\bar{p}} \left( \frac{\sigma_0}{\sqrt{\alpha}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \right) \right\}, \quad (16) $$

where the positivity of $d_*$ follows from condition (15) from our Assumption 2. Using the preceding constants, we prove:
**Theorem 1.** Let Assumptions 1-2 hold. Then for each initial state $\theta(0) \in \mathbb{R}^2$ such that $|\tilde{\theta}(0)| \le \sigma_0$, the corresponding solution $\tilde{\theta}(t)$ of (5) is such that the condition

$$|\tilde{\theta}(t)| \le \frac{9\lambda_l \sqrt{\alpha \bar{p}}}{8d_* \sqrt{\bar{p}}} \left( \frac{\sigma_0}{\sqrt{\alpha}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \right) e^{-c(\lambda_l - \lambda_s)t} + \frac{3\lambda_s}{2} \sqrt{\frac{\alpha}{2\bar{p}}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\alpha \epsilon}{2\pi}} \quad (17)$$

is satisfied for all $t \ge 0$. □

**Remark 1.** Given any $\sigma_0 > 0$, Assumption 2 is satisfied for small enough $\epsilon$ and $\bar{\delta}$, so our results are semiglobal (because $\sigma_0$ only appears in (15), and because $\lambda_l = O(1/\sqrt{\epsilon})$). The required positive constants $\bar{p}$, $\bar{p}$, and $q_0$ can be found by first noting that for each $H \in \mathcal{H}$, the unique solution of a Lyapunov equation provides a positive definite $2 \times 2$ matrix $P(H)$ such that $-P(H)H - HP(H) = -I$, and such that $P(H)$ is a continuous function of $H$ on its domain $\mathcal{H}$ (e.g., by [29, Chapter 5]). Then we choose $\bar{p}$ to be a positive lower bound for all of the eigenvalues of the set of matrices $S = \{P(H) : H \in \mathcal{H}\}$, $\bar{p}$ to be an upper bound for the norms of all matrices in this set $S$, and $q_0 = 1/(2\lambda_{\max})$ where $\lambda_{\max}$ is the largest eigenvalue of all matrices in $S$. This produces positive values $\bar{p}$, $\bar{p}$, and $q_0$, by the continuity of norms and eigenvalues of matrices as functions of the matrix entries and the compactness assumption on $\mathcal{H}$ from Assumption 1.

For instance, using the Mathematica computer program, $P(H)$ can be expressed using the Riccati Solve command, and then $\bar{p}$, $\bar{p}$, and $q_0$ can be computed as the minimum of the Mathematica function Min[Eigenvalues[P[H]]], the maximum of Norm[P[H]], and the minimum of 0.5/Max[Eigenvalues[P[H]]] respectively over all $2 \times 2$ matrices $H$ in the compact set $\mathcal{H} \subseteq \mathcal{M}_{2+}$. This does not generate a $P$ that satisfies (2) for all $H \in \mathcal{H}$, but it does provide positive constants $\bar{p}$, $\bar{p}$, and $q_0$ that are independent of the particular choice of the positive definite matrix $P$ that satisfies (2) where $P$ depends on the unknown $H \in \mathcal{H}$. This suffices, because our requirements on the extremum seeking parameters and our ultimate bound use $\bar{p}$, $\bar{p}$, and $q_0$ instead of a formula for $P$, so we do not need a formula for $P$ that satisfies (2) for our unknown $H$. □

**Remark 2.** Theorem 1 is new, even when $\bar{\delta} = 0$. See Section V for discussions of how Theorem 1 can allow larger $\epsilon$'s and smaller ultimate bounds in this special case, compared with earlier results. Also, in this case, the ultimate bound

$$B_U = \frac{3\lambda_s}{2} \sqrt{\frac{\alpha}{2\bar{p}}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\alpha \epsilon}{2\pi}} \quad (18)$$

from (17) has order $O(\sqrt{\epsilon})$, which is desirable for small $\epsilon$'s, and our conditions from Assumption 2 have the convex property that if they hold for a value $\epsilon > 0$ and for a given set of values for the other parameters, then they also hold for all smaller positive $\epsilon$ values as well. On the other hand, see Section IV below, where we explain how the $\delta(t)$ in (1) can represent the effects of time-varying unknown delays. □

**Remark 3.** Using (9)-(11), we can rewrite the decay rate $r_* = c(\lambda_l - \lambda_s)$ in (17) as

$$r_* = \sqrt{\frac{(k\alpha q_0)^2}{9} - \frac{81\sqrt{2\nu_1\bar{p}\bar{N}\hbar^2}\sqrt{\pi\epsilon(\ell+1)}}{32\bar{p}^2}},$$

where $\bar{N}$ and $\nu_1$ are from (8)-(9). It follows that smaller values of $\epsilon$, or larger $\alpha$ or $k$, can provide faster convergence of the estimation error. In Section V, we illustrate how changing some of the parameter values can speed up the convergence.

### III. PROOF OF THEOREM 1

The proof has four parts. First, we introduce a state transformation that makes the dynamics for $\tilde{\theta}$ amenable to our time-varying quadratic Lyapunov function analysis. In the second part, we build our quadratic Lyapunov function $W$ that satisfies a differential inequality with a $\sqrt{\bar{W}}$ on its right side. In the third part, we apply a comparison lemma argument to construct a suitable time-varying upper bound for $W$. Then in the fourth part, we use appropriate upper and lower bounds for $W$ to obtain the final bound (17) on the estimation error $|\tilde{\theta}(t)|$.

**First Part: State Transformation.** For the rest of the proof, we fix a positive definite matrix $P$ satisfying (2) as well as our conditions $|P| \le \bar{p}$ and $P \ge \bar{p}I$ for the fixed choices of $q_0$, $\bar{p}$, and $\bar{p}$ from the previous section. Then (3)-(5) give

$$\tilde{\theta}_i(t) = \sqrt{\alpha \omega_i} \cos \left( \omega_i t + kQ^* + \frac{k}{2} \tilde{\theta}^\top(t) H \tilde{\theta}(t) + k\delta(t) \right) \quad (19)$$

for $i = 1, 2$.

To simplify our analysis, we use the change of variables

$$\Theta = \begin{bmatrix} \hat{\theta}_1 \\ \hat{\theta}_2 \end{bmatrix}, \quad \text{where } \hat{\theta}_1 = \frac{1}{\sqrt{\alpha}} \tilde{\theta}_1 \text{ and } \hat{\theta}_2 = \frac{1}{\sqrt{\alpha}} \tilde{\theta}_2. \quad (20)$$

Since (20) gives

$$k \tilde{\theta}^{\top}(t) H \tilde{\theta}(t) = k \sqrt{\alpha} \Theta(t)^{\top} H \sqrt{\alpha} \Theta(t) = k \alpha \Theta^{\top}(t) H \Theta(t),$$

we can use (19) to get

$$\begin{aligned} \dot{\tilde{\theta}}_i(t) &= \frac{1}{\sqrt{\alpha}} \tilde{\theta}_i \\ &= \sqrt{\omega_i} \cos \left( \omega_i t + kQ^* + \frac{k}{2} \tilde{\theta}^\top(t) H \tilde{\theta}(t) + k\delta(t) \right) \quad (21) \\ &= \sqrt{\omega_i} \cos \left( \omega_i t + \hat{Y}(t) + k\delta(t) \right) \end{aligned}$$

for $i = 1, 2$, where

$$\hat{Y}(t) = kQ^* + \frac{k\alpha}{2} \Theta(t)^{\top} H \Theta(t). \quad (22)$$

We next use the function

$$\mathcal{V}(t) = \begin{bmatrix} \sqrt{\omega_1} \cos(\omega_1 t + \hat{Y}(t)) + k\delta(t) \\ \sqrt{\omega_2} \cos(\omega_2 t + \hat{Y}(t)) + k\delta(t) \end{bmatrix}. \quad (23)$$

Then (21) can be rewritten as

$$\dot{\Theta}(t) = \mathcal{V}(t). \quad (24)$$

We also use the state transformation

$$\mathcal{W}(t) = \begin{bmatrix} \frac{1}{\sqrt{\omega_1}} \sin(\omega_1 t + \hat{Y}(t)) \\ \frac{1}{\sqrt{\omega_2}} \sin(\omega_2 t + \hat{Y}(t)) \end{bmatrix} \quad (25)$$

and $\mathcal{X}(t) = \Theta(t) - \mathcal{W}(t)$.

Then (24) gives

$$\begin{aligned} \dot{\mathcal{X}}(t) &= \mathcal{V}(t) - \begin{bmatrix} \frac{1}{\sqrt{\omega_1}} \cos(\omega_1 t + \hat{Y}(t)) (\omega_1 + \dot{\hat{Y}}(t)) \\ \frac{1}{\sqrt{\omega_2}} \cos(\omega_2 t + \hat{Y}(t)) (\omega_2 + \dot{\hat{Y}}(t)) \end{bmatrix} \\ &= - \begin{bmatrix} \frac{1}{\sqrt{\omega_1}} \cos(\omega_1 t + \hat{Y}(t)) \\ \frac{1}{\sqrt{\omega_2}} \cos(\omega_2 t + \hat{Y}(t)) \end{bmatrix} \dot{\hat{Y}}(t) + \mathbb{H}(t), \end{aligned} \quad (26)$$
where

$$
\mathbb{H}(t) = \begin{bmatrix} \sqrt{\omega_1} [\cos(\omega_1 t + \hat{Y}(t) + k\delta(t)) - \cos(\omega_1 t + \hat{Y}(t))] \\ \sqrt{\omega_2} [\cos(\omega_2 t + \hat{Y}(t) + k\delta(t)) - \cos(\omega_2 t + \hat{Y}(t))] \end{bmatrix} . \quad (27)
$$

Now, observe that with the choices

$$
\begin{align*}
C_{i\delta}(t) &= \sqrt{\omega_i} \cos(\omega_i t + \hat{Y}(t) + k\delta(t)) \\
\text{and } C_i(t) &= \sqrt{\omega_i} \cos(\omega_i t + \hat{Y}(t))
\end{align*}
$$

for $i = 1, 2$, we can use (22) and (24) to get

$$
\begin{align}
\hat{Y}(t) &= k\alpha V(t)^{\top} H \Theta(t) = k\alpha [C_{1\delta}(t) C_{2\delta}(t)] H \Theta(t) \tag{28} \\
&= k\alpha [C_1(t) C_2(t)] H \Theta(t) + k\alpha \mathbb{H}(t)^{\top} H \Theta(t). \nonumber
\end{align}
$$

Consequently, (26) gives

$$
\dot{X}(t) = -k\alpha \left[ \frac{\frac{1}{\sqrt{\omega_1}} \cos(\omega_1 t + \hat{Y}(t))}{\frac{1}{\sqrt{\omega_2}} \cos(\omega_2 t + \hat{Y}(t))} \right] [C_1(t) \ C_2(t)] H\Theta(t) + \mathbb{L}(t), \quad (29)
$$

where

$$
\mathbb{L}(t) = \mathbb{H}(t) - k\alpha \left[ \frac{\frac{1}{\sqrt{\omega_1}} \cos(\omega_1 t + \hat{Y}(t))}{\frac{1}{\sqrt{\omega_2}} \cos(\omega_2 t + \hat{Y}(t))} \right] \mathbb{H}(t)^{\top} \mathbb{H} \Theta(t). \quad (30)
$$

Then, with $\mathcal{G}$ defined by

$$
\begin{equation}
\begin{split}
\mathcal{G}(t) = & \begin{bmatrix} -\cos^2(\omega_1 t + \hat{Y}(t)) & -\sqrt{\frac{\ell}{\omega_1 \omega_2}} C_1(t) C_2(t) \\ -\frac{1}{\sqrt{\ell \omega_1 \omega_2}} C_1(t) C_2(t) & -\cos^2(\omega_2 t + \hat{Y}(t)) \end{bmatrix},
\end{split}
\tag{31}
\end{equation}
$$

we can use (29) to obtain $\dot{X}(t) = k\alpha\mathcal{G}(t)H\Theta(t) + \mathbb{L}(t)$ and so also

$$
\begin{align}
\dot{X}(t) &= k\alpha \mathcal{G}(t) H X(t) + k\alpha \mathcal{G}(t) H W(t) + \mathbb{L}(t) \nonumber \\
&= -\frac{k\alpha H}{2} X(t) + \left[ \mathcal{G}(t) + \frac{1}{2} \right] k\alpha H X(t) \tag{32} \\
&\phantom{= -\frac{k\alpha H}{2} X(t) + \left[ \mathcal{G}(t) + \frac{1}{2} \right] k\alpha H} + k\alpha \mathcal{G}(t) H W(t) + \mathbb{L}(t) \nonumber
\end{align}
$$

since (25) gives $\Theta(t) = X(t) + W(t)$.

Also, with the choice

$$
\mathcal{M}(t) = \frac{1}{\omega_1} \begin{bmatrix} \frac{1}{4} \sin(2\omega_1 t + 2\hat{Y}(t)) & \frac{\sqrt{\ell}}{2} S_1(t) \\ \frac{1}{2\sqrt{\ell}} \frac{\cos(\omega_1(1+\ell)t + 2\hat{Y}(t))}{\omega_1(1+\ell)} & \frac{1}{4\sqrt{\ell}} \sin(2\omega_2 t + 2\hat{Y}(t)) \end{bmatrix} \quad (33)
$$

where

$$
S_1(t) = \frac{\sin(\omega_1(1+\ell)t + 2\hat{Y}(t))}{1+\ell} + \frac{\sin(\omega_1(1-\ell)t)}{1-\ell},
$$

and with the choice

$$
\mathcal{N}(t) = \begin{bmatrix} \frac{1}{2\omega_1} \cos(2\omega_1 t + 2\hat{Y}(t)) & \sqrt{\ell} \frac{\cos(\omega_1(1+\ell)t + 2\hat{Y}(t))}{\omega_1(1+\ell)} \\ \frac{1}{\sqrt{\ell}} \frac{\cos(\omega_1(1+\ell)t + 2\hat{Y}(t))}{\omega_1(1+\ell)} & \frac{1}{2\ell\omega_1} \cos(2\omega_2 t + 2\hat{Y}(t)) \end{bmatrix}, \quad (34)
$$

we deduce $\dot{\mathcal{M}}(t) = -(\mathcal{G}(t) + I/2) + \mathcal{N}(t)\hat{Y}(t)$ (by the double angle formula and sum rule for cosine) and so also

$$
\mathcal{G}(t) + \frac{I}{2} = -\dot{\mathcal{M}}(t) + \mathcal{N}(t)\hat{Y}(t). \quad (35)
$$

By substituting (35) into (32), we obtain

$$
\begin{align*}
\dot{X}(t) &= -\frac{k\alpha H}{2} X(t) + k\alpha \left[ -\dot{\mathcal{M}}(t) + \mathcal{N}(t) \hat{Y}(t) \right] H X(t) \\
&\quad + k\alpha \mathcal{G}(t) H W(t) + \mathbb{L}(t) \\
&= -\frac{k\alpha H}{2} X(t) - k\alpha \dot{\mathcal{M}}(t) H X(t) + \mathbb{L}(t) \\
&\quad + k\alpha \mathcal{N}(t) \hat{Y}(t) H X(t) + k\alpha \mathcal{G}(t) H W(t).
\end{align*}
$$

From the first equality in (28), it follows that

$$
\begin{equation}
\begin{split}
\dot{X}(t) ={}& -\frac{k\alpha H}{2} X(t) - k\alpha \dot{\mathcal{M}}(t) H X(t) + \mathbb{L}(t) \\
& + (k\alpha)^2 \mathcal{N}(t) \mathcal{V}(t)^{\top} H \Theta(t) H X(t) \\
& + k\alpha \mathcal{G}(t) H W(t).
\end{split}
\tag{36}
\end{equation}
$$

Moreover, since for all constants $a > 0$ and $b > 0$, the matrix

$$
S = \begin{bmatrix} a & b \\ b & a \end{bmatrix}
$$

satisfies $|S| = a + b$ which follows because the symmetry of $S$ gives $|S| = \rho(S)$ where $\rho(S)$ is the spectral radius of $S$, and because $S[1 \ 1]^\top = (a + b)[1 \ 1]^\top$ (which allows us to use the Perron-Frobenius theorem to get $\rho(S) = a + b$), and since $\ell \ge 1$ and (6) hold, we obtain the bounds

$$
\begin{equation}
\begin{aligned}
|\mathcal{G}(t)| &\le 1 + \sqrt{\ell}, \quad |\mathcal{M}(t)| \le \bar{M}\epsilon, \\
|\mathcal{N}(t)| &\le \bar{N}\epsilon, \quad |\mathcal{V}(t)| \le \frac{\sqrt{2\pi(1+\ell)}}{\sqrt{\epsilon}}
\end{aligned}
\tag{37}
\end{equation}
$$

for the function $\mathcal{G}$ from (31), the $\mathcal{M}$ in (33), the $\mathcal{N}$ in (34), and the $\mathcal{V}$ in (23), where the constants $\bar{M}$ and $\bar{N}$ are from (8), and where we again used the relation $\omega_2 = \ell\omega_1$ from (6). The bounds (37) will play an essential role in our Lyapunov analysis in the next part of the proof.

Second Part: New Lyapunov Function Approach. Let us introduce the candidate Lyapunov function

$$
\begin{equation}
\begin{aligned}
W(t, X) &= V(X) + 2k\alpha X^\top P \mathcal{M}(t) H X, \\
\text{where } V(X) &= X^\top P X.
\end{aligned}
\tag{38}
\end{equation}
$$

Then our choice of $\bar{h}$ in (7) and the bound on $\mathcal{M}(t)$ in (37)
give

$$
\begin{align*}
W(t, X) &\le V(X) + 2k\alpha |X|^2 \bar{p} \bar{M} |H| \epsilon \\
&\le \left[ 1 + \frac{2\bar{p} \bar{M} h \epsilon}{\bar{p}} \right] V(X)
\end{align*}
$$

and

$$
W(t, X) \ge V(X) - 2|X|^2 \bar{p} \bar{M} h \epsilon \ge \left(1 - \frac{2\bar{p} \bar{M} h}{\bar{p}}\right) V(X)
$$

for all $X \in \mathbb{R}^2$ and all $t \ge 0$, where we used the inequality
$V(X) \ge \underline{p}|X|^2$. Our assumption (13) then implies that

$$
W(t, X) \le \frac{9}{8} V(X) \text{ and } V(X) \le \frac{9}{8} W(t, X) \quad (39)
$$

are satisfied for all $X$ and all $t \ge 0$, since $2\bar{p}\epsilon\bar{M}h/\underline{p} \le 1/9$.
This motivates our referring to $W$ as a candidate Lyapunov
function.

Also, by multiplying (2) through by $k\alpha$, it follows that along
the dynamics (36), we have

$$
\begin{align}
\dot{V}(t) &\le -k\alpha q_0 V(X(t)) + 2X(t)^\top P \left[ -k\alpha \dot{\mathcal{M}}(t) H X(t) \right. \nonumber \\
&\qquad \left. + \mathcal{N}(t) \mathcal{V}(t)^\top H \Theta(t) (k\alpha)^2 H X(t) \right. \nonumber \\
&\qquad \left. + k\alpha \mathcal{G}(t) H W(t) + \mathbb{L}(t) \right] \tag{40} \\
&= -k\alpha q_0 V(X(t)) - 2k\alpha X(t)^\top P \dot{\mathcal{M}}(t) H X(t) \nonumber \\
&\quad + 2X(t)^\top P \left[ (k\alpha)^2 \mathcal{N}(t) \mathcal{V}(t)^\top H \Theta(t) H X(t) \right. \nonumber \\
&\qquad \left. + k\alpha \mathcal{G}(t) H W(t) + \mathbb{L}(t) \right]. \nonumber
\end{align}
$$
Then, along the solutions of (36), we can use our formula
$\dot{X}(t) = k\alpha G(t)H\Theta(t) + \mathbb{L}(t)$ to obtain

$$
\begin{equation}
\begin{aligned}
\dot{W}(t) \le {}& -k\alpha q_0 V(X(t)) \\
& +\{2(k\alpha)^2 X(t)^\top P N(t) V(t)^\top H \Theta(t) H X(t)\} \\
& +2k\alpha X(t)^\top P G(t) H W(t) + 2X(t)^\top P \mathbb{L}(t) \\
& +2(k\alpha)^2 X(t)^\top P M(t) H G(t) H \Theta(t) \\
& +2(k\alpha)^2 [G(t) H \Theta(t)]^\top P M(t) H X(t) \\
& +2k\alpha \mathbb{L}(t)^\top P M(t) H X(t) \\
& +2k\alpha X(t)^\top P M(t) H \mathbb{L}(t).
\end{aligned}
\tag{41}
\end{equation}
$$

where we used the time derivative of the second term
$2k\alpha X^\top P M(t) H X$ in the formula for $W$ in (38) to cancel
the term containing $\mathcal{M}(t)$ in (40) and then collected terms.

By using our bound $k\alpha|H| \le \bar{h}$ from (7) to upper bound
the term in curly braces in (41) and then the bounds from (37),
we obtain

$$
\begin{equation}
\begin{aligned}
\dot{W}(t) \le {}& -k\alpha q_0 V(X(t)) \\
& + \left\{ 2\bar{p}\bar{N}\epsilon \frac{\sqrt{2\pi(\ell+1)}}{\sqrt{\epsilon}}\bar{h}^2 |\Theta(t)||X(t)|^2 \right\} \\
& + 2\bar{p}(1 + \sqrt{\ell})\bar{h}\frac{\sqrt{\epsilon}}{\sqrt{2\pi}}\sqrt{\frac{\ell+1}{\ell}}|X(t)| \\
& + \left\{ 4\bar{p}\bar{M}\epsilon\bar{h}^2 (1 + \sqrt{\ell})|\Theta(t)||X(t)| \right. \\
& \qquad \left. + 2\bar{p}(1 + 2\bar{M}\epsilon\bar{h})|\mathbb{L}(t)||X(t)| \right\}.
\end{aligned}
\tag{42}
\end{equation}
$$

Also, using our assumption from (6) that $\omega_2 = \ell\omega_1$, we get

$$
|\Theta(t)| \le |X(t)| + |W(t)| \le |X(t)| + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}}, \quad (43)
$$

by (25). Moreover, the function $\mathbb{L}$ in (30) satisfies

$$
\begin{align}
|\mathbb{L}(t)| &\le |\mathbb{H}(t)| + \sqrt{\frac{1}{\omega_1} + \frac{1}{\omega_2}} |\mathbb{H}(t)| \bar{h} |\Theta(t)| \tag{44} \\
&= |\mathbb{H}(t)| \left(1 + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi} \bar{h}} |\Theta(t)|\right). \nonumber
\end{align}
$$

If we now let $\ell^\#$ denote the second right side term in (43) (i.e.,
the radical), then it follows from using (43) to upper bound
the $|\Theta(t)|$ in (44) and then using our formula (27) for $\mathbb{H}$ that

$$
\begin{align}
|\mathbb{L}(t)| &\le |\mathbb{H}(t)| \left[1 + \ell^\# \bar{h} \left(|X(t)| + \ell^\#\right)\right] \nonumber \\
&\le \sqrt{\omega_1 + \omega_2 k \bar{\delta}} \left[1 + \ell^\# \bar{h} \left(|X(t)| + \ell^\#\right)\right] \tag{45} \\
&\le \sqrt{1 + \ell} \sqrt{\frac{2\pi}{\epsilon}} k \left[1 + \ell^\# \bar{h} \left(|X(t)| + \ell^\#\right)\right] \bar{\delta}, \nonumber
\end{align}
$$

since cos has global Lipschitz constant 1. Using (43) and (45)
to upper bound the terms in curly braces in (42), we obtain

$$
\begin{equation}
\begin{aligned}
\dot{W}(t) \le {}& -k\alpha q_0 V(X(t)) \\
& + \sqrt{2\bar{p}}(1 + \sqrt{\ell})\bar{h} \frac{\sqrt{\epsilon}}{\sqrt{\pi}} \sqrt{\frac{\ell+1}{\ell}} |X(t)| \\
& + 2\bar{p}\bar{N}\sqrt{\epsilon}\sqrt{2\pi(\ell+1)}\bar{h}^2 (|X(t)| + \ell^\#) |X(t)|^2 \\
& + 4\bar{p}\bar{M}\epsilon\bar{h}^2 (1 + \sqrt{\ell}) (|X(t)| + \ell^\#) |X(t)| \\
& + 2\bar{p} \left(1 + 2\bar{M}\epsilon\bar{h}\right) \sqrt{1 + \ell} \sqrt{\frac{2\pi}{\epsilon}} k \left[1 + \ell^\# \bar{h} (|X(t)| + \ell^\#)\right] \bar{\delta} |X(t)|.
\end{aligned}
\tag{46}
\end{equation}
$$

Collecting powers of $|X(t)|$ from the right side of (46) gives

$$
\begin{equation}
\begin{aligned}
\dot{W}(t) \le {}& -k\alpha q_0 V(X(t)) + \nu_1 |X(t)| \\
& \nu_2 |X(t)|^2 + \nu_3 |X(t)|^3
\end{aligned}
\tag{47}
\end{equation}
$$

where the positive constants $\nu_1$, $\nu_2$, and $\nu_3$ are from (9).

Since condition (12) from Assumption 2 and (39) give

$$
\begin{align*}
\nu_1 |X| &\le \nu_1 \sqrt{\frac{9}{8\bar{p}}} W(t, X), \\
\nu_2 |X|^2 &\le \frac{\nu_2 V(X)}{\bar{p}} \le \frac{k\alpha q_0 V(X)}{4}, \text{ and} \\
\nu_3 |X|^3 &\le \nu_3 \left(\frac{9}{8\bar{p}} W(t, X)\right)^{3/2},
\end{align*}
$$

we can upper bound the last three right side terms in (47),
then combine terms in the result, and then use the lower bound
$V(X) \ge 8W(t, X)/9$ from (39) to upper bound the resulting
combined term $-3k\alpha q_0 V(X(t))/4$. This produces the result

$$
\begin{equation}
\begin{aligned}
\dot{W}(t) \le {}& -\frac{2k\alpha q_0}{3} W(t, X(t)) + \nu_1 \sqrt{\frac{9}{8\bar{p}}} W(t, X(t)) \\
& \nu_3 \left(\frac{9}{8\bar{p}} W(t, X(t))\right)^{3/2}.
\end{aligned}
\tag{48}
\end{equation}
$$

Third Part: Comparison Argument. Inspired by the use of
comparison principles in trajectory based approaches (e.g.,
[30, Lemma 1]), we introduce the comparison system

$$
\dot{\xi}(t) = -\frac{2k\alpha q_0}{3}\xi(t) + \nu_1 \sqrt{\frac{9}{8\bar{p}}} \xi(t) + \nu_3 \left(\frac{9}{8\bar{p}}\right)^{\frac{3}{2}} \xi(t) \quad (49)
$$

and consider its positive valued solutions. Then the function

$$
\lambda(t) = \sqrt{\xi(t)} \tag{50}
$$

satisfies

$$
\dot{\lambda}(t) = -\frac{k\alpha q_0}{3}\lambda(t) + \frac{3\nu_1}{4\sqrt{2\bar{p}}} + \frac{27\nu_3}{32\sqrt{2\bar{p}}}\lambda^2(t). \quad (51)
$$

By factoring the right hand side of (51) as a polynomial in
$\lambda(t)$ and recalling our condition (14) from Assumption 2, we
can use our formulas for $c$, $\lambda_s$, and $\lambda_l$ from (10)-(11) to obtain

$$
\dot{\lambda}(t) = c(\lambda(t) - \lambda_s)(\lambda(t) - \lambda_l). \quad (52)
$$

On the other hand, our assumption in Theorem 1 that $|\tilde{\theta}(0)| \le \sigma_0$ and our change of variables (20) give $|\Theta(0)| \le \sigma_0/\sqrt{\alpha}$. Hence, (6) and our choice of $X$ in (25) give

$$
|X(0)| \le \frac{\sigma_0}{\sqrt{\alpha}} + \sqrt{\frac{1}{\omega_1} + \frac{1}{\omega_2}} = \frac{\sigma_0}{\sqrt{\alpha}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}},
$$

so our upper bound on $W$ in (39) gives

$$
W(0, X(0)) \le \frac{9}{8\bar{p}} \left( \frac{\sigma_0}{\sqrt{\alpha}} + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \right)^2 < \lambda_l^2, \quad (53)
$$

where (53) followed from (15).

We can now use (53) to apply the comparison principle to
(48) and (49), by choosing the positive initial state

$$
\xi(0) = \max \{ \lambda_s^2, W(0, X(0)) \} \quad (54)
$$

for (49), as follows. By condition (54), we have $\xi(0) \ge W(0, X(0))$. Then the comparison principle gives

$$
W(t, X(t)) \le \xi(t) = \lambda^2(t). \quad (55)
$$

We next use (55) to obtain the final bound (17).

Fourth Part: Using Bounds on W to Prove (17). Since we
can apply (53)-(54) and the relation (50) between $\xi$ and $\lambda$
to check that $\lambda(0) \in [\lambda_s, \lambda_l)$, we can apply a separation of
variables argument to (52) to obtain

$$
\lambda(t) = \frac{\lambda_s(\lambda_l - \lambda(0)) + \lambda_l(\lambda(0) - \lambda_s)e^{c(\lambda_s - \lambda_l)t}}{\lambda_l - \lambda(0) + (\lambda(0) - \lambda_s)e^{c(\lambda_s - \lambda_l)t}}. \quad (56)
$$
Using (20), the formulas for $X$ and $V$ in (25) and (38), the second inequality in (39), the upper bound on $W(t, X(t))$ from (55), the relation $\omega_2 = \ell\omega_1$, the subadditivity of the square root, and (56), we obtain

$$
\begin{align}
|\tilde{\theta}(t)| &\le \sqrt{\alpha}|\Theta(t)| \le \sqrt{\alpha}\left(|X(t)| + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}}\right) \nonumber \\
&\le \sqrt{\frac{\alpha}{\ell}V(X(t))} + \sqrt{\alpha}\sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \nonumber \\
&\le \sqrt{\frac{9\alpha}{8\ell}}W(t, X(t)) + \sqrt{\alpha}\sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \tag{57} \\
&\le \frac{3}{2}\sqrt{\frac{\alpha}{2\ell}}\lambda_s + \frac{3}{2}\lambda_l\sqrt{\frac{\alpha}{2\ell}\frac{\lambda(0)-\lambda_s}{\lambda_l-\lambda(0)}}e^{c(\lambda_s-\lambda_l)t} \nonumber \\
&\quad + \sqrt{\alpha}\sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}}, \nonumber
\end{align}
$$

where $c, \lambda_s$, and $\lambda_l$ are in (10)-(11). Also, we can use (20), (25), the first inequalities in (39) and (53), and (54) to get

$$
\begin{align}
\frac{\lambda(0)-\lambda_s}{\lambda_l-\lambda(0)} &\le \frac{1}{d_*} \sqrt{W(0, X(0))} \le \frac{3}{2d_*} \sqrt{\frac{\ell}{2}} |X(0)| \nonumber \\
&\le \frac{3}{2d_*} \sqrt{\frac{\ell}{2}} \left( |\Theta(0)| + \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \right) \tag{58} \\
&\le \frac{3}{2d_*} \sqrt{\frac{\ell}{2\alpha}} |\tilde{\theta}(0)| + \frac{3}{2d_*} \sqrt{\frac{\ell}{2}} \sqrt{\left(1 + \frac{1}{\ell}\right) \frac{\epsilon}{2\pi}} \nonumber
\end{align}
$$

where $d_*$ is defined by (16), and where the first inequality in (58) follows by separately considering the two cases $\lambda(0) = \sqrt{W(0, X(0))}$ and $\lambda(0) = \lambda_s$. Hence, we can combine (57)-(58) to choose the coefficient of $e^{-r*t}$ with the convergence rate $r_* = c(\lambda_l - \lambda_s)$ to obtain the final bound (17).

#### IV. USING $\delta(t)$ TO INCORPORATE MEASUREMENT DELAYS

We can use the $\delta$ in (1) to represent the effects of uncertain measurement delays. To see how, note that if the available measurement is

$$y(t) = Q^* + \frac{1}{2} \tilde{\theta}^\top (t - \tau(t)) H \tilde{\theta}(t - \tau(t))$$

where the piecewise continuous nonnegative valued bounded function $\tau$ represents an input delay, and if we set $B_* = c_* + B_U$ where $c_* > 0$ is the coefficient of $e^{-r*t}$ in (17) with the exponential convergence rate $r_* = c(\lambda_l - \lambda_s)$ and $B_U > 0$ is from (18), then we can use (5) to express the time-delayed measurement in the form (3) for a $\delta$ satisfying

$$
\begin{align*}
|\delta(t)| &= \frac{1}{2} |(\tilde{\theta}(t - \tau(t)))^\top H \tilde{\theta}(t - \tau(t)) - \tilde{\theta}^\top(t) H \tilde{\theta}(t)| \\
&\le \frac{1}{2} |(\tilde{\theta}(t - \tau(t)))^\top H (\tilde{\theta}(t - \tau(t)) - \tilde{\theta}(t))| \\
&\quad + \frac{1}{2} |(\tilde{\theta}(t - \tau(t)) - \tilde{\theta}(t))^\top H \tilde{\theta}(t)| \\
&\le B_* |H| \sqrt{\alpha(\omega_1 + \omega_2)} |\tau|_\infty
\end{align*}
$$

for all $t \in [0, \bar{T})$, where $\bar{T} > 0$ is the supremum of all $t \ge 0$ such that $|\tilde{\theta}(s)| \le B_*$ for all $s \in [0, t)$, where we assume that the initial function for $\tilde{\theta}$ is constant on $[-|\tau|_\infty, 0]$. By the definition of the supremum and the continuity of $\tilde{\theta}$, it follows that if $\bar{T} < +\infty$, then $|\tilde{\theta}(\bar{T})| = B_*$. On the other hand, we can apply the same argument as in the proof of Theorem 1 except only on the time interval $[0, \bar{T})$ with the choice

$$\bar{\delta} = B_* |H| \sqrt{\alpha(\omega_1 + \omega_2)} |\tau|_\infty$$

to obtain (17) for all $t \in [0, \bar{T})$, provided Assumptions 1-2 are satisfied with this choice of $\bar{\delta}$, which holds for sufficiently small $|\tau|_\infty$ and $\epsilon > 0$. Then $\bar{T} > 0$, because the bound $\sigma_0$ on the norm of the initial state satisfies $\sigma_0 < B_*$. Hence, if $\bar{T} < +\infty$,

then (17) would give $|\tilde{\theta}(\bar{T})| < B_*$, since the exponential convergence rate in (17) is positive. This contradiction implies that $\bar{T} = +\infty$, so the conclusion of Theorem 1 remains true under the preceding assumptions when the delay $\tau$ is present.

The preceding analysis places a bound on the allowable suprema $|\tau|_\infty$ of the delays. We next show why one cannot prove an attractivity result for a neighborhood of the origin without a restriction on $|\tau|_\infty$ by instead choosing $\epsilon$ or $\alpha$ sufficiently small. This will demonstrate that the cases of unknown time-varying and known constant delays are very different; see, e.g., [31], which ensured attractivity of a neighborhood of the origin for arbitrarily large constant delays of the form $\tau = L\epsilon$ for any integer $L > 1$. This sheds light on the need for a smallness condition on $|\tau|_\infty$ that we specified above.

To show the impact of time-varying delays, consider a delay $\tau(t) = t - \eta \text{Floor}(t/\eta)$, of sawtooth type, where Floor($x$) is the largest integer in $[0, x]$ for each $x \ge 0$, and where the constant $\eta > 0$ is such that

$$\epsilon = \frac{\eta}{r} \qquad (59)$$

where $r$ is any integer larger than 2. We use the sequence

$$t_j = j\eta \qquad (60)$$

indexed by integers $j \ge 0$ and the corresponding first component of the extremum seeking error dynamics

$$\dot{\tilde{\theta}}_1(t) = \sqrt{\alpha\omega_1} \cos \left( \omega_1 t + kQ^* + \frac{k}{2} \tilde{\theta}^\top(t_j) H \tilde{\theta}(t_j) \right)$$

for all $t \in [t_j, t_{j+1})$ and $j \ge 0$. Then, for each $j \ge 0$, we get

$$
\begin{align}
\tilde{\theta}_1(t_{j+1}) - \tilde{\theta}_1(t_j) &= \sqrt{\alpha\omega_1} \int_{t_j}^{t_{j+1}} \cos(\omega_1 m + Q^\#(j)) \, dm \nonumber \\
&= \sqrt{\alpha\omega_1} \int_{t_j + \frac{Q^\#(j)}{\omega_1}}^{t_{j+1} + \frac{Q^\#(j)}{\omega_1}} \cos(\omega_1 m) \, dm \tag{61} \\
&= \sqrt{\frac{\alpha}{\omega_1}} \left[ \sin(\omega_1 \eta + \omega_1 t_j + Q^\#(j)) - \sin(\omega_1 t_j + Q^\#(j)) \right], \nonumber
\end{align}
$$

where $Q^\#(j) = kQ^* + \frac{k}{2}\tilde{\theta}^\top(t_j)H\tilde{\theta}(t_j)$ and where the last equality in (61) used (60) to get $t_{j+1} = t_j + \eta$. By (59), we get

$$\omega_1 \eta = \frac{2\pi\eta}{\epsilon} = 2\pi r,$$

so (61) implies that $\tilde{\theta}_1(t_{j+1}) - \tilde{\theta}_1(t_j) = 0$, which gives $\tilde{\theta}_1(t_j) = \tilde{\theta}_1(0)$ for each integer $j > 0$. Since $r$ in (59) is arbitrarily large, we conclude that for any value $\eta > 0$, there are arbitrarily small values of $\epsilon$ for which the system cannot admit an attractive neighborhood of the origin.

#### V. ILLUSTRATIONS

In many significant cases, we found that Theorem 1 from Section II above led to much larger bounds on the allowable values of $\epsilon$, and to decreases in the ultimate bounds on $\tilde{\theta}(t)$, as compared with the previous state of the art results that were based on a time delay approach to averaging, including cases covered by [21] where the measurement uncertainty $\delta$ was zero. Therefore, we can obtain closer approximations of the unknown parameter vector. Due to the importance of approximating parameter vectors, this confirms the usefulness of the work, from a practical point of view.
To illustrate this, we first consider an autonomous vehicle in an environment without GPS orientation from [14] and [21, Section 4.2]. The goal is to reach the location of the stationary minimum of $Q$ with the diagonal matrix $H = \text{diag}\{2, 2\}$ and $Q^* = 0$, which will be achieved by using bounded extremum seeking. Following [21, Section 4.2], we use $k = 11$, $\alpha = 0.0001$, $\ell = 2$, and $\delta = 0$, which allowed us to satisfy the assumptions of Theorem 1 above for small enough $\epsilon$ values. Following [21], we separately considered the cases

$$ \sigma_0 = 0.001 \qquad (62) $$

and

$$ \sigma_0 = \sqrt{2}. \qquad (63) $$

We summarize our comparison in Table 1 below, where the upper bounds on $\epsilon$ are denoted by $\epsilon^*$ and the ultimate upper bounds UB are for $|\hat{\theta}(t)|$, using the formulas from Theorem 1 above and the methods from [21], where ES (resp., BES) indicate the classical (resp., bounded) extremum seeking method from [21]. Since our method produced much bigger values of $\epsilon^*$ and smaller extremum seeking estimation errors, and since our ultimate bound has the additional desirable property of being $O(\sqrt{\epsilon})$ when $\bar{\delta} = 0$ (which is consistent with [20]), our method from Theorem 1 above can offer significant advantages, as compared with earlier extremum seeking methods. In Fig. 1, we show Mathematica simulations of the estimation

<table><caption>TABLE I COMPARISON OF $\epsilon$ BOUNDS AND ULTIMATE BOUNDS WITH $\delta = 0$</caption><thead><tr><th>Extremum Seeking Method</th><th>ε*</th><th>UB</th></tr></thead><tbody><tr><td>BES from [21] with (62)</td><td>0.0004</td><td>1.41</td></tr><tr><td>ES from [21] with (62)</td><td>3.33</td><td>1.26</td></tr><tr><td>Theorem 1 above with (62)</td><td>8.84</td><td>0.170205</td></tr><tr><td>BES from [21] with (63)</td><td>0.0001</td><td>2.05</td></tr><tr><td>ES from [21] with (63)</td><td>0.36</td><td>0.94</td></tr><tr><td>Theorem 1 above with (63)</td><td>1.48</td><td>0.0563498</td></tr></tbody></table>

error dynamics (5) using the above values, for the case (63) from Table I. We choose the initial state $\hat{\theta}(0) = [1, -1]^T$, with $\epsilon = 1.48$ and therefore also $\omega_1 = \frac{2\pi}{\epsilon} = 4.245$ and $\omega_2 = 2\omega_1 = 8.4908$.

We obtained faster convergence by increasing $\alpha$. For instance, in Fig. 2, we show Mathematica simulations for the error dynamics (5) obtained using the same parameter values that we used to generate the simulation in Fig. 1, except we increased $k$ from $k = 11$ to $k = 11.1$, we increased $\alpha$ from $\alpha = 0.0001$ to $\alpha = 0.00022$, and we reduced the $\epsilon$ value from $\epsilon = 1.48$ to $\epsilon = 0.65$. Compared with Fig. 1, the simulation in Fig. 2 achieved close approximation of the unknown minimum approximately 2 times faster than Fig. 1, illustrating the influence of the parameter values on the convergence rate. Moreover, our Assumptions 1-2 were satisfied in both cases, which illustrates the applicability of Theorem 1 above for different parameter values. With the choices of the parameters that we used in Fig. 2, our Theorem 1 provided the ultimate bound $B_U = 0.0553486$.

Another notable feature of our approach is that we allow $\ell$ to be irrational. We illustrate this in Fig. 3, where we used the same parameter values that we used in our simulation from Fig. 1, except we replaced $\ell = 2$ by $\ell = \sqrt{2}$. With these choices, the assumptions of Theorem 1 were again satisfied.

As is seen in Fig. 3, replacing $\ell$ by an irrational number produced slightly faster convergence of the estimation error to zero. Moreover, for the set of parameter values that we used in Fig. 3, the ultimate bound decreased to $B_U = 0.0548245$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Solutions of (5) with initial state $\tilde{\sigma}(0) = [1, -1]^T$ for (63) showing $\tilde{\theta}_1(t)$ (dashed red) and $\tilde{\theta}_2(t)$ (solid blue) with $k = 11$, $\ell = 2$, and $\epsilon = 1.48$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Solutions of (5) with initial state $\tilde{\sigma}(0) = [1, -1]^T$ for (63) showing $\tilde{\theta}_1(t)$ (dashed red) and $\tilde{\theta}_2(t)$ (solid blue) with $k = 11.1$, $\ell = 2$, and $\epsilon = 0.65$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Solutions of (5) with initial state $\tilde{\sigma}(0) = [1, -1]^T$ for (63) showing $\tilde{\theta}_1(t)$ (dashed red) and $\tilde{\theta}_2(t)$ (solid blue) with $k = 11$, $\ell = \sqrt{2}$, and $\epsilon = 1.48$</figcaption></figure>

With the same parameter values as we used in Table 1 except with $\delta = 0$ and $\epsilon = 0.0004$, Theorem 1 above again applied, and it gave the ultimate bound 0.0009 using (62) and with (63). The convergence we showed above behaves favorably, as compared to what we obtained when we instead used the classical extremum seeking dynamics

$$ \begin{aligned} \dot{\hat{\theta}}_1(t) &= \frac{2k_1}{a_1} \sin(\omega_1 t) J(t) \\ \dot{\hat{\theta}}_2(t) &= \frac{2k_2}{a_2} \sin(2\omega_1 t) J(t) \end{aligned} \qquad (64) $$

from [21] and [25, Section II.C] where $J(t) = (\hat{\theta}_1(t) + a_1 \sin(\omega_1 t))^2 + (\hat{\theta}_2(t) + a_2 \sin(2\omega_1 t))^2$ with the maximum value $\epsilon = 0.36$ from Table I and with the parameter values (63), $k_1 = k_2 = -0.01$, $a_1 = a_2 = 0.2$, and $\ell = 2$ that were used for (64) in [25]. In classical extremum seeking, one starts with $\theta_i = \hat{\theta}_i + a_i \sin(\omega_i t)$ for $i = 1, 2$, since the objective of extremum seeking is to find $\theta(t)$. We show
results from Mathematica simulations of the solutions from (64) in Fig. 4 using the preceding parameter values, where we plot $\theta_i(t) = \hat{\theta}_i(t) + a_i \sin(\omega_i t)$ for $i = 1, 2$ to allow a fair comparison with the plots of the $\hat{\theta}_i$'s in the previous figures, because $\theta^* = 0$ in this case (so $\tilde{\theta} = \theta$).

<figure><img src="image_5.png" /><figcaption>Fig. 4. Solutions using (64) with $\theta_1(t)$ (dashed red), $\theta_2(t)$ (solid blue), $k_1 = k_2 = -0.01$, $a_1 = a_2 = 0.2$, $\omega_1 = 2\pi/\epsilon$, $\omega_2 = 4\pi/\epsilon$, and $\epsilon = 0.36$</figcaption></figure>

Finally, when we instead used the bound $\bar{\delta} = 0.01$ on $|\delta|$, and using $\epsilon = 40$, $\alpha = 0.0004$, $H = \text{diag}\{2, 2\}$, $k = 0.02$, $\ell = \sqrt{2}$, and $\sigma_0 = \sqrt{2}$, Theorem 1 above provided the ultimate bound of $1.08658$ on $|\hat{\theta}(t)|$. In this final case, we reduced $k$ in order to satisfy our Assumption 2, since the assumption was not satisfied with our previous choice of $k = 11$. Hence, we cover measurement uncertainties that were beyond the scope of earlier work, and our Theorem 1 remains effective using the smaller $\epsilon$'s from [21] as well. In practice, the selection of the parameters $H$ and $k$ in the extremum seeking should be based on the specific application; see, e.g., [21] for a discussion of parameter choices in the special case of vehicle control. For the Mathematica code used in this section, see https://github.com/MichaelMalisoff/TAC2DES.git.

## VI. CONCLUSION

We advanced the state of the art for stability analysis of bounded gradient based extremum seeking for static quadratic maps, by using a new state transformation and a new time-varying Lyapunov function approach that allow us to achieve significantly smaller dither periods and smaller ultimate bounds on the estimation errors, compared with previous extremum seeking results that instead used averaging. Our approach allows measurement uncertainties that are conducive to modeling the effects of unknown time-varying measurement delays. We aim to provide analogs for Newton-based extremum seeking with measurement uncertainty, and to study applications to source seeking in aerial applications.

## REFERENCES

[1] F. Mazenc, M. Malisoff, and E. Fridman, "Bounded extremum seeking for single-variable static map using new change of coordinates," in *Proc. 63rd IEEE Conference on Decision and Control*, 2024, to appear.

[2] M. Leblanc, "Sur l'electrification des chemins de fer au moyen de courants alternatifs de fréquence élevee," *Revue Generale de l'Electricite*, vol. 12, no. 8, pp. 275-277, 1922.

[3] M. Krstic and H.-H. Wang, "Stability of extremum seeking feedback for general nonlinear dynamic systems," *Automatica*, vol. 36, no. 4, pp. 595-601, 2000.

[4] A. Scheinker, "100 years of extremum seeking: A survey," *Automatica*, vol. 161, no. 111481, 2024.

[5] A. Scheinker and M. Krstic, *Model-Free Stabilization by Extremum Seeking*. New York, NY: Springer, 2017.

[6] T. Oliveira and M. Krstic, *Extremum Seeking Through Delays and PDEs*. Philadelphia, PA: SIAM, 2022.

[7] Cao, Z. and Durr, H.-B. and Ebenbauer, C. and Allgower, F. and Gao, F., "Iterative learning and extremum seeking for repetitive time-varying mappings," *IEEE Transactions on Automatic Control*, vol. 62, no. 7, pp. 3339-3353, 2016.

[8] Guay, M. and Dochain, D., "A proportional-integral extremum-seeking controller design technique," *Automatica*, vol. 77, pp. 61-67, 2017.

[9] A. Ghaffari, M. Krstic, and D. Nesic, "Multivariable Newton-based extremum seeking," *Automatica*, vol. 48, no. 8, pp. 1759-1767, 2012.

[10] Guay, M., "A time-varying extremum-seeking control approach for discrete-time systems," *Journal of Process Control*, vol. 24, no. 3, pp. 98-112, 2014.

[11] Guay, M. and Dochain, D., "A time-varying extremum-seeking control approach," *Automatica*, vol. 51, pp. 356-363, 2015.

[12] Hazeleger, L. and Haring, M. and van de Wouw, N., "Extremum-seeking control for optimization of time-varying steady-state responses of nonlinear systems," *Automatica*, vol. 119, no. 109068, 2020.

[13] H.-B. Durr, M. Stankovic, C. Ebenbauer, and K. Johansson, "Lie bracket approximation of extremum seeking systems," *Automatica*, vol. 49, no. 6, pp. 1538-1552, 2013.

[14] A. Scheinker and M. Krstic, "Extremum seeking with bounded update rates," *Systems & Control Letters*, vol. 63, pp. 25-31, 2014.

[15] A. Scheinker and D. Scheinker, "Bounded extremum seeking with discontinuous dithers," *Automatica*, vol. 69, pp. 250-257, 2016.

[16] Y. Tan, D. Nesic, and I. Mareels, "On non-local stability properties of extremum seeking control," *Automatica*, vol. 42, no. 6, pp. 889-903, 2006.

[17] C. Labar, C. Ebenbauer, and L. Marconi, "ISS-like properties in lie-bracket approximations and application to extremum seeking," *Automatica*, vol. 136, no. 110041, 2022.

[18] Scheinker, A. and Scheinker, D., "Extremum seeking for optimal control problems with unknown time-varying systems and unknown objective functions," *International Journal of Adaptive Control and Signal Processing*, vol. 35, no. 7, pp. 1143-1161, 2021.

[19] R. Suttner and S. Dashkovskiy, "Robustness and averaging properties of a large-amplitude, high-frequency extremum seeking control scheme," *Automatica*, vol. 136, no. 110020, 2022.

[20] J. Zhang and E. Fridman, "Stabilization under unknown control directions and disturbed measurements via a time-delay approach to extremum seeking," *International Journal of Robust and Nonlinear Control*, vol. 34, no. 10, pp. 6513-6531, 2024.

[21] Y. Zhu and E. Fridman, "Extremum seeking via a time-delay approach to averaging," *Automatica*, vol. 135, no. 109965, 2022.

[22] J. Zhang and E. Fridman, "Lie-brackets-based averaging of affine systems via a time-delay approach," *Automatica*, vol. 152, no. 110971, 2023.

[23] T. Oliveira, M. Krstic, and D. Tsubakino, "Extremum seeking for static maps with delays," *IEEE Trans. Automatica Control*, vol. 62, no. 4, pp. 1911-1926, 2017.

[24] M. Malisoff and M. Krstic, "Multivariable extremum seeking with distinct delays using a one-stage sequential predictor," *Automatica*, vol. 129, no. 109462, 2021.

[25] X. Yang and E. Fridman, "A robust time-delay approach to continuous-time extremum seeking for multi-variable static map," in *Proc. 62nd IEEE Conference on Decision and Control*, 2023, pp. 6768-6773.

[26] ——, "A time-delay approach to extremum seeking with large measurement delays," *IFAC-PapersOnLine*, vol. 56, no. 2, pp. 168-173, 2023.

[27] F. Mazenc, "Strict Lyapunov functions for time-varying systems," *Automatica*, vol. 39, no. 2, pp. 349-353, 2003.

[28] M. Malisoff and M. Krstic, "Stabilization and robustness analysis for a chain of exponential integrators using strict Lyapunov functions," *Automatica*, vol. 68, pp. 184-193, 2016.

[29] E. Sontag, *Mathematical Control Theory: Deterministic Finite Dimensional Systems, Second Edition*. New York, NY: Springer Science & Business Media, 1998.

[30] F. Mazenc, M. Malisoff, and S.-I. Niculescu, "Stability and control design for time-varying systems with time-varying delays using a trajectory-based approach," *SIAM Journal on Control and Optimization*, vol. 55, no. 1, pp. 533-556, 2017.

[31] X. Yang, E. Fridman, and B. Zhao, "Bounded extremum seeking for single-variable static map with large measurement delay via time-delay approach to averaging," in *Proc. 62nd IEEE Conference on Decision and Control*, 2023, pp. 3664-3669.