Prediction-Correction Splittings for Nonsmooth Time-Varying Optimization

TL;DR

This paper introduces a prediction-correction splitting algorithm for solving nonsmooth, time-varying convex optimization problems, with proven convergence and demonstrated effectiveness in robotic formation control.

Contribution

It develops a novel online prediction-correction framework using splitting methods for nonsmooth, time-varying optimization, with convergence guarantees.

Findings

Convergence to a neighborhood of the optimal solution trajectory.

Effective performance demonstrated in robotic leader-following simulations.

Applicable to a broad class of nonsmooth, time-varying problems.

Abstract

We address the solution of time-varying optimization problems characterized by the sum of a time-varying strongly convex function and a time-invariant nonsmooth convex function. We design an online algorithmic framework based on prediction-correction, which employs splitting methods to solve the sampled instances of the time-varying problem. We describe the prediction-correction scheme and two splitting methods, the forward-backward and the Douglas-Rachford. Then by using a result for generalized equations, we prove convergence of the generated sequence of approximate optimizers to a neighborhood of the optimal solution trajectory. Simulation results for a leader following formation in robotics assess the performance of the proposed algorithm.

Figures

Equations

$$\text{mathbold}{x}^{\ast }(t)=\underset{\text{mathbold}x\in {\mathbb{R}}^n}{\mathrm{argmin}}\{f(\text{mathbold}x;t)+g(\text{mathbold}x)\}$$

$$\text{mathbold}{x}^{\ast }(t_k)=\underset{\text{mathbold}x\in {\mathbb{R}}^n}{\mathrm{argmin}}\{f(\text{mathbold}x;t_k)+g(\text{mathbold}x)\},k\in \mathbb{N}.$$

$$\text{mathbold}x\mapsto \left\{\text{mathbold}z\in {\mathbb{R}}^n|\forall \text{mathbold}y\in {\mathbb{R}}^n:⟨\text{mathbold}y-\text{mathbold}x,\text{mathbold}z⟩+\phi (\text{mathbold}x)\le \phi (\text{mathbold}y)\right\}.$$

$$\Vert {\nabla }_{\text{mathbold}x\text{mathbold}x\text{mathbold}x}f(\text{mathbold}x;t)\Vert$$

$$\Vert {\nabla }_{tt\text{mathbold}x}f(\text{mathbold}x;t)\Vert \le C_3.$$

$${\nabla }_{\text{mathbold}x}f(\text{mathbold}{x}^{\ast }(t_{k+1});t_{k+1})+\partial g(\text{mathbold}{x}^{\ast }(t_{k+1}))\ni 0.$$

$$\begin{array}{rl}\nabla h_k(\text{mathbold}x)& ={\nabla }_{\text{mathbold}x}f(\text{mathbold}{x}_k;t_k)+\\ & +{\nabla }_{\text{mathbold}x\text{mathbold}x}f(\text{mathbold}{x}_k;t_k)(\text{mathbold}x-\text{mathbold}{x}_k)+T_{\mathrm{s}}{\nabla }_{t\text{mathbold}x}f(\text{mathbold}{x}_k;t_k)\end{array}$$

$$\nabla h_k({\text{mathbold}x}_{k+1|k})+\partial g({\text{mathbold}x}_{k+1|k})\ni 0.$$

$${\nabla }_{t\text{mathbold}x}f(\text{mathbold}{x}_k;t_k)\simeq \left({\nabla }_{\text{mathbold}x}f(\text{mathbold}{x}_k;t_k)-{\nabla }_{\text{mathbold}x}f(\text{mathbold}{x}_k;t_{k-1})\right)/T_{\mathrm{s}},$$

$${\mathrm{prox}}_{\rho f}(\text{mathbold}x)=\underset{\text{mathbold}y\in {\mathbb{R}}^n}{\mathrm{argmin}}\left\{f(\text{mathbold}y)+\frac{1}{2\rho }{\Vert \text{mathbold}y-\text{mathbold}x\Vert }^2\right\},$$

$${\mathcal{T}}_{\mathrm{FB}}={\mathrm{prox}}_{\rho \gamma }\circ (I-\rho \nabla \phi )$$

$${\mathcal{T}}_{\mathrm{DR}}=\frac{1}{2}\left(I+{\mathrm{refl}}_{\rho \gamma }\circ {\mathrm{refl}}_{\rho \phi }\right)$$

$$\begin{array}{rl}& \text{mathbold}y(k)=\text{mathbold}x(k)-\rho \nabla \phi (\text{mathbold}x(k))\\ & \text{mathbold}x(k+1)={\mathrm{prox}}_{\rho \gamma }(\text{mathbold}y(k)).\end{array}$$

$$\begin{array}{rl}& \text{mathbold}x(k)={\mathrm{prox}}_{\rho \phi }(\text{mathbold}z(k))\\ & \text{mathbold}y(k)={\mathrm{prox}}_{\rho \gamma }(2\text{mathbold}x(k)-\text{mathbold}z(k))\\ & \text{mathbold}z(k+1)=\text{mathbold}z(k)+\text{mathbold}y(k)-\text{mathbold}x(k)\end{array}$$

$$\Vert \text{mathbold}x(k+1)-\text{mathbold}{x}^{\ast }\Vert \le {\zeta }_{\text{FB}}\Vert \text{mathbold}x(k)-\text{mathbold}{x}^{\ast }\Vert .$$

$${\zeta }_{\text{FB}}=\max \{|1-\rho m|,|1-\rho L|\},$$

$$\Vert \text{mathbold}x(k)-\text{mathbold}{x}^{\ast }\Vert \le {\zeta }_{\text{DR}}^k\frac{1+\rho L}{1+\rho m}\Vert \text{mathbold}x(0)-\text{mathbold}{x}^{\ast }\Vert ,$$

$${\zeta}_{\textrm{DR}}=\max\Big{\{}\frac{1}{1+\rho m},\frac{\rho L}{1+\rho L}\Big{\}}.$$

$$\nabla h_k({\text{mathbold}x}_{k+1|k})+\partial g({\text{mathbold}x}_{k+1|k})\ni 0.$$

$${\nabla }_{\text{mathbold}x}f(\text{mathbold}x;t_{k+1})+\partial g(\text{mathbold}x)\ni 0.$$

$$\zeta (C)\left[\zeta (P)+(\zeta (P)+1)\frac{2L}{m}\right]<1.$$

$$\underset{k\to \infty }{\mathrm{limsup}}\Vert \text{mathbold}{x}_k-\text{mathbold}{x}_k^{\ast }\Vert =O(\zeta (C)T_{\mathrm{s}}),$$

$$\underset{k\to \infty }{\mathrm{limsup}}\Vert \text{mathbold}{x}_k-\text{mathbold}{x}_k^{\ast }\Vert =O(\zeta (C)T_{\mathrm{s}}^2)+O(\zeta (P+C)T_{\mathrm{s}}).$$

$${\bar{T}}_{\mathrm{s}}=\frac{\tau -\zeta (P+C)}{\zeta (C)(\zeta (P)+1)}{\left(\frac{C_0{C}_1}{m^2}+\frac{C_2}{m}\right)}^{-1}$$

$$\bar{R}=\frac{2m}{C_1}\left(\frac{C_0{C}_1}{m^2}+\frac{C_2}{m}\right)({\bar{T}}_{\mathrm{s}}-T_{\mathrm{s}}).$$

$$\nabla \phi (\text{mathbold}x)+\partial \gamma (\text{mathbold}x)\ni \text{mathbold}p$$

$$S(\text{mathbold}p)=\{\text{mathbold}y|\nabla \phi (\text{mathbold}y)+\partial \gamma (\text{mathbold}y)\ni \text{mathbold}p\}.$$

$$\text{mathbold}{y}^{\ast }=\underset{\text{mathbold}}{\mathrm{argmin}}y\{\phi (\text{mathbold}y)+\gamma (\text{mathbold}y)-⟨\text{mathbold}p,\text{mathbold}y⟩\}.$$

$$⟨\text{mathbold}z-\text{mathbold}y,\text{mathbold}p-\nabla \phi (\text{mathbold}y)⟩$$

$$⟨\text{mathbold}z-\text{mathbold}w,\mathbf{q}-\nabla \phi (\text{mathbold}w)⟩$$

$$⟨\text{mathbold}w-\text{mathbold}y,\text{mathbold}p-\nabla \phi (\text{mathbold}y)⟩\le \gamma (\text{mathbold}w)-\gamma (\text{mathbold}y)\le ⟨\text{mathbold}w-\text{mathbold}y,\mathbf{q}-\nabla \phi (\text{mathbold}w)⟩$$

Full text

I Introduction

We are interested in the solution of time-varying optimization problems in the form

[TABLE]

where $f:\mathbb{R}^{n}\times\mathbb{R}_{+}\to\mathbb{R}$ is smooth and strongly convex, and $g:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is proper, closed and convex, but possibly non-differentiable. Since the solution $\mathbold{x}^{*}(t)$ – the trajectory – changes over time, our objective is to “track” such trajectory up to an error bound.

Problems conforming to the framework of (1) arise in many applications. The powerful and widely used model predictive control (MPC) requires that we solve an optimization problem that varies with time [1, 2, 3] in order to compute the control action. In signal processing, the estimation of time-varying signals on the basis of observations gathered online can be cast as the problem of solving a series of varying problems [4, 5, 6]. In robotics, path tracking and leader following problems can be cast in the framework of (1), see for example [7, 8, 9]. Other application domains are economics [10], smart grids [11], and non-linear optimization [12, 13].

In this paper, we are interested in solving problem (1) in a discrete-time framework, which is amenable to direct implementation with digital hardware. In particular, we discretize the problem with a sampling period $T_{\mathrm{s}}:=t_{k+1}-t_{k}$, which results in the sequence of time-invariant problems

[TABLE]

Due to the time-varying nature of the problem, the sequence of problems is not known ahead of time, and each problem $k$ is revealed when sampled at $t_{k}$.

The smaller the sampling time, the higher the accuracy of the solution computed will be. However, the time-invariant problems may be too difficult to solve inside a sampling period, hence there exists a natural trade-off between accuracy and practical implementation constraints.

We are interested in the solution of (2) using a prediction-correction scheme. Prediction-correction algorithms have been studied both in the discrete-time framework that we employ [14, 15, 16] and in a continuous-time setup [17, 18, 19].

In the cited works, however, the proposed algorithms can handle only smooth optimization problems. Our goal in this paper is to extend the proposed approach to the more general nonsmooth (2). In particular, in [15] $g(\mathbold{x})$ could be interpreted as the indicator function of a convex set, while in this work we are interested in accepting general, nonsmooth convex functions.

While the prediction-correction scheme will be described in detail in the following section, we mention here that two optimization problems will be needed to be solved at each time $t_{k}$. The solution of these problems will be obtained by using splitting methods, and in particular the forward-backward and the Douglas-Rachford splitting. For background and review of these methods we refer the interested reader to [20, 21], and for some examples of application to convex optimization problems [22, 23].

Organization

The paper is structured as follows. Section II introduces the proposed prediction-correction scheme with splitting methods. Section III states the convergence results for the proposed algorithm, and reports a theoretical finding instrumental in proving them. Section IV describes the numerical simulations carried out to evaluate the performance of the algorithm. Section V concludes the paper. Proofs can be found in the Appendix.

Basic definitions

We say that a function $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ is $m$-strongly convex for a constant $m\in\mathbb{R}_{+}$ iff $\varphi(\mathbold{x})-\frac{m}{2}\|\mathbold{x}\|^{2}$ is convex. The function $\varphi$ is said to be $L$-smooth if its gradient is $L$-Lipschitz continuous, or equivalently $\varphi(\mathbold{x})-\frac{L}{2}\|\mathbold{x}\|^{2}$ is concave. We denote the class of $m$-strongly convex and $L$-smooth functions with $\mathcal{S}_{m,L}(\mathbb{R}^{n})$.

A function $\varphi:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is said to be closed if for any $a\in\mathbb{R}$ the set $\{\mathbold{x}\in\operatorname{dom}(f)\ |\ \varphi(\mathbold{x})\leq a\}$ is closed. A function is said to be proper if it does not attain $-\infty$. We denote the class of closed, convex, and proper functions with $\Gamma_{0}(\mathbb{R}^{n})$.

Given $\varphi\in\Gamma_{0}(\mathbb{R}^{n})$ we define its subdifferential as the set-valued operator $\partial\varphi:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ such that

[TABLE]

II Prediction-Correction and Splittings

We describe the prediction-correction scheme that will be applied to solve (2), and we introduce the splitting methods that will be employed in combination with it. We will use the following assumptions.

Assumption 1

The function $f:\mathbb{R}^{n}\times\mathbb{R}_{+}\to\mathbb{R}$ belongs to $\mathcal{S}_{m,L}(\mathbb{R}^{n})$ uniformly in time. The function $g:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ belongs to $\Gamma_{0}(\mathbb{R}^{n})$ and it is in general nonsmooth.

Assumption 2

The function $f$ has bounded time derivative of its gradient derivative as: $\left\lVert\nabla_{t\mathbold{x}}f(\mathbold{x};t)\right\rVert\leq C_{0}$.

Assumption 3

The function $f$ is at least three times differentiable and has bounded derivatives w.r.t. $\mathbold{x}\in\mathbb{R}^{n}$ and $t\in\mathbb{R}_{+}$ as:

[TABLE]

Assumption 1 ensures by strong convexity that the solution to the problem is unique at each time, and that the gradient of $f$ is Lipschitz continuous. These conditions are common for time-varying scenarios, e.g., see [14, 10]. Assumption 2 further guarantees that the temporal variability of the gradient of $f$ is bounded, and thus that the problem does not vary too widely between time instants. Finally, Assumption 3 – which is not strictly necessary to prove convergence of our algorithm – imposes boundedness of the tensor $\nabla_{\mathbold{x}\mathbold{x}\mathbold{x}}f(\mathbold{x};t)$, which is a typical assumption for second-order algorithms. Moreover, it bounds the variability of the Hessian of $f$ over time, which guarantees the possibility of performing more accurate predictions of the optimal trajectory.

II-A Prediction-correction

The scheme is characterized by two phases, the prediction step, which at time $t_{k}$ computes an approximate solution to the problem at time $t_{k+1}$; and the correction step, during which the predicted solution is refined using the information that becomes available at time $t_{k+1}$.

The problem that we are interested to approximately solve during the prediction step is $\mathbold{x}^{*}(t_{k+1})=\operatorname*{arg\,min}_{\mathbold}{x}\{f(\mathbold{x};t_{k+1})+g(\mathbold{x})\}$, which is equivalently characterized by the generalized equation

[TABLE]

However, at time $t_{k}$ the function $f(\mathbold{x};t_{k+1})$ is not known and therefore we need to approximate it using the information available, i.e., $f(\mathbold{x};t_{k})$ and the current state $\mathbold{x}(t_{k})=:\mathbold{x}_{k}$. This setup models for instance problems in which the cost function depends on observations that are gathered online.

In order to compute the prediction – that we will denote $\tilde{\mathbold{x}}_{k+1|k}$ – we approximate the gradient of $f(\mathbold{x};t_{k+1})$ by using a Taylor expansion around $(\mathbold{x}_{k};t_{k})$, which is given by

[TABLE]

where the subscript indicates that $h$ is computed with information available at time $t_{k}$.

Remark 1

By Assumption 1, it follows that $h_{k}$ belongs to $\mathcal{S}_{m,L}(\mathbb{R}^{n})$, since $\nabla^{2}h_{k}(\mathbold{x})=\nabla_{\mathbold{x}\mathbold{x}}f(\mathbold{x}_{k};t_{k})$.

An approximate solution to (3) can now be computed solving the new generalized equation

[TABLE]

Thus the prediction step requires the solution of this approximated generalized equation with initial condition $\mathbold{x}_{k}$.

At time $t_{k+1}$ the new information $f(\mathbold{x};t_{k+1})$ is made available and therefore the correction step can be computed. The correction step approximately solves (3) by using the prediction as initial condition.

Remark 2

In order to compute the approximation $\nabla h_{k}(\mathbold{x})$ of the gradient of $f(\mathbold{x};t_{k+1})$ we assumed to know the derivative $\nabla_{t\mathbold{x}}f(\mathbold{x};t_{k})$. However in practice this is not always possible, and thus we can think of approximating the derivative using a first-order backward finite difference

[TABLE]

which introduces an error bounded by $O(T_{\mathrm{s}}^{2})$ [14].

II-B Splitting methods

We have defined two problems that need solving, represented by the generalized equations (3) and (5). In order to do so, we can apply splitting methods, that are well suited to minimizing the sum of two (possibly non-differentiable) convex functions. In the following we introduce some necessary background on operator theory and the two most widely used splitting methods, the forward-backward and the Douglas-Rachford. For a treatment of operator theory and its uses for convex optimization we refer the reader to [20, 21].

Definition 1

A mapping $\mathcal{T}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is said to be nonexpansive if it has unitary Lipschitz constant, that is if $\left\lVert\mathcal{T}\mathbold{x}-\mathcal{T}\mathbold{y}\right\rVert\leq\left\lVert\mathbold{x}-\mathbold{y}\right\rVert$ for any two $\mathbold{x},\mathbold{y}\in\mathbb{R}^{n}$.

Definition 2

Let $f:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ be a closed, proper and convex function, and let $\rho>0$. The corresponding proximal operator is defined as

[TABLE]

and the reflective operator is $\operatorname{refl}_{\rho f}(\mathbold{x})=2\operatorname{prox}_{\rho f}(\mathbold{x})-\mathbold{x}$.

Definition 3

Let $\mathcal{T}:\mathbb{R}^{n}\to\mathbb{R}^{n}$, the fixed points of this operator are all the points $\bar{\mathbold{x}}\in\mathbb{R}^{n}$ such that $\bar{\mathbold{x}}=\mathcal{T}\bar{\mathbold{x}}$.

Operator theory can be employed to solve convex optimization problems by designing suitable operators whose fixed points are the solutions of the original problem. Consider in particular the problem $\min_{\mathbold}{x}\{\varphi(\mathbold{x})+\gamma(\mathbold{x})\}$, where $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ belongs to $\mathcal{S}_{m,L}(\mathbb{R}^{n})$, while $\gamma:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ belongs to $\Gamma_{0}(\mathbb{R}^{n})$. This is equivalent to asking that we solve the generalized equation $\nabla\varphi(\mathbold{x})+\partial\gamma(\mathbold{x})\ni 0$111Where recall that $\gamma(\mathbold{x})$ is not necessarily smooth.. Notice that both (3) and (5) conform to this class of problems.

Since the problem of interest has a separable structure, we introduce now two examples of splitting operators that exploit this fact. The first is the forward-backward (FB) operator, defined as

[TABLE]

where $\rho<2/L$, and the second is the Douglas-Rachford (DR)

[TABLE]

where it must be $\rho>0$.

In order to derive the fixed points of the operators above we can employ the Banach-Picard iteration, which consists of the repeated application of the chosen operator [21]. The resulting algorithms are called splitting methods, and in particular the forward-backward spitting (FBS) is $\mathbold{x}(k+1)=\mathcal{T}_{\mathrm{FB}}\mathbold{x}(k)$, $k\in\mathbb{N}$, or equivalently

[TABLE]

The fixed points of $\mathcal{T}_{FB}$ computed via FBS coincide with the solutions to the original problem $\min_{\mathbold}{x}\{\varphi(\mathbold{x})+\gamma(\mathbold{x})\}$.

The Douglas-Rachford splitting (DRS) is instead applied to the auxiliary variable $\mathbold{z}$ and not – as is the FBS – directly to $\mathbold{x}$. The DRS is thus described by the update $\mathbold{z}(k+1)=\mathcal{T}_{\mathrm{DR}}\mathbold{z}(k)$, $k\in\mathbb{N}$, or equivalently

[TABLE]

where the additional variable $\mathbold{y}$ is introduce to break the computations in smaller pieces. Let now $\mathbold{z}^{*}$ be a fixed point of $\mathcal{T}_{DR}$ computed by the DRS, the corresponding solution to the original problem can be computed as $\mathbold{x}^{*}=\operatorname{prox}_{\rho\varphi}(\mathbold{z}^{*})$.

Under the assumptions above, for the $\varphi$ and $\gamma$ functions, it is possible to prove that the two splitting methods converge to a fixed point of the respective operator, and thus that they reach a solution to the original problem.

In particular, for the FBS it is possible to show that [24]

[TABLE]

where

[TABLE]

and recursively, $\left\lVert\mathbold{x}(k)-\mathbold{x}^{*}\right\rVert\leq\zeta_{\textrm{FB}}^{k}\left\lVert\mathbold{x}(0)-\mathbold{x}^{*}\right\rVert$

For the DRS, instead, the following Lemma holds.

Lemma 1

The sequence $\{\mathbold{x}(k)\}_{k\in\mathbb{N}}$ generated by the DRS according to (9) satisfies the following inequality

[TABLE]

[TABLE]

$\square$

Proof:

Given in the Appendix. ∎

In the following, we will use $\zeta(k)$ as a function of $k$ to indicate either $\zeta^{k}_{\textrm{FB}}$ or $\zeta^{k}_{\textrm{DR}}\frac{1+\rho L}{1+\rho m}$.

II-C Solving generalized equations

As mentioned above, we are now interested in solving the prediction and correction steps, i.e., the generalized equations (5) and (3). That is, we need to minimize $h_{k}(\mathbold{x})+g(x)$ for the prediction step, and $f(\mathbold{x};t_{k+1})+g(\mathbold{x})$ for the correction step. For this, we can apply the operator splitting methods described in the previous part of this section.

However, the splitting methods are guaranteed to converge asymptotically to the desired solution, but for practical reasons (finite length of the sampling period) we stop the algorithms after a fixed number of steps: $P$ for the prediction and $C$ for the correction phases. Therefore we obtain in both steps an approximation of the optimal solution.

Using the convergence results reported above, it is therefore possible to characterize the convergence rate of both splitting methods and quantify the reduction in the error after the prescribed number of steps.

We have now all the ingredients to fully describe the proposed prediction-correction algorithm, which is reported in Algorithm 1. At every time $t_{k}$, we perform $P$ step of either the FBS or the DRS on the predicting generalized equation,

[TABLE]

constructed for the prediction problem (5) [cf. line 4]; this yields an approximate predictor $\tilde{\mathbold{x}}_{k+1|k}$.

At time $t_{k+1}$, we observe the new cost function $f(\cdot;t_{k+1})$ [cf. line 6], and we perform $C$ steps of either the FBS or the DRS on the correcting generalized equation,

[TABLE]

constructed for the correction problem (3) [cf. line 7]; this yields the approximate optimizer $\mathbold{x}_{k+1}$.

Remark 3

We choose to use the same step-sizes for the prediction and correction steps, since both $f(\mathbold{x};t_{k+1})$ and $h_{k}(\mathbold{x})$ belong to $\mathcal{S}_{m,L}(\mathbb{R}^{n})$ and thus the step-sizes can be chosen minimizing (11) and (13).

III Convergence Analysis

Algorithm 1 entails the approximated solution of two optimization problems, therefore we do not expect the algorithm to converge exactly to the optimal trajectory. Instead, this section aims at proving the convergence of Algorithm 1 to a neighborhood of $\mathbold{x}^{*}(t_{k})=:\mathbold{x}_{k}^{*}$ which depends on the sampling period $T_{\mathrm{s}}$.

We state the following results, which will be proved in the remainder of this Section and in the Appendix.

Theorem 1

Let Assumptions 1-2 hold, and choose the parameters of Algorithm 1 such that

[TABLE]

Then the trajectory $\{\mathbold{x}_{k}\}_{k\in\mathbb{N}}$ generated by Algorithm 1 converges Q-linearly to a neighborhood of the optimal trajectory $\{\mathbold{x}_{k}^{*}\}_{k\in\mathbb{N}}$, and in particular

[TABLE]

where function $\zeta(k)$ is either $\zeta^{k}_{\textrm{FB}}$ or $\zeta^{k}_{\textrm{DR}}\frac{1+\rho L}{1+\rho m}$, depending if one uses FBS or DRS in Algorithm 1. $\square$

Theorem 2

Let Assumptions 1, 2, and 3 hold, and select the prediction and correction horizons $P$ and $C$, and $\tau\in(0,1)$, such that $\zeta({P+C})<\tau$.

Then there exist an upper bound for the sampling time $\bar{T}_{\mathrm{s}}$ and a convergence region $\bar{R}$ such that if $T_{\mathrm{s}}<\bar{T}_{\mathrm{s}}$ and $\left\lVert\mathbold{x}_{0}-\mathbold{x}_{0}^{*}\right\rVert<\bar{R}$, then

[TABLE]

In particular, the bound for the sampling time and the convergence region are characterized by

[TABLE]

and

[TABLE]

*Function $\zeta(k)$ is either $\zeta^{k}_{\textrm{FB}}$ or $\zeta^{k}_{\textrm{DR}}\frac{1+\rho L}{1+\rho m}$, depending if one uses FBS or DRS in Algorithm 1. *$\square$

Remark 4

Notice that in Theorem 1 the asymptotical error bound depends linearly on the sampling time; this result agrees with the fact that the smaller $T_{\mathrm{s}}$ is, the better the algorithm can track the optimal trajectory.

Moreover, if bounds on the third order derivatives of $f$ are available, then we can characterize the convergence to a smaller neighborhood of the optimal trajectory. Indeed Theorem 2 shows that the asymptotical error depends on the square of $T_{\mathrm{s}}$ and thus is tighter than the result of Theorem 1.

Remark 5

Theorem 1 holds globally, while Theorem 2 holds only if the initial condition $\mathbold{x}_{0}$ is selected in a neighborhood $\bar{R}$ of the optimum. Moreover, notice that the larger we choose $T_{\mathrm{s}}$, the smaller $\bar{R}$ becomes, again in accordance with the fact that the closer we sample the problem the better we approximate the optimal trajectory.

III-A Proof of Theorems 1 and 2

The full proof of Theorems 1 and 2 follows closely along the lines of [15, Appendix B], and it is reported in the Appendix for completeness. Here however, we give a sketch of the reasoning and a result on implicit solution mappings, which enables us to generalize [25] to our setting.

Recall from Algorithm 1 that the following sources of error are present: (i) the Taylor approximation used during the prediction step, and (ii-iii) the early termination of the splitting method during both the prediction and correction steps. Thus, in order to evaluate the overall error committed by the proposed algorithm, i.e., $\left\lVert\mathbold{x}_{k}-\mathbold{x}_{k}^{*}\right\rVert$, we need to derive a bound on the errors coming from (i), (ii), and (iii).

The proof proceeds as follows. We bound the Taylor approximation used during the prediction step by using a result in implicit solution mappings (see Theorem 3). We then bound the error coming from the early termination of the splitting method during prediction, by using the linear convergence rate of the splitting method (cf. Eq.s (10) and (12)). The next step is the derivation of a bound on the early termination of the splitting method during correction, which is analogous to the one during prediction. Finally, we put everything together and derive an error recursion, from which the claims.

Since the main novelty of the proof is the mentioned result in implicit solution mappings, we present this here. First of all, consider the parametrized generalized equation222In which the differentiation is taken w.r.t. to the $\mathbold{x}$.

[TABLE]

where $\mathbold{p}\in\mathbb{R}^{d}$ is a parameter, and define the solution mapping of (15) as

[TABLE]

We can characterize the Lipschitz continuity of the solution mapping (16) using the following Theorem 3333This result extends Theorem 2F.6 in [25], which holds only if $\gamma$ is the indicator function of a convex set, and thus $\partial\gamma$ is a normal cone – that is, it holds only for variational inequalities and not generalized equations. Theorem 3 here resembles [26, Th. 1], although the latter is applied to a slightly different setting. .

Theorem 3

Let $\varphi\in\mathcal{S}_{m,L}(\mathbb{R}^{n})$ and $\gamma\in\Gamma_{0}(\mathbb{R}^{n})$, then the solution mapping (16) of the parameterized equation (15) is single-valued and $m^{-1}$-Lipschitz continuous. $\square$

Proof:

The proof is divided in two parts.

Single-valuedness

The solutions of (15) are the solutions of

[TABLE]

But since $\varphi$ is strongly convex, so is the whole objective function once the parameter $\mathbold{p}$ is fixed, which implies that the problem has a unique solution for each $\mathbold{p}$.

Lipschitz continuity

Let $\mathbold{y}=S(\mathbold{p})$ and $\mathbold{w}=S(\mathbf{q})$ for $\mathbold{p},\mathbf{q}\in\mathbb{R}^{n}$. In order to prove Lipschitz continuity we need to show that $\left\lVert\mathbold{y}-\mathbold{w}\right\rVert\leq m^{-1}\left\lVert\mathbold{p}-\mathbf{q}\right\rVert$. By the definition of solution mapping it holds $\mathbold{p}-\nabla\varphi(\mathbold{y})\in\partial\gamma(\mathbold{y})$ and $\mathbf{q}-\nabla\varphi(\mathbold{w})\in\partial\gamma(\mathbold{w})$; therefore by the definition of the subdifferential it must hold

[TABLE]

for any $\mathbold{z}\in\mathbb{R}^{n}$. In particular, inequality (17) is satisfied for $\mathbold{z}=\mathbold{w}$, and (18) for $\mathbold{z}=\mathbold{y}$; changing sign in the latter, the following inequalities hold

[TABLE]

or, equivalently,

[TABLE]

Since $\nabla\varphi$ is the derivative of an $m$-strongly convex function, then it is strongly monotone [21], in the sense that for any $\mathbold{x},\mathbold{y}\in\mathbb{R}^{n}$ it holds

[TABLE]

With this result and the fact that $\langle\mathbold{w}-\mathbold{y},\mathbf{q}-\mathbold{p}\rangle\leq\left\lVert\mathbold{w}-\mathbold{y}\right\rVert\left\lVert\mathbf{q}-\mathbold{p}\right\rVert$, it follows then

[TABLE]

which proves the desired result. ∎

IV Simulations

To assess the practical performance of the proposed methods, we present some simple yet realistic numerical simulations for a leader following case problem.

IV-A Problem formulation

We are interested in solving a leader following problem for a group of $N$ robots in a rigid formation. In particular, the leader (or [math]-th robot) follows freely a trajectory $\{\mathbold{x}_{k}^{0}\}_{k\in\mathbb{N}}$ on the plane, and the follower robots (labeled $1$ through $N$) must: (i) follow the leader, while (ii) keeping a fixed formation. Hereafter, we denote with $\mathbold{x}^{i}\in\mathbb{R}^{2}$ the position of the $i$-th robot, and define $\mathbold{x}=[(\mathbold{x}^{0})^{\top},(\mathbold{x}^{1})^{\top},\ldots,(\mathbold{x}^{N})^{\top}]^{\top}\in\mathbb{R}^{2(N+1)}$.

The followers, however, have access only to a noisy and partial measurement of the leader’s position, which is given by $z_{k+1}^{i}=\mathbf{v}_{i}^{\top}\mathbold{x}_{k+1}^{0}+n^{i}_{k+1}$ where $\mathbf{v}_{i}\in\{[1,0]^{\top},[0,1]^{\top}\}$ and the noise is Gaussian $n^{i}\sim\mathcal{N}(0,\sigma_{i})$.

Inspired by the formulation presented in [9], the idea is then to estimate $\mathbold{x}^{0}$ by solving a least squares problem as

[TABLE]

The formation is to be kept constant, with the leader robot in the center and the followers at a distance $d$ from it, with equal angles between them. By defining a reference frame centered in the leader, with $x$-axis pointing in the direction of motion, the position of the $i$-th robot must satisfy

[TABLE]

These constraints are collected in the linear system $\mathbold{A}\mathbold{x}=\mathbold{b}$:

[TABLE]

and imposed by choosing $g(\mathbold{x})=\iota_{\mathbold{A}\mathbold{x}=\mathbold{b}}(\mathbold{x})$.

Notice that the problem is solved by a fusion center that collects the measurements and then dictates the trajectory that the followers need to perform.

IV-B Numerical results

The simulations where carried out for a formation of $N=10$ followers, placed on the unitary circle around the leader ($d=1$). The noise variance was equal to $\sigma_{i}=0.1$ for all robots, and six followers had $v_{i}=[1,0]^{\top}$ while the other four had $v_{i}=[0,1]^{\top}$. The regularization constant was $\lambda=10$, for the FBS we set $\rho=1/L$, while for DRS we set $\rho=0.08$ (both to maximize the respective contraction ratio). The sampling time was set to $T_{\mathrm{s}}=0.1$\mathrm{s}. The trajectory of the leader is a Lissajous curve with ratio 1:3 and spans an area of $[-3,3]^{2}$ with a period of about $40$\mathrm{s}. Our performance metric is the tracking error, defined as $E_{k}=\left\lVert\mathbold{x}_{k}-\mathbold{x}_{k}^{*}\right\rVert$.

First of all, we evaluate the performance of FBS vs. DRS with $P=0,C=5$ (Figure 2-top); as we see, FBS outperforms DRS for this problem setting (but this is far from a general statement).

Second, we look at FBS for different values of the prediction horizon $P$, and with fixed correction horizon $C=5$ (Figure 2-middle). The presence of at least one prediction step improves the accuracy of the algorithm, thus justifying the choice of a prediction-correction scheme. We see also that, since we compute the prediction by using a backward derivative approximation of the time derivative of the gradient, and noise is added, the case $P=5$ hits the noise “floor”.

Figure 2-bottom depicts instead the error for different values of the correction steps for FBS, with $P=0$.

The smaller the number of correction steps is, the worse the asymptotical error. This result shows that augmenting $C$ seems to be more effective in reducing the tracking error than augmenting $P$, in this particular case. Note that however, correction steps are more computational expensive than prediction steps (since the latter are computed on a quadratic cost function).

We then evaluate the effect of the sampling time on the performance of the proposed algorithm by computing the asymptotical error for different values of $T_{\mathrm{s}}$, $P$ and $C$, which is computed as $\max_{k>2D/3}E_{k}$ with $D=100$\mathrm{s}$$ the duration of the simulation. In this case (and for verification purposes only) we let the noise be depended on the sampling period as $\sigma_{i}=.01T_{\mathrm{s}}$, so to remove its effect from the numerical results. (Naturally this is not the case in practice and $T_{\mathrm{s}}$ needs to be chosen to limit the effect of the noise on the prediction, as well as maintain an acceptable update frequency).

In Figure 3 shows that the larger the sampling time is, the worse the performance of the proposed algorithm is, since the sequence of problems (2) “tracks” the original problem worse. Notice that, as expected from the theorems, we obtain the $O(T_{\mathrm{s}})$ and $O(T_{\mathrm{s}}^{2})$ dependencies for the asymptotical error, further justifying the use of prediction steps.

V Conclusions and Future Directions

In this work, we presented an algorithm to solve time-varying optimization problems that applies splitting operators in conjunction with a prediction-correction scheme. We described two theoretical results that guarantee the convergence of the sequence computed by the algorithm to a neighborhood of the optimal solutions that depends on the sampling time, or its square. Finally we described some numerical results obtained by applying the algorithm to a leader following problem in robotics.

Future work will focus on extending the algorithm to solve problems characterized by weaker conditions, such as removing strong convexity. Moreover, the application to distributed optimization problems will be explored, as well as further acceleration via e.g., [27].

-A Taylor approximation error

By making use of Theorem 3, we derive an upper bound to the error introduced by the Taylor approximation. Consider the generalized equation (5) that characterizes the prediction step, and define the following functions:

[TABLE]

and

[TABLE]

With these definitions the prediction problem (5) is equivalent to solving $(\psi+\Psi)(\mathbold{x}_{k+1|k})\ni 0$, while the problem (3) is equivalent to $\Psi(\mathbold{x}_{k+1}^{*})\ni 0$, and therefore it holds $(\psi+\Psi)(\mathbold{x}_{k+1}^{*})\ni\psi(\mathbold{x}_{k+1}^{*})$.

Consider now the parametrized generalized equation

[TABLE]

with solution mapping $S(\mathbold{p})$. By the results above we have $S(\mathbold{p})=\mathbold{x}_{k+1|k}$ when $\mathbold{p}=0$, and $S(\mathbf{q})=\mathbold{x}_{k+1}^{*}$ when $\mathbf{q}=\psi(\mathbold{x}_{k+1}^{*})$. Employing now Theorem 3 it follows

[TABLE]

If Assumptions 1-2 hold, then following the steps detailed in [15, Appendix B] it is possible to derive the bound

[TABLE]

for the error due to the Taylor approximation at the prediction step.

If moreover Assumption 3 holds then the following quadratic bound can be derived instead

[TABLE]

-B Early termination error

In Algorithm 1 both during the prediction and the correction steps, the selected splitting method is halted after $P$ and $C$ steps, respectively.

FBS

Making use of the convergence results reported in (10) for the FBS, the prediction error after $P$ steps will be bounded as

[TABLE]

while the correction error after $C$ steps will be bounded as

[TABLE]

DRS

Applying instead Lemma 1 for the DRS it follows that the prediction error is bounded as

[TABLE]

where $\zeta(P)=\zeta_{\textrm{DR}}^{P}(1+\rho L)/(1+\rho m)$, while the correction error as

[TABLE]

In the following, we will indicate with function $\zeta(k)$, either $\zeta_{\textrm{FB}}^{k}$ for FBS, or $\zeta_{\textrm{DR}}^{P}(1+\rho L)/(1+\rho m)$ for DRS.

-C Overall error bound

Recall that our aim is to derive a bound for the error $\left\lVert\mathbold{x}_{k+1}-\mathbold{x}_{k+1}^{*}\right\rVert$ that depends on the error at the previous time instant $\left\lVert\mathbold{x}_{k}-\mathbold{x}_{k}^{*}\right\rVert$.

In order to do so we need one last result. Under Assumptions 1-2, it is possible to see that the solution mapping $t_{k}\mapsto\mathbold{x}^{*}(t_{k})$ satisfies Theorem 3, and thus that

[TABLE]

We are now ready to find the desired error bound.

Using the results above, we can now derive the following chain of inequalities

[TABLE]

which gives a bound for the prediction error. Moreover, employing (24) it holds

[TABLE]

and, under Assumptions 1-2, we can write

[TABLE]

where

[TABLE]

If moreover Assumption 3 holds as well, then we have

[TABLE]

-D Convergence proof

As the previous section showed, under Assumptions 1-2 we can characterize the dynamics of the error with the inequality

[TABLE]

Therefore, for the error to converge, it is necessary to choose the algorithm parameters in such a way that $\eta_{1}<1$. In particular, given the definition of $\eta_{1}$ and the rate $\zeta$, we need to make sure that the prediction and correction horizons, $P$ and $C$, and the step-size $\rho$, satisfy the condition $\eta_{1}<1$.

If this is the case, than it is clear that by the definition of $\eta_{0}$, Theorem 1 holds.

Suppose now that Assumption 3 holds as well. We have convergence of the algorithm if the error performed at each iteration does not exceed the error of the previous iteration, that is if there exists $\tau<1$ such that

[TABLE]

Therefore the error converges if $\eta_{1}<\tau<1$ which implies

[TABLE]

and if the initial condition $\mathbold{x}_{0}$ is chosen such that

[TABLE]

Using now the fact that

[TABLE]

we can compute

[TABLE]

which proves Theorem 2 when $k\to\infty$. $\blacksquare$