# Repetitive Control Tools for an Original Approach to Convex Optimization Problems under Affine Periodic Constraints

Daniele Astolfi, Cristiano Verrelli

► To cite this version:

Daniele Astolfi, Cristiano Verrelli. Repetitive Control Tools for an Original Approach to Convex Optimization Problems under Affine Periodic Constraints. Systems and Control Letters, 2025, 201, pp.106095. 10.1016/j.sysconle.2025.106095. hal-05042219

HAL Id: hal-05042219

https://hal.science/hal-05042219v1

Submitted on 22 Apr 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
# Repetitive Control Tools for an Original Approach to Convex Optimization Problems under Affine Periodic Constraints

Daniele Astolfi<sup>a</sup>, Cristiano M. Verrelli<sup>b</sup>

## Abstract

This paper provides a solution to the online convex optimization problem under a class of affine constraints, periodic with a known period. Functions whose minimizer vector exhibits a constant component within the kernel space of the constraint horizontal matrix are considered. By resorting to the latest developments in the repetitive control (RC) theory, two algorithms are originally presented: the first one resorting to the point-wise use of the delay as a universal periodic signal generator, the second one relying on the PDE (Partial Differential Equation) transport-equation-based theory. Both of them naturally extend the standard primal-dual algorithm acting in the constant constraint scenario, while guaranteeing global asymptotic convergence properties. Indeed, the two main different RC approaches in the literature are applied to the same optimization problem, while drawing original conclusions under the adoption of a common view. The derivation of an internal-model-based finite-dimensional (spectral) approximation for the latter introduces a further connection with the renowned adaptive learning control.

## 1. Introduction

This paper focuses on convex optimization problems with time-varying equality constraints that are available at run-time. These problems arise when an optimization algorithm has to continuously adjust the optimal operating point based on constraints that evolve over time and are known to the user at run-time. A key example is given by the power grid management [18], where distributed energy resources have to be optimized while considering voltage constraints and economic objectives in a decentralized framework. Unlike traditional centralized approaches, modern methods enable partitioned systems to autonomously coordinate, leading to time-varying equality constraints. Similar challenges arise: i) in wireless communications, where beamforming has to adapt to fluctuating channels, and ii) in industrial automation, where scheduling and resource allocation follow constraints repeated over time. These examples highlight the need for efficient optimization tools in dynamic environments. The reader is referred to [28, 29, 41, 25, 40] and references therein for further related examples of applications.

It is well-acknowledged that tools from control theory can be used to analyze optimization algorithms that ensure stability and convergence guarantees. Optimization algorithms are typically dynamical systems directly operating in the discrete-time domain or obtained from a time-discretization of a primary algorithm originally designed in the continuous-time domain. A typical continuous-time algorithm that has been extensively

studied in light of different approaches borrowed from control theory is the primal-dual algorithm or Arrow-Hurwicz-Uzawa flow [1], which is used to solve convex optimization problems in the presence of constraints. In particular, such an algorithm is obtained by resorting to the theory of Lagrange multipliers (or, more generically, from the Karush-Kuhn-Tucker conditions in the presence of inequality constraints): it relies on the gradient of the function be optimized and the use of a suitable integral action (see Section 3 of the present paper for further details). Nevertheless, several additional tools (e.g., Lyapunov analysis, contraction, dissipativity) might be used to study and design optimization algorithms certifying global (or semi-global) convergence in both continuous and discrete-time frameworks. The reader is referred, for instance, to [21, 30, 35, 11] for designs in discrete-time and to [17, 32, 19, 15] for designs in continuous-time. However, these are just non-exhaustive results owing to the problem's wide nature. Anyway, recent results [19] have interestingly shown that the classical primal-dual algorithm cannot achieve asymptotic convergence to the optimal (time-varying) point in the presence of time-varying constraints. Moreover, most of the existing results in such a context ensure only the boundedness of the solution provided by the optimization algorithm but not its asymptotic convergence to the optimum (see, e.g., [13, 31, 40, 26]). An exception is given for instance in [22], but with a methodology that is not computationally easy to apply. On the other hand, some preliminary results have recently shown that the use of internal-model theory plays a key role in the construction of algorithms with asymptotic convergence guarantees, see, e.g. [29, 8]. The specific contribution of this paper relies on designing an original optimization algorithm solving a convex optimization problem in the presence of periodic constraints with a known period.

<footnote><sup>a</sup>Université Lyon 1, Villeurbanne, France – CNRS, UMR 5007, LAGEPP, France.</footnote>

<footnote><sup>b</sup>University of Rome Tor Vergata, Electronic Engineering Department, Via del Politecnico 1, 00133, Rome, Italy.</footnote>

<footnote>*Research partially funded by the ANR project Alligator (ANR-22-CE48-0009-01).</footnote>

<footnote>Email addresses: daniele.astolfi@univ-lyon1.fr (Daniele Astolfi), verrelli@ing.uniroma2.it (Cristiano M. Verrelli)</footnote>
The design is restricted to the class of functions whose mini-
mizer vector exhibits a constant component within the kernel
space of the horizontal matrix characterizing the affine con-
straint. Even though similar ideas have also been recently pub-
lished in the context of quadratic optimization problems – see
[7] –, a discrete-time algorithm is designed there within a dif-
ferent framework. Here, the design does not rely on a finite-
dimensional internal model principle [10, 34, 39]. Its use is just
hinted, in the last part of the paper, as an approximation of the
infinite-dimensional system implementing the delay operator.
Instead, the design genuinely resorts to the RC tools<fnref n="1" />, see, e.g.,
[12, 3, 9, 36, 37, 27, 38], where a delay is used as a universal pe-
riodic signal generator (see the related discussion in [3])<fnref n="2" />. Two
algorithms are originally presented, which refer to currently
separate branches of the RC theory: the former resorts to the RC
theory of [9, 37]<fnref n="3" /> with a grain of structural originality regarding
the use of the saturation function, whereas the latter resorts to
the PDE transport-equation-based theory of RC of [3]. Both of
them, though exhibiting an increasing level of implementation
complexity while passing from a point-wise to a distributed use
of the delay, turn out to be natural extensions of the standard
primal-dual algorithm acting in the constant constraint scenario
and being typically described by a simple integral action. In
particular, the dual variable dynamics are modified into a single
delay in the first case and into a PDE transport-equation in the
second one. Global asymptotic convergence to the optimum is
analyzed via a Lyapunov analysis, while simulation results il-
lustrate the effectiveness of the proposed approaches. In light
of the previous considerations, an additional contribution of the
paper, which is hard to underestimate, relies on the simultane-
ous use of the two main different RC approaches in the litera-
ture to the same optimization problem. Indeed, the same field of
application allows for a direct comparison of the convergence
properties, the algorithm features, and the conditions at the root
of the different designs. Nevertheless, the crucial role of the
saturation function in its two variants (suitably included to ad-
dress robustness issues in real applications) is highlighted when
a common view is adopted for both the point-wise-delayed so-
lution and the distributed-delayed solution to the problem. Fi-
nally, the derivation of a finite-dimensional spectral approxi-
mation of the distributed-delayed solution introduces a further
connection with the renowned adaptive learning control.

Notation

We set $\mathbb{R}_+ = [0, \infty)$. Given a vector $x \in \mathbb{R}^l$, its components
are denoted by $x_1, \dots, x_l$. On the other hand, given a positive

real $T$, the vector $\pi(t)$ is $T$-periodic if

$$
\pi(t+T) - \pi(t) = 0 \quad (1)
$$

for any $t \ge -T$. Given an $m \times n$ matrix of rank $m \le n$, we denote with $a_{\max}$, $a_{\min} > 0$ the positive real numbers satisfying $a_{\min}I \le AA^\top \le a_{\max}I$. Furthermore, the scalar saturation function $sat_M(q)$ is the continuous odd increasing function of the scalar $q \in \mathbb{R}$ ($M > 0$) that satisfies $sat_M(q) = q$, $q \in [0, M]$, $sat_M(q) = M$, $q \ge M$, while the saturated vector $sat_M(x)$ ($x = [x_1, \dots, x_p]^\top$) is the vector whose elements are $sat_M(x_j)$, $j = 1, \dots, p$. Given a function $w : (t, x) \in \mathbb{R}^+ \times [0, T] \mapsto w(t, x) \in \mathbb{R}$, the notation $\partial_t w$ (resp. $\partial_x w$) denotes the partial derivative of $w$ with respect to the variable $t$ (resp. with respect to the variable $x$). We keep the notation for the weak and the strong definition of partial derivatives. When a function $w$ depends only on the variable of the time $t$, we denote its derivative $\dot{w}$. Finally, we denote as $L^2(0, 1)$ the Hilbert space of real-valued square-integrable functions over the interval $(0, 1)$ and as $H^1(0, 1)$ the Hilbert space of real-valued absolutely continuous functions over $[0, 1]$ with square-integrable derivative.

2. Problem formulation

In this section, the main optimization problem to be ad-
dressed in the paper, namely, problem $\mathcal{P}_t$, is formulated.

2.1. Definitions

For the sake of clarity and exhaustiveness, we start from the
classical problem (see problem $\mathcal{P}$, of which problem $\mathcal{P}_t$ con-
stitutes a time-varying generalization). Now, given a function
$f: \mathbb{R}^n \to \mathbb{R}^n$, the column vector $\nabla f$ denotes its gradient, that
is, $\nabla f(x) := \frac{\partial f}{\partial x}(x)$. We first recall the definition of strongly
convex function.

**Definition 1.** The function $f: \mathbb{R}^n \to \mathbb{R}$ is $\rho$-strongly convex if

$$
(\nabla f(x) - \nabla f(\hat{x}))^{\top} (x - \hat{x}) \geq \rho \|x - \hat{x}\|^2 \quad (2)
$$

for any pair $x, \hat{x} \in \mathbb{R}^n$.

Being inspired by [19], we first formulate the following
equality-constrained optimization problem $\mathcal{P}$.

**Definition 2. (Problem P):** Dynamically reconstruct, without statically computing Lagrange multipliers and/or the pseudo-inverse $A^+$, the point that minimizes a sufficiently smooth (at least continuously differentiable), time-invariant, known strongly-convex vector function $f: \mathbb{R}^n \to \mathbb{R}$:

$$
\min_{x \in \mathbb{R}^n} f(x) \tag{3}
$$

subject to the affine (convex) constraint:

$$
Ax = b \tag{4}
$$

with $A \in \mathbb{R}^{m \times n}$ being a known constant horizontal $m$-rank matrix ($m \le n$) and $b \in \mathbb{R}^m$ being a known constant vector.

<footnote n="1">¹When linear-parameterization-inducing functional approximations are adopted to reproduce the uncertain periodic signal to be rejected at steady-state, the RC is often referred to as *adaptive learning control* [34] (see also the related references therein).</footnote>

<footnote n="2">²The RC is often referred to as *repetitive learning control* (see [9, 37, 38] for instance) when such a delay is viewed as a part of a learning estimation scheme that can provide the estimate of the unknown periodic signal to be rejected at steady-state.</footnote>

<footnote n="3">³The saturated nature of the algorithm is crucial to preserve robustness in real practice and thus avoid asymptotical drifts: it plays the same role as the projection algorithm within the adaptive control scenario (see [16] for the design and the resulting convergence properties with no saturations).</footnote>
Under the presented restrictions, Problem $\mathcal{P}$ has a global unique minimum over $\mathbb{R}^n$. Furthermore, any $x_*$ satisfying (3) can be decomposed as

$$x_* = A^\top \eta_1 + \eta_2 \quad (5)$$

where $\eta_2 \in \text{Ker}(A)$ and $\eta_1$ is given by

$$\eta_1 = (AA^\top)^{-1}b. \quad (6)$$

Expression (5) is thus equivalent to

$$x_* = A^+b + (I - A^+A)w \quad (7)$$

for any $w \in \mathbb{R}^n$, where $A^+ = A^\top(AA^\top)^{-1}$ is the (Moore-Penrose) pseudo-inverse of $A$. The $f$-minimizer point is thus the point $x_m = A^+b + (I - A^+A)w_m$ with $w_m$ satisfying

$$w_m \in \left\{ \arg\min_{w \in \mathbb{R}^n} f(A^+b + (I - A^+A)w) \right\}. \quad (8)$$

Now, the following time-varying variant $\mathcal{P}_t$ of the above Problem $\mathcal{P}$ is formulated and solved hereafter. Such a variant of Problem $\mathcal{P}$ is characterized by the presence of a periodic time function being added to the constant vector $b$. The technical goal is then to dynamically reconstruct the time-varying vector $x_m(t)$ that minimizes a special class of strongly convex functions. Such a minimizing vector $x_m(t)$ might be viewed, in applications, as the reference provided by a supervisor to a lower control system.

**Definition 3.** (Problem $\mathcal{P}_t$): Given $A \in \mathbb{R}^{m \times n}$, $b \in \mathbb{R}^m$ as in Problem $\mathcal{P}$ and given a sufficiently smooth (at least continuously differentiable) $T$-periodic time function $\pi(\cdot) : \mathbb{R} \to \mathbb{R}^m$ (available at run-time), dynamically reconstruct, without statically computing Lagrange multipliers and/or the pseudo-inverse $A^+$, the $T$-periodic $f$-minimizer of a sufficiently smooth (at least continuously differentiable), time-invariant, known strongly-convex vector function $f: \mathbb{R}^n \to \mathbb{R}$:

$$\forall t \in \mathbb{R}_+, \quad \min_{x(t) \in \mathbb{R}^n} f(x(t)) \quad (9)$$

subject to the affine (convex) time-varying constraint:

$$Ax(t) = b + \pi(t) \quad (10)$$

under the following assumption $\mathcal{A}_k$: the $A$-Kernel-component of the $f$-minimizer vector is constant.

First notice how $\mathcal{P}_t$ reduces to $\mathcal{P}$ when $\pi(t) = 0$ and constraint (10) reduces, accordingly, to be time-invariant, i.e., as in (4). As a consequence, the algorithm solving $\mathcal{P}_t$ and identifying the $f$-minimizer will also solve $\mathcal{P}$ in such a case.

## 2.2. Class of functions and constraints

On the other hand, assumption $\mathcal{A}_k$ restricts the class of admissible constraints and $f$-functions in this paper. Indeed, in the reminder of this paper we will consider the class of functions: $f(x) = f_y(y) + f_z(z)$, where $f_y$ and $f_z$ are $\rho$-strongly convex,

$y = [x_1, \dots, x_r]^\top$ ($r < m$), $z = [x_{r+1}, \dots, x_n]^\top$ and of the affine (convex) time-varying constraint:

$$Ax = \left( \frac{A_1}{\mathbb{O}} \middle| \mathbb{O} \right) \quad x = \left( \frac{b_{p1} + \pi_p(t)}{b_{p2}} \right) = b + \pi(t) \quad (11)$$

where $\pi_p(\cdot) : \mathbb{R} \to \mathbb{R}^r$ is a sufficiently smooth (at least continuously differentiable) $T$-periodic time function (available at run-time); the rank of the $m \times n$ matrix $A$ is $m$, the $r \times r$ matrix $A_1$ is of full rank, the rank of the $(m-r) \times (n-r)$ matrix $A_2$ is $m-r$, $b = [b_{p1}^\top, b_{p2}^\top]^\top$ ($b_{p1} \in \mathbb{R}^r$, $b_{p2} \in \mathbb{R}^{m-r}$) is constant and $\pi(t) = [\pi_p(t)^\top, 0^\top]^\top$. The following illustrative scenarios shed light on such a framework.

### 2.2.1. Illustrative scenario I

Minimize the $\rho$-strongly convex function $f(x) = f_y(y) + f_z(z)$ in the sense of (9)-(10), with $y = [x_1, \dots, x_r]^\top$ ($r < m$), $z = [x_{r+1}, \dots, x_n]^\top$ and assume that there exists a non-singular $m \times m$ matrix $T_c$ such that

$$T_c(b + \pi(t)) = \left( \frac{b_{p1} + \pi_p(t)}{b_{p2}} \right), \quad T_c A = \left( \frac{A_1}{\mathbb{O}} \middle| \mathbb{O} \right) \quad (12)$$

with $A_1, A_2, b_{p1}, b_{p2}, \pi_p(\cdot)$ owning the previously introduced properties. Then the minimization problem is the one described by (9)-(11). As an example consider $f(x(t)) = \|x(t)\|^2$ ($n = 3$) and the constraint

$$\left( \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 0 & 1 \end{array} \right) x(t) = \left( \begin{array}{c} \sin(t) \\ 1 + \sin(t) \end{array} \right)$$

under $T_c = [1, 0; -1, 1]$. Notice how (invertible) changes of coordinates disjointly regarding $y$ and $z$ (thus preserving the block diagonal form of $A$ in (11)) and letting $f_y(\cdot), f_z(\cdot)$ be written in special forms are also allowed. This might be, for instance, the case of quadratic forms in $f_y(\cdot), f_z(\cdot)$ to be diagonalized.

### 2.2.2. Illustrative scenario II

Adopt the notation of this section and consider a vector $x$ of $n$ variables, a subset of them being constrained to satisfy the periodic invertible $A_1y = b_{p1} + w_p(t)$, and a subset of them being constrained to satisfy $A_2z = b_{p2}$ while having to simultaneously minimize a $\rho$-strongly convex function $f_z(z)$. Set any $\rho$-strongly convex function $f_y(y)$ and then determine, on the basis of the procedure solving Problem $\mathcal{P}_t$, the unique point $x_m(t) = [y_m^\top(t), z_m^\top]$ without resorting to the explicit computation of $A_1^{-1}$ and $A_2^+$ when the knowledge of $\nabla f_z(z)$, $A$, $b + \pi(t)$ (at runtime) and $T$ is available.

## 2.3. Minimizer form

For functions belonging to the class introduced in Section 2.2, a constant component belonging to the Kernel of $A$ is allowed to characterize the $f$-minimizer vector. Indeed, under the presented restrictions, any $x_*$ satisfying (10) and minimizing $f$ of (9) can be decomposed as

$$x_m(t) = A^+(b + \pi(t)) + (I - A^+A)w_m \quad (12)$$
with constant $w_m$ satisfying

$$w_m \in \left\{ \arg \min_{w \in \mathbb{R}^n} f(A^+(b+\pi(t)) + (\mathbb{I}-A^+A)w) \right\}. \quad (13)$$

In particular,

$$x_m(t) = \left( \begin{array}{c|c} A_1^{-1} & \mathcal{O} \\ \hline \mathcal{O} & A_2^{+} \end{array} \right) \left( \begin{array}{c} b_{p1} + \pi_p(t) \\ \hline b_{p2} \end{array} \right) + \left( \begin{array}{c|c} \mathcal{O} & \mathcal{O} \\ \hline \mathcal{O} & \mathbb{I} - A_2^{+} A_2 \end{array} \right) w_m.$$

### 3. Solution to Problem $\mathcal{P}$

In this section, with the aim of extending it to the time-varying scenario of Problem $\mathcal{P}_t$, we first recall the standard primal-dual algorithm (also called Arrow-Hurwicz-Uzawa flow [1]) providing Problem $\mathcal{P}$ with the solution:

$$\begin{aligned} \dot{\hat{x}} &= -\nabla f(\hat{x}) - A^\top \hat{\lambda} \\ \hat{\lambda} &= \gamma(A\hat{x} - b) \end{aligned} \quad (14)$$

where $\gamma$ is a positive design parameter. Its convergence properties are depicted by the following theorem.

**Theorem 1.** For any $\gamma > 0$ and for any initial conditions $(\hat{x}(0), \hat{\lambda}(0))$, the solutions to the algorithm (15) exponentially converge to the solution of Problem $\mathcal{P}$, namely, $\lim_{t \to \infty} \hat{x}(t) = x_m$ exponentially with a uniform rate over compact sets.

*Proof.* Consider the Lagrangian map $\mathcal{L}$ given by

$$\mathcal{L}(x, \lambda) = f(x) + \lambda^\top(Ax - b).$$

Recognize that $x_m$ satisfies, along with a constant $\lambda_m$, the system:

$$\begin{aligned} \nabla^\top f(x_m) + \lambda_m^\top A &= 0 \\ Ax_m - b &= 0 \end{aligned}$$

so that (14) can be rewritten in terms of $\tilde{x} = \hat{x} - x_m$ and $\tilde{\lambda} = \lambda_m - \hat{\lambda}$ as

$$\begin{aligned} \dot{\tilde{x}} &= -\nabla f(\tilde{x}) - A^\top \tilde{\lambda} + \nabla f(x_m) + A^\top \lambda_m \\ &= -(\nabla f(\tilde{x} + x_m) - \nabla f(x_m)) + A^\top \tilde{\lambda} \\ \dot{\tilde{\lambda}} &= -\gamma[(A\tilde{x} - b) - (Ax_m - b)] \\ &= -\gamma A \tilde{x}. \end{aligned}$$

Consider the quadratic Lyapunov function (and its time derivative along the trajectories of the above error system)

$$V(\tilde{x}, \tilde{\lambda}) = \tilde{x}^\top \tilde{x} + \frac{1}{\gamma} \tilde{\lambda}^\top \tilde{\lambda} \quad (15)$$

and use the property (2) of $f$ to recognize the uniform asymptotic stability of the origin of the above error system. In particular, $\tilde{x}$ converges to $x_m$, exponentially with a uniform rate over compact sets (see, e.g., [19] or adapt [34]). □

We remark, at this stage, that several modifications of this algorithm have been proposed, but it is out of the scope of this article to include a precise discussion in this section. We refer, for instance, to [19, 15] or references therein.

### 4. Point-wise-delayed solution to Problem $\mathcal{P}_t$

In this section, we extend the solution of the previous section to the periodic case characterizing Problem $\mathcal{P}_t$ in which the constraint is now time-varying and periodic. To this end, inspired by the literature on RC, we modify the scheme (14) by introducing an element able to generate $T$-periodic signals, and not only constant as in the case of an integral action. As a consequence, the $\hat{\lambda}$-dynamics in (14) is modified by introducing a delay equation as follows

$$\hat{\lambda}(t) = \hat{\lambda}(t-T) + \gamma(A\hat{x} - b - \pi(t)) \quad (16)$$

in which we recognize the classical delay element of the RC given by $\hat{\lambda}(t) = \hat{\lambda}(t-T)$ which generalizes the integral action $\hat{\lambda} = 0$ to any periodic signal (see [33]). In the subsequent section, we show how to adapt such an idea to the context of an optimization algorithm.

#### 4.1. A partially saturated version

By resorting to the RC tools (in light of the technical remark 3 of [37]), naturally modify (16) into $(\gamma$ is again a positive design parameter, here called learning gain, while $\text{sat}_M(\cdot)$ is the $M$-saturation function of [20] recalled in Section 1 and $\varphi(\cdot) : \mathbb{R}_+ \to [0, 1]$ is the continuity-preserving (continuous increasing) function in [37] such that $\varphi(0) = 0$, $\varphi(t) = 1$ for any $t \ge T)$:

$$\begin{aligned} \dot{\hat{x}} &= -\nabla f(\hat{x}) - A^\top \hat{\lambda} \\ \hat{\lambda}(t) &= \text{sat}_M(\hat{\lambda}(t-T)) + \gamma\varphi(t)(A\hat{x}(t) - b - \pi(t)) \quad (17) \\ \hat{\lambda}(r) &= 0, \quad \forall r \le 0. \end{aligned}$$

Then we recall from [20] the following:

**Lemma 1.** Given $M \in \mathbb{R}_+$, for any $\tau \in \mathbb{R}$ and for any $\sigma \in \mathbb{R}$ such that $|\sigma| \le M$:

$$(\text{sat}_M(\tau) - \sigma)^2 - (\tau - \sigma)^2 \le 0,$$

whose proof is also reported hereafter (it will constitute a base for a further modification).

*Proof.* Just use the following straightforward computation

$$(\text{sat}_M(\tau) - \sigma)^2 - (\tau - \sigma)^2 = (\text{sat}_M(\tau) - \tau)(\text{sat}_M(\tau) - 2\sigma + \tau), \\ \text{starting from which the three cases: } \tau > M, \tau \in [-M, M], \tau < -M \text{ can be examined separately. } \square$$

A variant of Lemma 1 is then introduced at this stage. It will be used in Section 5 to prove Theorem 4.

**Lemma 2.** Given $M \in \mathbb{R}_+$, for any $\tau \in \mathbb{R}$ and for any $\sigma \in \mathbb{R}$ such that $|\sigma| \le M$:

$$(\tau - \sigma)(\tau - \text{sat}_M(\tau)) \ge 0.$$
Proof. We are in the same scenario of Lemma 1. Directly from

$$(\mathrm{sat}_M(\tau) - \sigma)^2 - (\tau - \sigma)^2 \le 0,$$

add and subtract $\tau$ within the first term, take the square of the resulting binomial, and definitely get

$$(\mathrm{sat}_M(\tau) - \tau)^2 - 2(\tau - \sigma)(\tau - \mathrm{sat}_M(\tau)) \le 0$$

which the thesis follows on.

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

On the basis of the Lemma 1, we are able to prove the following theorem, not until we have introduced, for the sake of compactness, the following positive real:

$$M_0 := \max_{j=1,\dots,m} \left\{ \sup_{\tau \in [0,T]} \left| [\lambda_m(\tau) - (AA^\top)^{-1}\hat{\pi}(\tau)]_j \right| \right\}, \quad (18)$$

where $[\cdot]_j, j = 1, \dots, m$, denote the components of an $m$-dimensional vector. Such a positive real shall be actually used in the statements of Theorems 2-5 in the remainder of this paper.

**Theorem 2.** For any $\gamma > 0$, for any $M \ge M_0$, and for any initial condition $\hat{x}(0)$, the solutions to the algorithm (17) asymptotically (not uniformly w.r.t. the initial conditions) converges to the solution of Problem $\mathcal{P}_t$, i.e., $\lim_{t \to \infty} \hat{x}(t) = x_m(t)$.

Proof. Consider the Lagrangian map (at each $t$)

$$\mathcal{L}(x(t), \lambda(t)) = f(x(t)) + \lambda^\top(t)(Ax - b - \pi(t)).$$

Recognize that the $T$-periodic $x_m(t)$ [recall (1)] satisfies, along with the $T$-periodic $\lambda_m(t)$, the system:

$$\begin{aligned} \nabla^\top f(x_m(t)) + \lambda_m(t)^\top A &= 0 \\ Ax_m(t) - b - \pi(t) &= 0 \end{aligned} \quad (19)$$

so that (17), on the basis of (12), can be rewritten in terms of $\tilde{x} = \hat{x} - x_m$ and $\tilde{\lambda} = \lambda_m - (AA^\top)^{-1}\hat{\pi} - \tilde{\lambda}$ as (t-dependence is omitted where its use is obvious)

$$\begin{aligned} \dot{\tilde{x}} &= -\nabla f(\tilde{x}) - A^\top \tilde{\lambda} + \underbrace{\nabla f(x_m)}_{0} + A^\top \tilde{\lambda}_m - \underbrace{A^+ \pi}_{\tilde{x}_m} \\ &= -(\nabla f(\tilde{x} + x_m) - \nabla f(x_m)) + A^\top \tilde{\lambda} \\ \tilde{\lambda}(t) &= \mathrm{sat}_M(\tilde{\lambda}(t-T)) + \gamma \varphi(t) A \tilde{x}(t). \end{aligned} \quad (20)$$

Consider now the quadratic Lyapunov-like functional

$$\begin{aligned} V_p(t) &= \tilde{x}(t)^\top \tilde{x}(t) + \frac{1}{\gamma} \int_{t-T}^{t} \tilde{\lambda}_S(\tau)^\top \tilde{\lambda}_S(\tau) d\tau \\ \tilde{\lambda}_S &= \mathrm{sat}_M(\tilde{\lambda}) - \lambda_m + (AA^\top)^{-1} \hat{\pi}. \end{aligned} \quad (21)$$

By taking its time derivative along the trajectories of system (20), by using the property (2) of $f$, and by employing Lemma 1 with $M \ge M_0$ defined as in (18), finally establish the asymptotic convergence of $\hat{x}(t)$ to $x_m$ through Barbalat's lemma (see similar arguments in [27, 37, 9, 20] and notice the crucial use of the saturation function to get the boundedness of $\tilde{\lambda}(t)$). □

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

## Technical Remarks

1. The above result holds even in the case of no signal $\pi(t)$ affecting $b$ since any constant function is also $T$-periodic [recall (1)]. Indeed, a constant function is periodic of any positive period.

2. The Lyapunov-like functional (21), when no saturation is used, is a direct generalization of the Lyapunov function (15).

3. Assumption $\mathcal{A}_k$ has allowed $(AA^\top)^{-1}\hat{\pi}$ to appear within (20) multiplied by $A^\top$ so that its contribution can be successfully compensated by the learning estimate $\tilde{\lambda}$.

4. The arguments within the above proof of Theorem 2 are still valid even if the quadratic Lyapunov-like functional (21) – classically used in the literature of partially saturated learning estimation schemes – is replaced by the same Lyapunov-like functional

$$V_{ps}(t) = \tilde{x}(t)^\top \tilde{x}(t) + \frac{1}{\gamma} \int_{t-T}^{t} \tilde{\lambda}(\tau)^\top \tilde{\lambda}(\tau) d\tau$$

that the next subsection will resort to. In such a case, it suffices, when computing $\hat{V}_{ps}(t)$ to i) add and subtract $[\mathrm{sat}_M(\hat{\lambda}(t-T))]_j$ within the squared term $[(\lambda(t) - \hat{\lambda}_j(t))]_j^2$, ii) perform the square of the resulting binomial, iii) add and subtract $[\hat{\lambda}(t)]_j$ within twice the product of the two terms of the binomial and then use the explicit expression of $[\hat{\lambda}(t)]_j$ of (20), iv) use Lemma 1.

## 4.2. A fully saturated variant

The previous solution, though originally designed for the optimization problem under investigation, inherits from the quoted literature the drawback of saturating just the delayed estimate while leaving the error – driving the estimation process – unsaturated. This may lead to estimates $\hat{\lambda}(t)$ far beyond the saturation limits, especially during the transient, though the $\tilde{x}$-term within $\hat{\lambda}(t)$ acts as a stabilizing term in the $\tilde{x}$-dynamics. To overcome such a drawback, a fully saturated variant of the estimation scheme in (17), leading to the algorithm

$$\begin{aligned} \dot{\tilde{x}} &= -\nabla f(\tilde{x}) - A^\top \tilde{\lambda} \\ \tilde{\lambda}(t) &= \mathrm{sat}_M(\tilde{\lambda}(t-T) + \gamma \varphi(t)(A\hat{x}(t) - b - \pi(t))) \\ \tilde{\lambda}(r) &= 0, \quad \forall r \le 0, \end{aligned} \quad (22)$$

is innovatively presented as well, which positively answers the question in [20] regarding the existence of the proof for such a variant. To provide it, we originally state (and resort to) the following lemma, which extends Lemma 1 to the case of any $\sigma$.

**Lemma 3.** Given $M \in \mathbb{R}_+$, for any $\tau \in \mathbb{R}$ and any $\sigma \in \mathbb{R}$:

$$(\mathrm{sat}_M(\tau) - \sigma)^2 - (\tau - \sigma)^2 \le \max\{\sigma - M, 0\}^2 + \max\{-\sigma - M, 0\}^2.$$
Proof. Use the compact notation

$$
G(\tau, \sigma) = (\mathrm{sat}_M(\tau) - \sigma)^2 - (\tau - \sigma)^2.
$$

Then, for $|\sigma| \le M$, $G(\tau, \sigma) \le 0$ in accordance with Lemma 1.
On the other hand, for $\tau > M$, $G(\tau, \sigma)$ becomes

$$
G(\tau, \sigma)|_{\tau>M} \doteq G_{\tau M>} = (M - \sigma)^2 - (\tau - \sigma)^2
$$

so that, for $\sigma > M$, $G_{\tau M>} \le (M - \sigma)^2 = \max\{\sigma - M, 0\}^2$, while, for $\sigma < -M$, once $\tau$ is rewritten as $\tau = M + \epsilon$ with $\epsilon > 0$, $G_{\tau M>} = (M - \sigma)^2 - (M - \epsilon - \sigma)^2 = -2\epsilon(M - \sigma) - \epsilon^2 < 0$. Finally, for $\tau < -M$, $G(\tau, \sigma)$ becomes

$$
G(\tau, \sigma)|_{\tau<-M} \doteq G_{\tau M<} = (-M - \sigma)^2 - (\tau - \sigma)^2
$$

so that, for $\sigma < -M$, $G_{\tau M<} \le (-M - \sigma)^2 = \max\{-\sigma - M, 0\}^2$, while, for $\sigma > M$, once $\tau$ is rewritten as $\tau = -M - \epsilon$ with $\epsilon > 0$, $G_{\tau M<} = (-M - \sigma)^2 - (-M - \epsilon - \sigma)^2 = 2\epsilon(-M - \sigma) - \epsilon^2 < 0$.

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

Lemma 3 is at the root of the proof of the following theorem.

**Theorem 3.** For any $\gamma > 0$, for any $M \ge M_0$, and for any initial condition $\hat{x}(0)$, the solutions to the algorithm (22) asymptotically (not uniformly w.r.t. the initial conditions) converges to the solution of Problem $\mathcal{P}_t$, i.e., $\lim_{t \to \infty} \hat{x}(t) = x_m(t)$.

Proof. The proof is analogous to the proof of Theorem 2, except for the use of the quadratic Lyapunov-like functional (no saturation action is used within it as highlighted by Remark 4 of the previous subsection)

$$
V_{ps}(t) = \tilde{x}(t)^{\top} \tilde{x}(t) + \frac{1}{\gamma} \int_{t-T}^{t} \tilde{\lambda}(\tau)^{\top} \tilde{\lambda}(\tau) d\tau, \quad (23)
$$

which, for sufficiently large $M \ge M_0$, with $M_0$ given as in (18),
has a time derivative satisfying for $t \ge T$:

$$
\begin{aligned}
\dot{V}_{ps}(t) \leq & -2\rho\tilde{x}(t)^{\top}\tilde{x}(t) - 2\tilde{x}(t)^{\top}A^{\top}\tilde{\lambda}(t) + \frac{1}{\gamma}\tilde{\lambda}(t)^{\top}\tilde{\lambda}(t) \\
& - \frac{1}{\gamma}\tilde{\lambda}(t-T)^{\top}\tilde{\lambda}(t-T).
\end{aligned}
$$

Next, by adding and subtracting the term

$$
\frac{1}{\gamma}(\tilde{\lambda}(t) - \gamma A \tilde{x}(t))^{\top} (\tilde{\lambda}(t) - \gamma A \tilde{x}(t))
$$

to the previous inequality, we obtain

$$
\begin{aligned}
\dot{V}_{ps}(t) &\le -2\rho\tilde{x}(t)^{\top}\tilde{x}(t) - 2\tilde{x}(t)^{\top}A^{\top}\tilde{\lambda}(t) + \frac{1}{\gamma}\tilde{\lambda}(t)^{\top}\tilde{\lambda}(t) \\
&\quad - \frac{1}{\gamma}(\tilde{\lambda}(t) - \gamma A\tilde{x}(t))^{\top}(\tilde{\lambda}(t) - \gamma A\tilde{x}(t)) + \frac{1}{\gamma}E(t) \\
&\le -\tilde{x}(t)^{\top}(2\rho\mathbb{I} + \gamma A^{\top}A)\tilde{x}(t) + \frac{1}{\gamma}E(t)
\end{aligned}
$$

where the term $E(t)$ is defined as

$$
\begin{aligned}
E(t) = & (\tilde{\lambda}(t) - \gamma A \tilde{x}(t))^{\top} (\tilde{\lambda}(t) - \gamma A \tilde{x}(t)) \\
& - \tilde{\lambda}(t-T)^{\top} \tilde{\lambda}(t-T).
\end{aligned}
$$

Using the expression for the estimate $\hat{\lambda}(t)$ in (22), each
component $E_j(t)$ of such a term $E(t)$ reads $E_j(t) =$

$(\mathrm{sat}_M(s_j(t) + B_j(t)) - B_j(t))^2 - s_j(t)^2$ in terms of the com-
pact notation $s_j(t) = [\tilde{\lambda}(t - T)]_j$ and $B_j(t) = [-\lambda_m(t) +
(AA^\top)^{-1}\hat{\pi}(t) - \gamma\varphi(t)A\tilde{x}(t)]_j$. As a consequence, by repeatedly
using Lemma 3 (use $\tau$ in place of $s_j(t) + B_j(t)$ and $\sigma$ in place
of $B_j(t))$ and accounting for the action of the max-function, we
obtain $E(t) \le \gamma^2 \tilde{x}^\top A^\top A \tilde{x}(t)$, leading to $\dot{V}_{ps}(t) \le -2\rho \tilde{x}(t)^\top \tilde{x}(t)$.
From now on, arguments similar to the ones used in the proof
of Theorem 2 can be used to conclude the desired result.
$\square$

Notice how, differently from the arguments depicted in the
technical remark 4 of this section (as well as in the proof of
Theorem 2), the fully saturated nature of the estimation scheme
is proved at the expense of the negative term $-\gamma^2 \tilde{x}^\top A^\top A \tilde{x}(t)$
within the derivative of $V_{ps}(t)$.

**5. Distributed-delayed solution to Problem** $\mathcal{P}_t$

A particular feature of the schemes (17), (22), is that only the
information at $\hat{\lambda}(t)$ and $\hat{\lambda}(t-T)$ is needed for implementation.
However, in the recent literature on RC, different solutions have
been proposed, which resort to all the distributed knowledge of
$\hat{\lambda}(\tau)$ with $\tau \in [t, t+T]$, see, e.g. [3]. This can be formally intro-
duced by adopting a different notation for the delayed equation
(16), which is used to describe such an equation via a PDE,
and, more precisely, via a transport equation. In this light, the
designer can be provided with additional degrees of freedom
for the overall scheme, with well-known spectral properties of
transport equations naturally leading to finite-dimensional ap-
proximations on the basis of Fourier representations (the reader
is referred to Section 5.3 below).

5.1. A short historical detour on RC literature

As a matter of fact, a delayed equation $\hat{\lambda}(t) = \hat{\lambda}(t-T)$ can
be equivalently represented by the following PDE:

$$
\partial_t \hat{\lambda}(t,s) = -\partial_s \hat{\lambda}(t,s), \quad \hat{\lambda}(t,0) = \hat{\lambda}(t,T),
$$

with $\hat{\lambda}(t,s) \in \mathbb{R}^m$, for all $(t,s) \in \mathbb{R}_+ \times [0,T]$. Based on such a
representation, we present in the next subsection a distributed-
delayed version of an algorithm solving the problem $\mathcal{P}_t$.

5.2. Distributed solution

By resorting to a distributed use of the memory [2], in place
of (17) (or (22)) relying on point-wise use of the delay, we de-
sign the scheme:

$$
\begin{align}
\dot{\hat{x}}(t) &= -\nabla f(\hat{x}(t)) - A^{\top} \int_{0}^{T} e^{sT} \hat{\lambda}(t,s) ds \nonumber \\
\partial_t \hat{\lambda}(t,s) &= -\partial_s \hat{\lambda}(t,s) + \gamma e^{sT} (A\hat{x}(t) - b(t) - \pi(t)) \tag{24} \\
&\quad - \mu \, dz_M(\hat{\lambda}(t,s)) \nonumber
\end{align}
$$

$$
\hat{\lambda}(t,0) = \hat{\lambda}(t,T), \quad \hat{\lambda}(0,x) = \eta_0(x),
$$

for all $(t,s) \in \mathbb{R}_+ \times [0,T]$, $\hat{\lambda}(t,s) \in \mathbb{R}^m$, where the $dz_M$ is the dead-zone function defined as

$$
\mathrm{d}z_M(s) = - \mathrm{sat}_M(s) + s \quad (25)
$$
which is added to obtain an anti-wind-up effect (similarly to saturation in the equation (17)), and with $\mu > 0$ and $M > 0$. The following theorem is then established.

**Theorem 4.** For any $\gamma > 0$, $\mu > 0$, for any $M \ge M_0$, and for any initial conditions $(\hat{x}_0, \eta_0) \in \mathbb{R}^n \times \mathbb{H}^1(0, T)$, solutions to the algorithm (24) asymptotically converge to the solution of Problem $P_w$, namely, $\lim_{t \to \infty} \hat{x}(t) = x_m(t)$.

*Proof.* Since $\lambda_m, \pi$ are T-periodic signals, we can represent those trajectories as the solution to

$$ 
\begin{aligned}
\partial_t \Lambda(t,s) &= -\partial_s \Lambda(t,s) & \partial_t \Pi(t,s) &= -\partial_s \Pi(t,s) \\
\Lambda(t,0) &= \Lambda(t,T) & \Pi(t,0) &= \Pi(t,T) \\
\Lambda(0,s) &= \Lambda_0(s) & \Pi(0,s) &= \Pi_0(s)
\end{aligned}
\quad (26) $$

with $\lambda_m(t) = \int_0^T e^{sT} \Lambda(t,s)ds$ and $\pi(t) = \int_0^T e^{sT} \Pi(t,s)ds$. This can be shown by using Lemma 1 in [3]. Furthermore, $\pi$ being $C^1$, for solutions $\Pi$ in $H^2(0, T)$ we also get $\partial_{tt}\Pi(t,s) = -\partial_{ts}\Pi(t,s)$, $\partial_t\Pi(t,0) = \partial_t\Pi(t,T)$, $\pi(t) = \int_0^T e^{sT} \partial_t\Pi(t,s)ds$. As a consequence, in accordance with the first part of the proof of Theorem 2 (e.g., (18)-(19)), we write

$$ \dot{x}_m(t) = -\nabla f(x_m(t)) - A^\top \int_0^T e^{sT} \Lambda(t,s)ds + A^+ \int_0^T e^{sT} \partial_t \Pi(t,s)ds. $$

Now consider the following change of coordinates

$$ 
\begin{aligned}
\hat{x} \mapsto \tilde{x} &:= \hat{x} - x_m, \\
\tilde{\lambda}(\cdot, s) \mapsto \tilde{\lambda}(\cdot, s) &:= \tilde{\lambda}(\cdot, s) - \Lambda(\cdot, s) + (AA^\top)^{-1} \partial_t \Pi(\cdot, s),
\end{aligned}
\quad (27) $$

from which we obtain, following steps similar to the previous proofs,

$$ 
\begin{aligned}
\dot{\tilde{x}} &= -(\nabla f(\tilde{x} + x_m) - \nabla f(x_m)) - A^\top \int_0^T e^{sT} \tilde{\lambda}(t,s)ds \\
\partial_t \tilde{\lambda}(t,s) &= -\partial_s \tilde{\lambda}(t,s) + \gamma e^{sT} A \tilde{x} - \mu \, dz_M(\tilde{\lambda}(t,s) + c(t,s)) \\
\tilde{\lambda}(t,0) &= \tilde{\lambda}(t,T)
\end{aligned}
\quad (28) $$

with the compact notation $c(t,s) = \Lambda(t,s) + (AA^\top)^{-1}\partial_t\Pi(t,s)$. It can be observed that system (28) can be seen as the interconnection between two subsystems with no throughout terms like the ones in [12]. Now consider the Lyapunov functional

$$ V(\tilde{x}, \tilde{\lambda}) = \tilde{x}^\top \tilde{x} + \frac{1}{\gamma} \int_0^T \tilde{\lambda}(\cdot, s)^\top \tilde{\lambda}(\cdot, s) ds $$

which is equivalent to the standard norm in $\mathbb{R}^n \times L^2(0, T)$. Using the convexity property of $f$, the computation of its derivative along the solutions of (28) gives

$$ 
\begin{aligned}
\dot{V} &= -2\rho \tilde{x}^\top \tilde{x} - 2\tilde{x}^\top A^\top \int_0^T e^{sT} \tilde{\lambda}(t,s)ds \\
&\quad + \frac{2}{\gamma} \int_0^T \tilde{\lambda}(t,s)^\top (-\partial_s \tilde{\lambda}(t,s)ds + \gamma e^{sT} A \tilde{x})ds \\
&\quad - \frac{2\mu}{\gamma} \int_0^T \tilde{\lambda}(t,s)^\top dz_M(\tilde{\lambda}(t,s) + c(t,s))ds \\
&= -2\rho |\tilde{x}(t)|^2 - \frac{1}{\gamma} \left[ \tilde{\lambda}(t,s)^\top \tilde{\lambda}(t,s) \right]_0^T \\
&\quad - \frac{2\mu}{\gamma} \int_0^T \tilde{\lambda}(t,s)^\top dz_M(\tilde{\lambda}(t,s) + c(t,s))ds.
\end{aligned}
$$

Using the boundary condition we have $[\tilde{\lambda}(t,s)^\top \tilde{\lambda}(t,s)]_0^T = 0$. Furthermore, the selection of $M$ as in the statement of the theorem implies $\lambda(t,s)^\top dz_M(\tilde{\lambda}(t,s) + c(t,s)) \ge 0$ for all $t,s$, in accordance with Lemma 2. As a consequence, we finally obtain $\dot{V} \le -2\rho |\tilde{x}(t)|^2$. This implies also that $\dot{V} \le 0$. Note that owing to the dead-zone function, $\tilde{\lambda}(t,s)$ is always bounded and so is $\hat{x}(t)$. Furthermore, considering initial conditions in $H^1(0, T)$ allows us to conclude that $\hat{x}$-solutions are also uniformly continuous. Barbalat's lemma allows us to establish the thesis. Additional technical details can be found in [3]. □

### 5.3. Towards Adaptive Learning Control & Harmonic Regulation through Spectral Approximation

One interesting aspect is that the proposed algorithms are infinite-dimensional and therefore are required to be approximated via a numerical scheme (no asymptotic convergence to zero can be established in such a case). While a Padé-based approximation of the point-wise delay of Section 4 has been already dealt with in [33, 9], for what the algorithm of this section is concerned a possible approximating manner (here neglect the saturation) is by using a finite-dimensional ODE approximating the delay equation. Since the spectrum of a delay equation

$$ \partial_t \tilde{\lambda}(t,s) = -\partial_s \tilde{\lambda}(t,s), \quad \lambda(t,0) = \lambda(t,T) $$

is given by $0 \cup \{\pm ik\frac{2\pi}{T}\}$, $k \in \mathbb{N}$, one can approximate such an operator by projecting on the Fourier basis, as explained in [6], and already explored in [23]. Indeed, a continuous-time approximate solution to the infinite-dimensional schemes proposed so far can be also provided in light of [34] and of the technical comments within [39] (see also references therein) connecting the Adaptive Learning Control (recall footnote 1) to the internal-model-based design in periodic frameworks (see also [37]), also denoted as harmonic regulation [4]. In particular, we introduce the following set of matrices

$$ 
\begin{aligned}
S &= \frac{2\pi}{T} \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}, & C &= \begin{pmatrix} 1 & 0 \end{pmatrix} \\
\bar{S} &= \text{blkdiag}(0, S, 2S, \dots, NS), & \bar{C} &= \begin{pmatrix} 1, & \dots, & C \end{pmatrix}_{\substack{N \text{ times} \\ \vdots}} \\
\Phi &= \text{blkdiag}(\bar{S}, \dots, \bar{S}), & \Gamma &= \text{blkdiag}(\bar{C}, \dots, \bar{C})
\end{aligned}
$$

for some $N \in \mathbb{N}$, $N \ge 1$. Note that $\Phi \in \mathbb{R}^{(2N+1)m \times (2N+1)m}$ and $\Gamma \in \mathbb{R}^{m \times (2N+1)m}$. Then, we propose the following algorithm:

$$ 
\begin{aligned}
\dot{\hat{x}} &= -\nabla f(\hat{x}) - A^\top \Gamma \hat{\lambda} \\
\hat{\lambda} &= \Phi \hat{\lambda} + \gamma \Gamma^\top (A\hat{x} - b - \pi(t)) - \mu \, dz_M(\hat{\lambda})
\end{aligned}
\quad (29) $$

where $\lambda \in \mathbb{R}^{(2N+1)m}$ and where $\gamma$ is the tuning parameter. The algorithm (29) thus provides an approximation of the previous algorithm (17) and (24) in the harmonic sense described below.

**Theorem 5.** For any $\gamma > 0$, $\mu > 0$, for any $M \ge M_0$, the algorithm (29) converges exponentially to a periodic solution satisfying

$$ \lim_{t \to \infty} \int_t^{t+T} (A\hat{x}(t) - b - \pi(t)) \exp\left(\frac{2k\pi}{T}t\right) dt = 0, $$
for all $k = 0, 1, \dots, N$.

*Proof.* Let $M$ be chosen large enough as in the previous theorems. Introduce the following compact notation $\hat{z} = (\hat{x}, \hat{\lambda})$ and rewrite system (29) as $\hat{z} = F(\hat{z}) + G(b + \pi(t))$ with

$$F(\hat{z}) = \begin{pmatrix} -\nabla f(\hat{x}) - A^{\top}\Gamma\hat{\lambda} \\ \Phi\hat{\lambda} + \gamma\Gamma^{\top}A\hat{x} - \mu \, \mathrm{d}z_M(\hat{\lambda}) \end{pmatrix}, \quad G = \begin{pmatrix} 0 \\ -\gamma\Gamma^{\top} \end{pmatrix}.$$

It can be noticed that the $\hat{x}$ and $\hat{\lambda}$ dynamics have a (incrementally) passive interconnection. We will show that this system is incrementally input-to-state stable ($\delta$ISS) with respect to the input $b + \pi(t)$ by constructing a $\delta$ISS Lyapunov function. To this end consider any two pair of trajectories $\hat{z}' = (\hat{x}', \hat{\lambda}')$, and $\hat{z}'' = (\hat{x}'', \hat{\lambda}'')$ of system (29), corresponding to different inputs $d'(t) = b' + \pi'(t)$ and $d''(t) = b'' + \pi''(t)$, and let $\tilde{z} = (\tilde{x}, \tilde{\lambda}) = (\hat{x}' - \hat{x}'', \hat{\lambda}' - \hat{\lambda}'')$, $\tilde{d} = d' - d''$. With some abuse of notation, we also denote $\tilde{F}(\tilde{z}) = F(\hat{z}') - F(\hat{z}'')$. Then, let $P = \mathrm{diag}(I_n, \gamma^{-1}I_{(2N+1)m})$ and consider first the incremental Lyapunov function $V = \tilde{z}^{\top}P\tilde{z}$. Note that such a Lyapunov function is “equivalent” to the Lyapunov functions used in the previous sections. Using the fact that $\Phi + \Phi^{\top} = 0$ and the monotonicity of the dead-zone function, it can be verified that

$$\dot{V} = 2\tilde{z}^{\top}P\tilde{F}(\tilde{z}) + G\tilde{d}(t) \leq -2\rho|\tilde{x}|^2 + 2\tilde{z}^{\top}PG\tilde{d}(t).$$

This Lyapunov function is not a strict $\delta$ISS Lyapunov function, yet it shows that whenever $\tilde{d} = 0$, then $\tilde{x}$ converges to zero. Furthermore, one can show that also $\tilde{\lambda}$ converges to zero by noting that the pair $(S_N, G_N)$ is observable by construction. As a consequence, we strictify the Laypunov function $V$ by following the same construction as in [14, Theorem 2]. In particular, select any $\Upsilon$ such that $(\Phi + \Upsilon\Gamma)$ is Hurwitz and let $R = R^{\top} \succ 0$ satisfy

$$R(\Phi + \Upsilon\Gamma) + (\Phi + \Upsilon\Gamma)^{\top}R \leq -I_{(2N+1)m}.$$

By defining $Q = \mathrm{blkdiag}(0_{n \times n}, R)$, consider the (strict) Lyapunov function $W(\tilde{z}) = V(\tilde{z}) + \varepsilon\tilde{z}^{\top}Q\tilde{z}$. Simple computations show that

$$\dot{W} = 2\tilde{z}^{\top}(P + \varepsilon Q)(\tilde{F}(\tilde{z}) + G\tilde{d}(t)) \leq -\varepsilon|\tilde{z}|^2 + \tilde{z}^{\top}(P + \varepsilon Q)G\tilde{d},$$

showing that the system is incrementally input-to-state stable with respect to $\tilde{d}$. The claim of the theorem follows by direct application of Proposition 1 in [24]. □

### Technical Remarks

1. If $\pi$ is $C^1$ with finite spectrum, that is, its Fourier series is given by a finite number $n_\pi$ of components $\pi(t) = \sum_{k=1}^{n_\pi} \pi_k \cos(2k\pi t + \varphi_k)$, then the result of Theorem 5 implies that $A\hat{x} - b - \pi = 0$ at steady-state if $N \geq n_\pi$, and therefore, by uniqueness of the solutions we obtain $\hat{x} = x_m$.

2. In the general case in which the signal $\pi$ has a Fourier series with an infinite number of harmonics, the result of Theorem 5 states that all the harmonics up to $N$ in the signal $A\hat{x} - b - \pi = 0$ are compensated at steady-state. If all signals are smooth enough, this further implies that $\hat{x}$ converge arbitrarily close to $x_m$ in a $L^2$ sense by increasing the

parameter $N$, namely, there exists $c > 0$ independent of $N$ such that, for any $N \geq 1$

$$\lim_{t \to \infty} \int_t^{t+T} |\hat{x}(t) - x_m(t)|^2 dt \leq \frac{c}{N}$$

for any initial condition $(\hat{x}(0), \hat{\lambda}(0)) \in \mathbb{R}^n \times \mathbb{R}^{2N+1}$. This can be shown combining arguments similar to [5, 6]. In that case, the bound is also uniform over compact sets, see, e.g. [23, 5, 39].

3. Due to the integral action incorporated in the $\hat{\lambda}$-dynamics, the algorithm (29) still recovers the properties of the standard dual-primal dynamics (14) whenever $\pi(t) = 0$.

## 6. A theoretical comparison

The theoretical comparison of this section concerns the structure of the solutions presented in the previous two sections and provides, for the first time in the literature, the tools to connect them. This is an additional contribution of this paper, which allows the reader to successfully flip through different branches of the literature regarding the RC. Specifically, the point-wise characterization issue is here addressed in the context of PDEs. Accordingly, as discussed above, the main link relies on recognizing that the delayed equation in (16) is actually a PDE, and more precisely a transport equation. Hence, using the notation proposed in [3], such a delayed equation is equivalently represented as (omit the saturation for the sake of simplicity)

$$\begin{aligned} \partial_t \hat{\lambda}(t,s) &= -\partial_s \hat{\lambda}(t,s) \\ \hat{\lambda}(t,0) &= \hat{\lambda}(t,T) + \gamma(A\hat{x} - b - \pi(t)). \end{aligned}$$

This allows us to highlight the special link existing between the two solutions so far presented in the paper. If the saturation function and $\varphi(t)$ are neglected, then the solution of Section 5 and the related stability arguments reduce to the one of Section 4 when the distributed use of the memory becomes point-wise. In other words, by preserving the notation of Section 5, algorithm (22) or (17) reduces to

$$\begin{aligned} \dot{\hat{x}}(t) &= -\nabla f(\hat{x}(t)) - A^{\top}\hat{\lambda}(t,0) \\ \partial_t \hat{\lambda}(t,s) &= -\partial_s \hat{\lambda}(t,s) \\ \hat{\lambda}(t,0) &= \hat{\lambda}(t,T) + \gamma(A\hat{x} - b - \pi(t)) \end{aligned} \quad (30)$$

with $\hat{\lambda}(0, s) = \hat{\lambda}_0(s)$. On the other hand, since the computations of the previous section lead to the condition (26) with $\lambda_m(t) = \Lambda(t, 0)$, $\pi(t) = \Pi(t, 0)$, and since, $\pi$ being $C^1$, for solutions $\Pi$ in $H^2(0, T)$ the following identities $\partial_t\Pi(t, s) = -\partial_s\Pi(t, s)$, $\partial_t\Pi(t, 0) = \partial_t\Pi(t, T)$, $\dot{\pi}(t) = \partial_t\Pi(t, 0)$ hold, the error equations become

$$\begin{aligned} \dot{\hat{x}}(t) &= -(\nabla f(\tilde{x} + x_m) - \nabla f(x_m)) - A^{\top}\tilde{\lambda}(t, 0) \\ \partial_t \tilde{\lambda}(t, s) &= -\partial_s \tilde{\lambda}(t, s) \\ \tilde{\lambda}(t, 0) &= \tilde{\lambda}(t, T) + \lambda_m(t, T) + (AA^{\top})^{-1}\dot{\pi}(t) + \gamma A\tilde{x}(t) \\ &\quad - \lambda_m(t, T) - (AA^{\top})^{-1}\dot{\pi}(t) \end{aligned}$$
<figure><img src="image_5.png" /><figcaption>Figure 1: Point-wise-delayed solution: time profiles for the first two components of $x_m$ and $\hat{x}_m$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2: Point-wise-delayed solution: time profiles ($t \le 5$ s) for the last three components of $x_m$ and $\hat{x}_m$.</figcaption></figure>

which coincides, within the theoretical framework of this section, with the ones obtained for the algorithms of Section 4. On the other hand, the Lyapunov-like functional

$$V(\tilde{x}, \tilde{\lambda}) = \tilde{x}^{\top} \tilde{x} + \frac{1}{\gamma} \int_{0}^{T} \tilde{\lambda}(\cdot, s)^{\top} \tilde{\lambda}(\cdot, s) ds$$

of Section 5 is nothing but (21) or (23).

### Technical remarks

1. As noted in [2], the algorithm (30) cannot be generically guaranteed to be input-to-state (ISS) with respect to additional perturbations in the boundary conditions (e.g. measurement noise affecting the term $b + \pi(t)$). Nonetheless, the fully saturated scheme (22) can partially mitigate such an issue because it guarantees always boundedness of the trajectories thanks to the saturation function. Finally, the distributed scheme (24) satisfies the aforementioned ISS properties because possible measurement noise affect in this case the dynamics and not the boundary conditions.

2. If a narrow $M_0$ in (18) is not available, the statements of Theorems 2-5 reduce to mere existence statements in terms of a sufficiently large design parameter $M$. In such a case, overly large values of $M$ may increase algorithm sensitivity to noise.

## 7. Simulation results

The simulation results of this section concern the partially saturated algorithm of Section 4.1 and the spectral approximation of Section 5.3 (in its Adaptive Learning Control-like version [39, 37] with no dead-zone action), respectively. Specifically, we consider the case described in Section 2.2, with

<figure><img src="image_7.png" /><figcaption>Figure 3: Distributed-delayed solution: time profiles for the first two components of $x_m$ and $\hat{x}_m$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4: Distributed-delayed solution: time profiles ($t \le 5$ s) for the last three components of $x_m$ and $\hat{x}_m$.</figcaption></figure>

$n = 5$, $m = 4$, $r = 2$, $A_1 = [1, 1; 3, 4]$, $A_2 = [1, 0, 2; 0, 2, 1]$, $b = [2; 3; 1; 2]$, $\pi(t) = [\sin(2\pi t); \cos(2\pi t); 0; 0]$, $f_y(y) = \sum_{j=1}^{2} x_j^2$, $f_z(z) = \sum_{j=3}^{5} x_j^2$, $T = 1$ s. By straightforward computations, $A_1^{-1} = [4, -1; -3, 1]$, $A_2^+ = 1/21[5, -2; -4, 10; 8, 1]$, $I - A_2^+ A_2 = 1/21[16, 4, -8; 4, 1, -2; -8, -2, 4]$ result, while $x_m$, on the basis of (13) [recall that $A_2(I - A_2^+ A_2) = \mathbb{O}$], takes the form: $x_m = [5 + 4\sin(2\pi t) - \cos(2\pi t); -3 - 3\sin(2\pi t) + \cos(2\pi t); 1/21; 16/21; 10/21]$. According to the theoretical developments of Sections 4.1 and 5.3, the two algorithms there presented guarantee $\hat{x}$-convergence to $x_m$. By using $\varphi(t) = t$ over $[0, 1]$ and $\gamma = 30$ (zero initial conditions are set for $\hat{x}$, $M = 300$), the convergence for the former is actually achieved, as Figures 1 and 2 properly illustrate; by setting $\gamma = 60$ (zero initial conditions are set for all the estimates of the Fourier coefficients) and restricting the design to the first 5 components (frequencies: $0, 2\pi/T, 2*(2\pi/T)$) for each Fourier series expansion of the uncertain $T$-periodic vector components, the convergence for the latter is exponentially achieved in light of the technical remark 1 of Section 5 and the Persistency of Excitation arguments [34], as Figures 3 and 4 illustrate. Finally, the components of $\hat{\lambda}(t)$ generated by both the algorithms – not reported for the sake of brevity –, converge to the uncertain $T$-periodic $\lambda_m(t) - (AA^\top)^{-1}\hat{\pi}(t)$ (see [9] for further details). Notice how the different estimation strategy adopted by the two algorithms leads to different estimation behaviors in the transient.

## 8. Conclusions and Future Perspectives

The problem of the online convex optimization over a class of affine periodic constraints (Section 2) has been investigated. By relying on RC techniques, original infinite-dimensional algorithms assessing asymptotic convergence to the optimum have been proposed. Two main routes have been paved, which turn
out to be natural extensions of the approach followed in the con-
stant constraint scenario (Section 3): one resorting to a delayed
equation with point-wise (or boundary) interconnection (Sec-
tion 4), and one resorting to a delayed equation with distributed
interconnection (Section 5). Differences and common struc-
tures have been discussed in detail. A finite-dimensional spec-
tral approximation has been also derived (Section 5), which ex-
hibits deep connections with the adaptive learning control (also
denoted as harmonic regulation) in light of its interpretation
in terms of the internal-model-based design in periodic frame-
works. Future work might concern a possible relaxation of the
technical assumption on the components of the time-varying
constraints in order to handle more general classes of optimiza-
tion problems.

References

[1] K.J. Arrow, L. Hurwicz, and H. B. Chenery. Studies in linear and non-linear programming. 1958.

[2] D. Astolfi. Output-feedback repetitive control for minimum-phase non-linear systems with arbitrarily relative degree. *IFAC-PapersOnLine*, 54(14):464–469, 2021.

[3] D. Astolfi, S. Marx, and N. van de Wouw. Repetitive control design based on forwarding for nonlinear minimum-phase systems. Automatica, 129:109671, 2021.

[4] D. Astolfi, L. Praly, and L. Marconi. Harmonic internal models for struc-
turally robust periodic output regulation. *Systems & Control Letters*, 161:105154, 2022.

[5] D. Astolfi, L. Praly, and L. Marconi. Nonlinear robust periodic output regulation of minimum phase systems. Mathematics of Control, Signals, and Systems, 34(1):129–184, 2022.

[6] M. Bajodek and D. Astolfi. Comparison between different continuous-time realizations based on the tau method of a nonlinear repetitive control scheme. In *62nd IEEE Conference on Decision and Control*, pages 1897–1902, 2023.

[7] N. Bastianello, R. Carli, and S. Zampieri. Internal model-based online optimization. *IEEE Transactions on Automatic Control*, 2023.

[8] G. Bianchin and B. Van Scoy. The internal model principle of time-varying optimization. *arXiv preprint arXiv:2407.08037*, 2024.

[9] S. Bifaretti, P. Tomei, and C.M. Verrelli. Establishing improved convergence and robustness properties for the repetitive learning control. *Applied Mathematics and Computation*, 218:11311–11322, 2012.

[10] M. Bin, J. Huang, A. Isidori, L. Marconi, M Mischiati, and E. Sontag. Internal models in control, bioengineering, and neuroscience. *Annual Review of Control, Robotics, and Autonomous Systems*, 5:55–79, 2022.

[11] M. Bin, I. Notarnicola, and T. Parisini. Semiglobal exponential stability of the discrete-time Arrow-Hurwicz-Uzawa primal-dual algorithm for constrained optimization. *Mathematical Programming*, pages 1–32, 2024.

[12] F. Califano, M. Bin, A. Macchelli, and C. Melchiorri. Stability analysis of nonlinear repetitive control schemes. *IEEE control systems letters*, 2(4):773–778, 2018.

[13] X. Cao and KJ R. Liu. Online convex optimization with time-varying constraints and bandit feedback. *IEEE Transactions on automatic control*, 64(7):2665–2680, 2018.

[14] A. Cecilia, D. Astolfi, M. Bin, and R. Costa-Castelló. Canceling output disturbances in observer design through internal model filters. Automatica, 162:111529, 2024.

[15] V. Cerone, SM. Fosson, S. Pirrera, and D. Regruto. A new framework for constrained optimization via feedback control of lagrange multipliers. *arXiv preprint arXiv:2403.12738*, 2024.

[16] P. Chen, M. Sun, Q. Yan, and X. Fang. Adaptive asymptotic rejection of unmatched general periodic disturbances in output-feedback nonlinear systems. *IEEE Transactions on Automatic Control*, 57:1056–1061, 2012.

[17] A. Cherukuri, E. Mallada, S. Low, and H. Cortés. The role of convexity in saddle-point dynamics: Lyapunov function and robustness. *IEEE Transactions on Automatic Control*, 63(8):2449–2464, 2017.

[18] E. Dall'Anese and A. Simonetto. Optimal power flow pursuit. *IEEE Transactions on Smart Grid*, 9(2):942–952, 2018.

[19] A. Davydov, V. Centorrino, A. Gokhale, G. Russo, and F. Bullo. Contracting dynamics for time-varying convex optimization. *arXiv preprint arXiv:2305.15595*, 2023.

[20] W.E. Dixon, E. Zergeroglu, D.M. Dawson, and B.T. Costic. Repetitive learning control: a Lyapunov-based approach. *IEEE Transactions on Systems, Man, and Cybernetics - PART B: Cybernetics*, 32(4):538–545, 2002.

[21] A. Falsone, I. Notarnicola, G. Notarstefano, and M. Prandini. Tracking-admm for distributed constraint-coupled optimization. Automatica, 117:108962, 2020.

[22] M. Fazlyab, S. Paternain, V.M Preciado, and A. Ribeiro. Prediction-correction interior-point method for time-varying convex optimization. *IEEE Transactions on Automatic Control*, 63(7):1973–1986, 2017.

[23] J. Ghosh and B. Paden. Nonlinear repetitive control. *IEEE Transactions on Automatic Control*, 45(5):949–954, 2000.

[24] M. Giaccagli, D. Astolfi, V. Andrieu, and L. Marconi. Incremental stabilization of cascade nonlinear systems and harmonic regulation: a forwarding-based design. *IEEE Transactions on Automatic Control*, 2024.

[25] A. Hauswirth, Z. He, S. Bolognani, G. Hug, and F. Dörfler. Optimization algorithms as robust feedback controllers. *Annual Reviews in Control*, 57:100941, 2024.

[26] Y. Lin, I. Shames, and D. Nešić. Online convex optimization using coordinate descent algorithms. Automatica, 165:111681, 2024.

[27] R. Marino, P. Tomei, and C.M. Verrelli. Learning control for nonlinear systems in output feedback form. *Systems & Control Letters*, 61:1242–1247, 2012.

[28] A. Simonetto, E. Dall'Anese, S. Paternain, G. Leus, and G.B. Giannakis. Time-varying convex optimization: Time-structured algorithms and applications. *Proceedings of the IEEE*, 108(11):2032–2048, 2020.

[29] A. Simonetto and P. Massioni. Nonlinear optimization filters for stochastic time-varying convex optimization. *International Journal of Robust and Nonlinear Control*, 34(12):8065–8089, 2024.

[30] J.W. Simpson-Porco. Low-gain stability of projected integral control for input-constrained discrete-time nonlinear systems. *IEEE Control Systems Letters*, 6:788–793, 2021.

[31] S. Sun, J. Xu, and W. Ren. Distributed continuous-time algorithms for time-varying constrained convex optimization. *IEEE Transactions on Automatic Control*, 2022.

[32] Y. Tang, G. Qu, and N. Li. Semi-global exponential stability of augmented primal-dual gradient dynamics for constrained convex optimization. Systems & Control Letters, 144:104754, 2020.

[33] P. Tomei and C.M. Verrelli. Linear repetitive learning controls for nonlinear systems by padé' approximants. *International Journal of Adaptive Control and Signal Processing*, 29(6):783–804, 2015.

[34] Patrizio Tomei and Cristiano M. Verrelli. Advances on adaptive learning control: The case of non-minimum phase linear systems. Systems & Control Letters, 115(2):55–62, 2018.

[35] B. Van Scoy, J.W. Simpson-Porco, and L. Lessard. Automated Lyapunov analysis of primal-dual optimization algorithms: An interpolation approach. In *62nd IEEE Conference on Decision and Control*, pages 1306–1311, 2023.

[36] C.M. Verrelli. A larger family of nonlinear systems for the repetitive learning control. Automatica, 71:38–43, 2016.

[37] C.M. Verrelli. PI-generalizing saturated repetitive learning control for a class of nonlinear uncertain systems: Robustness w.r.t relative degree zero or one. *Systems & Control Letters*, 164(105248), 2022.

[38] C.M. Verrelli. Exponential output tracking via a saturated rlc for a class of nonlinear systems under psd of periodic references. *Systems & Control Letters*, 192:10588, 2024.

[39] C.M. Verrelli and P. Tomei. Adaptive learning control for nonlinear systems: A single learning estimation scheme is enough. Automatica, 149:110833, 2023.

[40] Y. Zhang, E. Dall'Anese, and M. Hong. Online proximal-admm for time-varying constrained convex optimization. *IEEE Transactions on Signal and Information Processing over Networks*, 7:144–155, 2021.

[41] T. Zheng, J.W. Simpson-Porco, and E. Mallada. Closed-loop motion planning for differentially flat systems: A time-varying optimization framework. *arXiv preprint arXiv:2310.13090*, 2023.