Stochastic Online Feedback Optimization for Networks of Non-Compliant Agents

TL;DR

This paper extends online feedback optimization to stochastic settings with non-compliant agents, providing theoretical error bounds and demonstrating applicability to power systems.

Contribution

It introduces a stochastic OFO framework accommodating non-compliant agents, with proven mean-square error bounds and practical application to power networks.

Findings

Derived mean-square error bounds for the stochastic OFO algorithm.

Validated the approach through application to power system scenarios.

Demonstrated robustness of the method with non-compliant agents.

Abstract

In several applications of online optimization to networked systems such as power grids and robotic networks, information about the system model and its disturbances is not generally available. Within the optimization community, increasing interest has been devoted to the framework of online feedback optimization (OFO), which aims to address these challenges by leveraging real-time input-output measurements to empower online optimization. We extend the OFO framework to a stochastic setting, allowing the subsystems comprising the network (the $\textit{agents}$) to be $\textit{non-compliant}$. This means that the actual control input implemented by the agents is a random variable depending upon the control setpoint generated by the OFO algorithm. Mean-square error bounds are obtained for the general algorithm and the theory is illustrated in application to power systems.

Tables

Table 1. Table 1: Comparison between stochastic OFO (S-OFO) and deterministic OFO (D-OFO) algorithms for various choices of { Φ n } \{\Phi_{n}\} .

$\left\{{\Phi }_n\right\}$	Support	S-OFO		D-OFO
		PC $\left(k{W}^2\right)$	VD $(p.u.)$	PC $\left(k{W}^2\right)$	VD $(p.u.)$
Beta$(4,2)$	$[0,1]$	$110$	$4\times {10}^{-4}$	$85.8$	$4.6\times {10}^{-4}$
Beta$(2,4)$	$[0,1]$	$273.7$	$1.8\times {10}^{-4}$	$216.6$	$2\times {10}^{-4}$
Uniform	$[0,1]$	$201$	$2.7\times {10}^{-4}$	$161.6$	$3.4\times {10}^{-4}$
Beta$(4,2)$	$[-0.5,1]$	$194.3$	$2.7\times {10}^{-4}$	$155.4$	$3.4\times {10}^{-4}$
Beta$(2,4)$	$[-0.5,1]$	$435.6$	$2.4\times {10}^{-4}$	$460.7$	$3.1\times {10}^{-4}$
Uniform	$[-0.5,1]$	$364.6$	$2.3\times {10}^{-4}$	$331.1$	$3.4\times {10}^{-4}$
Beta$(4,2)$	$[-1,1]$	$303.3$	$2.2\times {10}^{-4}$	$255.1$	$3\times {10}^{-4}$
Beta$(2,4)$	$[-1,1]$	$436.3$	$2.4\times {10}^{-4}$	$811.5$	$10.7\times {10}^{-4}$
Uniform	$[-1,1]$	$442$	$2.7\times {10}^{-4}$	$553.9$	$6.7\times {10}^{-4}$

Equations

$$y_n=h^{(n)}(x_n):=h(x_n,r_n)$$

$$x_n=\phi (u_n,{\Phi }_{n+1})$$

$$x_n$$

$$y_n$$

$$u_n^{\ast }\in {\mathrm{argmin}}_{u\in {\mathcal{U}}^{(n)}}\mathsf{E}[f^{(n)}(u,{\Phi }_{n+1})]$$

$$f^{(n)}(u,\Phi )=f_x^{(n)}(u,\Phi )+f_y^{(n)}(u,\Phi )$$

$$u_{n+1}={\mathsf{Proj}}_{{\mathcal{U}}^{(n)}}\{u_n-\alpha \nabla {\widehat{\text{Upgamma}}}^{(n)}\}$$

$$\mathop{\rm lim{\,}sup}_{N\to\infty}{\sf E}[\|u_{N}-u^{*}_{N}\|^{2}]\leq\alpha\upepsilon_{a}+\upepsilon_{b}+\upepsilon_{c}+\frac{b^{\text{\rm\tiny\ref{t:online_tracking_result_limsup}}}}{\alpha^{3/2}}\sqrt{\epsilon_{c}[\alpha\upepsilon_{a}+\upepsilon_{b}]}$$

$$\mathsf{E}[\Vert X{\Vert }^p]=\Vert X{\Vert }_p^p,\mathsf{E}[\Vert X{\Vert }^p\mid {\mathcal{F}}_n]=\Vert X{\Vert }_{p,n}^p$$

$${\nabla }_u{\tilde{f}}^{(n)}(u_n,{\Phi }_{n+1}):={\nabla }_u{f}^{(n)}(u_n,{\Phi }_{n+1})-{\nabla }_u{f}^{(n)}(u_n^{\ast },{\Phi }_{n+1})$$

$$\Vert \nabla g_y^{(n)}(y)-\nabla g_y^{(n)}(y^{\prime })\Vert$$

$$\Vert \nabla g_x^{(n)}(x)-\nabla g_x^{(n)}(x^{\prime })\Vert$$

$$\mathsf{E}[\Vert A_{n+1}{\Vert }^4\mid {\mathcal{F}}_n]\le {\sigma }_{\Delta }^4,\mathsf{E}[\Vert b_{n+1}{\Vert }^4\mid {\mathcal{F}}_n]\le {\sigma }_{\Delta }^4$$

$$\mathsf{E}[\Vert {\hat{y}}_n-y_n{\Vert }^4\mid {\mathcal{F}}_n]\le {\epsilon }_m^4,\mathsf{E}[\Vert {\hat{x}}_n-x_n{\Vert }^4\mid {\mathcal{F}}_n]\le {\epsilon }_m^4$$

$$\mathsf{E}[\Vert A_{n+1}^{\circ }-A_{n+1}{\Vert }^4\mid {\mathcal{F}}_n]\le {\epsilon }_m^4$$

$$u_{n+1}={\mathsf{Proj}}_{{\mathcal{U}}^{(n)}}\{u_n-\alpha {\nabla }_u{f}^{(n)}(u_n,{\Phi }_{n+1})\}$$

$${\nabla }_u{f}^{(n)}(u_n,{\Phi }_{n+1})={A_{n+1}}^{⊺}[C^{⊺}\nabla g_y^{(n)}(y_n)+\nabla g_x^{(n)}(x_n)]$$

$$\mathsf{E}[\Vert {\nabla }_u{\tilde{f}}^{(n)}(u_n,{\Phi }_{n+1}){\Vert }^2\mid {\mathcal{F}}_n]\le L_f^2\Vert {\tilde{u}}_n{\Vert }^2$$

$$\mathsf{E}[{\nabla }_u{f}^{(n)}(u_n,{\Phi }_{n+1}{)}^{⊺}(u_n-u_n^{\ast })\mid {\mathcal{F}}_n]\ge {\overline{\mu }}_f\Vert {\tilde{u}}_n{\Vert }^2$$

$$\mathsf{E}[\Vert {\nabla }_u{f}^{(n)}(u_n^{\ast },{\Phi }_{n+1})-\mathsf{E}[{\nabla }_u{f}^{(n)}(u_n^{\ast },{\Phi }_{n+1})]{\Vert }^2\mid {\mathcal{F}}_n]\le {\sigma }_f^2$$

$$\nabla {\widehat{\text{Upgamma}}}^{(n)}={A_{n+1}^{\circ }}^{⊺}{C}^{⊺}\nabla g_y^{(n)}({\hat{y}}_n)+{A_{n+1}^{\circ }}^{⊺}\nabla g_x^{(n)}({\hat{x}}_n)$$

$$\Vert {\tilde{u}}_N{\Vert }_2^2\le {{\rm Y}}_{\alpha }^N\Vert {\tilde{u}}_0{\Vert }_2^2+\sum _{k=0}^{N-1}{{\rm Y}}_{\alpha }^k[{\psi }_{N-k-1}^2+q_{\alpha }+2{\beta }_{N-k-1}]$$

$${\beta }_{n-1}$$

$$q_{\alpha }$$

$${\text{upepsilon}}_a$$

$${\text{upepsilon}}_b$$

$$u_n^{r\ast }\in {\mathrm{argmin}}_{u\in {\mathcal{U}}^{(n)}}\mathsf{E}[f^{(n)}(u,{\Phi }_{n+1})+\frac{1}{2}\eta \Vert u{\Vert }^2]$$

$$\nabla {\widehat{\text{Upgamma}}}^{(n)}={A_{n+1}^{\circ }}^{⊺}[C^{⊺}\nabla g_y^{(n)}({\hat{y}}_n)+\nabla g_x^{(n)}({\hat{x}}_n)]+\eta u_n$$

$$\mathop{\rm lim{\,}sup}_{N\to\infty}{\sf E}[\|u_{N}-u^{r*}_{N}\|^{2}]\leq\alpha\upepsilon^{r}_{a}+\upepsilon^{r}_{b}+\upepsilon^{r}_{c}+\frac{b^{\text{\rm\tiny\ref{t:online_reg}}}}{\alpha^{3/2}}\sqrt{\epsilon^{r}_{c}[\alpha\upepsilon^{r}_{a}+\upepsilon^{r}_{b}]}$$

$${\text{upepsilon}}_a^r$$

$${\text{upepsilon}}_b^r$$

$$f_y^{(n)}(u_n,{\Phi }_{n+1})$$

$$y_n$$

$$x_n$$

Full text

1. 1 Introduction

2. 2 Main Results

3. 2.1 Preliminaries

4. 2.2 Feedback-Based SGD and Tracking with Partial Information

5. 3 Numerical Experiments

6. 3.1 Toy Problems

7. 3.2 Real-Time Optimal Power Flow

8. 4 Conclusions

9. A Technical Proofs

1 Introduction

This paper aims to leverage stochastic gradient descent (SGD) algorithms to address optimization problems that typically arise in applications to networked systems such as communication systems, power grids and robotic networks [5, 10, 7].

In general, such systems are comprised of several interconnected entities (or agents) that are either controlled by a single central controller or multiple local controllers. The goal is for these controllers to adjust the system’s inputs so that its outputs are steered to an operating point that minimizes a possibly time-varying performance criterion.

Throughout this paper, the behavior of a networked system of ${\cal A}$ agents with a discretized timescale is modeled as

[TABLE]

in which $n$ is the time index, $\{y_{n}\}\subseteq\mathbb{R}^{m}$ represent the system’s output variables, $\{x_{n}\}\subseteq\mathbb{R}^{d}$ its input and $\{r_{n}\}\subseteq\mathbb{R}^{d}$ an exogenous disturbance process. The superscript on $h$ indicates that the system map is time varying; in this case, the time dependency is inherited from ${\bf r}:=\{r_{n}\}$.

We consider the setting in which the agents are not necessarily compliant with their controllers, in the sense that given a desired input sequence $\{u_{n}\}\subseteq\mathbb{R}^{d}$, the actual input sequence implemented by the agents is a random variable given by

[TABLE]

where $\varphi$ is a function that maps the desired input $u$ to the actual input $x$ and $\{\Phi_{n}\}\subseteq{\sf X}$ an independent and identically distributed (i.i.d.) sequence of random variables. Most of this paper is focused on the linear (or linearized) setting, wherein both the input-output mapping $h$ and the compliance mapping $\varphi$ are linear (or affine). In particular,

[TABLE]

where $A_{n+1}=A(\Phi_{n+1})\in\mathbb{R}^{d\times d}$, $b_{n+1}=b(\Phi_{n+1})\in\mathbb{R}^{d}$ and $C,D\in\mathbb{R}^{m\times d}$. A simple example is when $A_{n+1}=\hbox{\rm diag\thinspace}(\Phi_{n+1})$ and $b_{n+1}=0$, with $\Phi_{n+1}\in[0,1]^{d}$. This represents the case when the agents choose to implement any control input between [math] and the desired setpoint $u_{n}$, with probability determined by the distribution of $\Phi_{n+1}$.

The constrained optimization problem that is the object of study in this paper takes the form

[TABLE]

where the above expectation is taken with respect to the distribution of $\Phi_{n+1}$. The set ${\cal U}^{(n)}$ represents the system’s input constraints (e.g. physical or engineering constraints) at time $n$, while $f^{(n)}$ represents the system’s performance objectives. In this paper, we quantify the system performance via the additive model

[TABLE]

with $f_{x}$ and $f_{y}$ quantifying the performance in terms of the actual input and output of the system, respectively. Specifically, $f_{x}$ is a composition of functions $f^{(n)}_{x}(u,\Phi)=(g^{(n)}_{x}\circ\varphi)(u,\Phi)$, with $g^{(n)}_{x}(x)$ measuring the performance in terms of the system actual input; whereas $f^{(n)}_{y}(u,\Phi)=(g^{(n)}_{y}\circ h^{(n)}\circ\varphi)(u,\Phi)$, with $g^{(n)}_{y}(y)$ measuring the performance in terms of the system output. Throughout the paper, we consider the convex optimization case, where ${\cal U}^{(n)}$ is a convex set and $f^{(n)}(u,\Phi)$ is a convex function of $u$.

In typical engineering applications, obtaining solutions or critical points of (2) in real time might be infeasible due to several challenges. In general, one has no knowledge of the distributions of $\{\Phi_{n}\}$, the disturbance process ${\bf r}$ or even the full system model (1). Moreover, even if all of these pieces of information were available, computing the solution to (2) in real-time would be computationally expensive.

Past research has focused on addressing some of the above challenges through the framework of online feedback optimization (OFO) [3, 5, 8, 15]. Within this framework, the system’s input and ouput measurements are used to drive optimization algorithms in real-time in order to approximately solve (2). However, these previous works only considered compliant networks, in which $A(\Phi_{n})=I$ and $b(\Phi_{n})=0$ for each $n$. In this special case, (2) reduces to deterministic optimization: $\min_{u\in{\cal U}^{(n)}}f^{(n)}(u)$. A prototypical OFO approach to solving this problem is the approximate projected gradient scheme:

[TABLE]

where $\alpha>0$ is a constant step-size parameter and $\{\nabla\widehat{\Upgamma}^{(n)}\}$ is a sequence constructed from output measurements $\{{\hat{y}}_{n}\}$ with the goal of approximating $\{\nabla_{u}f^{(n)}(u_{n})\}$.

The extensive number of success stories of algorithms based upon stochastic approximation (SA) such as stochastic gradient descent (SGD) motivates for a stochastic extension of (3) to address non-compliant networks. In this case, the sequence $\{\nabla\widehat{\Upgamma}^{(n)}\}$ is constructed from input-output measurements $\{\hat{x}_{n},{\hat{y}}_{n}\}$ and aims to approximate $\{\nabla_{u}f^{(n)}(u_{n},\Phi_{n+1})\}$.

Contributions:  This paper extends the OFO framework to a stochastic setting, addressing situations in which the agents of the networked system of interest are non-compliant. For the general algorithm (3), we obtain the following upper bound on the mean-squared error (MSE) of estimates:

[TABLE]

in which b^{\text{\rm\tiny\ref{t:online_tracking_result_limsup}}} is a constant, $\upepsilon_{a}$ depends upon the volatility of $\{A_{n+1}\}$, $\upepsilon_{b}$ is related to the error in approximating $\{\nabla_{u}f^{(n)}(u_{n},\Phi_{n+1})\}$ through input-output measurements and $\upepsilon_{c}$ depends on the time-variability of the optimization problem (2).

In the setting of static optimization with full gradient information (i.e., $\upepsilon_{b}=\upepsilon_{c}=0$), the bound (4) coincides with the MSE bound expected for a SGD algorithm with constant step-size [6, 14]. If in addition the sequence $\{A_{n+1}\}$ is time-invariant (i.e., $\Phi_{0}\equiv\Phi_{n}$ for all $n$), it follows that $\upepsilon_{a}=0$, leading to convergence as expected from [5].

Relevant Literature:

*Online feedback optimization (OFO). *  Interest in the OFO framework has been growing within the controls community, particularly motivated by applications to power systems [5, 3, 8]. We point the reader to [3] for a more comprehensive survey of OFO.

The papers [5, 8] tackle the compliant case with additional inequality constraints. While [5] requires partial knowledge about the system map, [8] eliminates this limitation by leveraging gradient-free optimization methods. Another approach to overcome this challenge may be found in [15]. Similarly to the setting of the present paper, many of the references in OFO are restricted to convex optimization; see [11] for an extension to non-convex problems.

The reference [18], which is perhaps most closely related to our work, studies an application of primal-dual methods to static optimization without using input-output measurements. Their problem is much more challenging as they allow for the distribution of $\Phi_{n+1}$ to depend upon $u_{n}$. For a constant step-size $\alpha$, they obtain moment bounds independent of the step-size between $\{u_{n}\}$ and a performatively stable point $u^{*p}$. This generally differs from the desired optimal point $u^{*}$, but is optimal for the distribution that it induces on $\{\Phi_{n}\}$. Bounds on the error $\|u^{*p}-u^{*}\|$ depend upon the problem’s parameters (e.g., Lipschitz constants and strong convexity parameters) and are also not related to the step-size $\alpha$. Convergence of $\{u_{n}\}$ to $u^{*p}$ in the mean squared sense is only obtain for an adaptive step-size schedule $\{\alpha_{n}\}$ of the form $\alpha_{n}=\alpha_{0}/n$ with $\alpha_{0}>0$ a constant.

*Stochastic Approximation. *  Many machine learning and optimization algorithms are built upon stochastic approximation. This framework was born in the 1950s in the seminal work [16], and theoretical refinements have appeared ever since [14, 6, 13, 9, 17]. Most relevant to the present work is the paper [14] which obtains finite-time MSE bounds of order $O(\alpha)$ for SGD. Although the assumptions of this prior work are similar to the ones in the present paper, it concerned a static unconstrained optimization problem and assumed that $\{\nabla_{u}f^{(n)}(u_{n},\Phi_{n+1})\}$ could be measured directly.

Organization:  The remainder of the paper is organized in three additional sections. Section 2 sets the stage for analysis by going over the main assumptions and notation imposed in the paper as well as presenting the main results. These results are then illustrated through simple numerical experiments in Section 3. Moreover, Section 3 contains an application of the stochastic OFO algorithm in the context of power systems. Conclusions and directions for future research are contained in Section 4, while technical proofs can be found in the Appendix.

2 Main Results

2.1 Preliminaries

Notation:  We use $\|\,\cdot\,\|$ to denote the Euclidean norm for vectors and the induced operator norm for matrices. For a random variable (be it vector or matrix-valued) $X$, we write

[TABLE]

to denote its $L_{p}$ moments and conditional moments with respect to the filtration of the history generated by (3) up to time $n$: ${\cal F}_{n}\mathbin{:=}\{u_{0},A_{1},b_{1},A_{2},b_{2}\cdots,A_{n},b_{n}\}$. Moreover, we denote ${\bar{A}}\mathbin{:=}{\sf E}[A(\Phi_{n})]$ and ${\overline{b}}\mathbin{:=}{\sf E}[b(\Phi_{n})]$.

For each $n$, we use the following short-hand notation for tracking errors:   $\tilde{u}_{n}\mathbin{:=}u_{n}-u^{*}_{n}$,

[TABLE]

Assumptions:

(A1)  The functions $g^{(n)}_{x},g^{(n)}_{y}$ are strongly convex with parameter $\mu$ and continuously differentiable. Moreover, their gradients are Lipschitz continuous: there is $L_{g}<\infty$ such that for each $n$,

[TABLE]

for all $y,y^{\prime}\in\mathbb{R}^{m}$ and $x,x^{\prime}\in\mathbb{R}^{d}$.

(A2)  For each $n$, the set ${\cal U}^{(n)}\subset\mathbb{R}^{d}$ is convex and compact. Moreover, the sequence $\{{\cal U}^{(n)}\}$ is uniformly bounded: $b_{{\cal U}}\mathbin{:=}\sup_{n\geq 1}\sup_{u\in{\cal U}^{(n)}}\|u\|<\infty$.

(A3)  The sequence $\{\Phi_{n}\}$ is i.i.d. and there exists a constant $\sigma_{\Delta}<\infty$ such that for all $n$,

[TABLE]

Moreover, $\sup_{n}\|r_{n}\|\leq\sigma_{\Delta}$.

(A4)  The sequences of measurements $\{{\hat{y}}_{n},\hat{x}_{n}\}$ admit the bounds: for a constant $\varepsilon_{m}$ and all $n$,

[TABLE]

Moreover, it is assumed that the sequence $\{A_{n+1}\}$ can be recovered from the observations $\{\hat{x}_{n}\}$ with precision $\varepsilon_{m}$: there exists a matrix-valued sequence $\{A^{\circ}_{n+1}\}$ constructed from $\{\hat{x}_{n},u_{n}\}$ such that for all $n$,

[TABLE]

(A5)  The matrices $C{\bar{A}}$ and ${\bar{A}}$ are of full column rank.

Assumptions (A1) and (A2) are common in OFO [5, 3], while (A3) is imposed to ensure $\{\nabla_{u}f^{(n)}(u_{n},\Phi_{n+1})-{\sf E}[\nabla_{u}f^{(n)}(u_{n},\Phi_{n+1})]\}$ is a martingale difference sequence, a standard assumption in a vast part of the SA/SGD literatures [6, 14].

The conditions in (A4) are a slight strengthening of the measurement error conditions in the OFO literature: while a $L_{2}$ bound is typically assumed, we require a $L_{4}$ bound due to the presence of multiplicative noise in the algorithm [5, 8]. One special case in which $\{A_{n+1}\}$ can be recovered exactly from $\{u_{n},\hat{x}_{n}\}$ is when $A_{n+1}=\hbox{\rm diag\thinspace}(\Phi_{n+1})$, $b_{n+1}=0$ and $x_{n}=\hat{x}_{n}$, leading to $\{\hat{x}^{i}_{n}\}$ independent of $\{u^{j}_{n}:j\neq i\}$. For each $n$, $\{A_{n+1}\}$ is obtained via $A^{i,i}_{n+1}={A^{\circ}}^{i,i}_{n+1}={\hat{x}}^{i}_{n}/u^{i}_{n}$ for $1\leq i\leq d$. This choice is employed in the experiments surveyed in Section 3.

Assumption (A5) may seem strong at first but is imposed so that $\nabla_{u}f$ is strongly monotone in its first variable. In applications where this assumption is not satisfied, one could implement a regularized version of the algorithm (see the discussion following Thm. 2.3).

2.2 Feedback-Based SGD and Tracking with Partial Information

In the scenario in which full information is available at each time instant $n$, a projected stochastic gradient descent algorithm can be formulated as follows:

[TABLE]

Under the assumptions of this paper, the above choice satisfies assumptions that have been previously imposed within the SGD literature [14], and we have the following lemma.

Lemma 2.1**.**

Suppose that (A1)–(A3) hold.

(i)*  Then, $\nabla_{u}f^{(n)}$ is Lipschitz continuous in quadratic mean in its first variable: there is $L_{f}<\infty$ such that for each $n$,*

[TABLE]

(ii)*  If in addition (A5) holds, $\nabla_{u}f^{(n)}$ is strongly monotone in its first variable in conditional mean: there is a constant ${\overline{\mu}}_{f}>0$ such that for each $n$,*

[TABLE]

(iii)*  There exists a constant $\sigma_{f}$ such that the following holds for each $n$:*

[TABLE]

The proof of Lemma 2.1 is deferred to the Appendix.

In the setting of this paper, however, full information about the system model, its disturbances, and randomness is not available. Instead, we apply the approximate feedback-based SGD algorithm (3), in which

[TABLE]

Similarly to [5], we required exact knowledge of the Jacobian of $h$ for each $n$. Other implementations may approximate or estimate this quantity by linearization of the system map around an operating point or application of machine learning and gradient-free optimization techniques [4, 15, 8].

Theorem 2.2**.**

Suppose that (A1)–(A5) hold. Suppose in addition that the update rule (3) is implemented and that $\nabla\widehat{\Upgamma}^{(n)}$ satisfies (5) for each $n$. Then, the mean square tracking error admits the bound:

[TABLE]

where $\Upsilon_{\alpha}=1-2\alpha{\overline{\mu}}_{f}+\alpha^{2}L_{f}^{2}$, $\psi_{n}=\|u^{*}_{n+1}-u^{*}_{n}\|$,

[TABLE]

in which b^{\text{\rm\tiny\ref{t:one_step_couple_haobj}}} is a constant depending upon $b_{\cal U}$ and $L_{f}$.

In view of the finite-time bound in Thm. 2.2, the MSE bound (4) is obtained under additional conditions on the step-size and the time-variability of the problem.

Theorem 2.3**.**

Suppose the assumptions of Thm. 2.2 hold. Suppose, in addition, that there exists ${\bar{\gamma}}<\infty$ such that $\psi_{n}\leq{\bar{\gamma}}$ for all $n$. Then, the upper bound (4) holds for $\alpha<\frac{{\overline{\mu}}_{f}}{2L^{2}_{f}}$, in which

[TABLE]

The proofs of Thms. 2.2 and 2.3 are given in the Appendix.

Importance of (A5):  While assumption (A5) is necessary for the bounds in Thm. 2.3 to hold, it is likely that this assumption will not be satisfied in several applications. Instead, one could consider a regularized version of (2): for a small constant $\eta>0$,

[TABLE]

Analogous bounds as in Thm. 2.3 are obtained for the MSE, by solving (7) via (3) with

[TABLE]

Corollary 2.4**.**

Suppose that (A1)–(A4) hold and that there exists ${\bar{\gamma}}<\infty$ such that $\psi_{n}\leq{\bar{\gamma}}$ for all $n$. Then, the following holds when implementing (3) with (8) and $\alpha<\frac{\eta}{2(L_{f}+\eta)^{2}}$:

[TABLE]

in which

[TABLE]

and $b^{r}$ is a constant depending upon $b_{\cal U}$, $L_{f}$ and $\eta$.   \sqcap$$\sqcup

3 Numerical Experiments

3.1 Toy Problems

The first numerical example investigated was designed to test the tracking capability of the algorithm as well as the MSE bounds in Thm. 2.3. Consider (2) with

[TABLE]

For each $n$, $A_{n+1}=\hbox{\rm diag\thinspace}(\Phi_{n+1})$ in which $\Phi_{n+1}\in\mathbb{R}^{2}$ with entries $\Phi^{i}_{n+1}\sim\text{Unif}[0,1]$. The matrix $C$ is of the form $C=\hbox{\rm diag\thinspace}(\nu)$ where $\nu\in\mathbb{R}^{2}$ has entries $\nu^{i}\sim\text{Unif}[-5,0]$ and $D\in\mathbb{R}^{2\times 2}$ has entries sampled independently from $\text{Unif}[0,1]$. We note that (A5) is satisfied for this model so no regularization is needed.

To mimic real-life output measurements, the following observation model was adopted: ${\hat{y}}_{n}=y_{n}+w^{\bullet}_{n+1}$, in which $w^{\bullet}_{n}\sim N(0,I)$ for each $n$.

Tracking:  For a fixed simulation runlength $N=3\times 10^{5}$, the sequence $\{r_{n}\}$ was chosen as

[TABLE]

where the notation $\Pi,\Delta$ denotes the unit square and triangle waves, respectively. The frequencies and amplitudes were $\omega^{\circ}=[5,7]^{\intercal}\times 10^{-4}$ and $a_{n}=[10,\zeta_{n}]^{\intercal}$, in which

[TABLE]

The algorithm (3) was applied without regularization and with $\alpha=2\times 10^{-3}$. The projection sets were chosen to be the time-invariant and of the form ${\cal U}^{(n)}\equiv{\cal U}\{u:u^{\intercal}u\leq 9\}$.

Fig. 1 (a) shows plots of the sequences $\{u^{i}_{n},{u^{*}}^{i}_{n}:1\leq i\leq 2\}$ as functions of $n$. As expected from Thm. 2.3, the sequence $\{u_{n}\}$ reaches a neighborhood around the moving optimizer $\{u^{*}_{n}\}$ after a transient period and is able to track this optimal trajectory with a bounded error.

Estimation Error:  To test the estimation error bounds in Thm. 2.3, the sequence $\{r_{n}\}$ was chosen as $r_{n}\equiv r=[2,1]^{\intercal}$, so that ${\bar{\gamma}}=0$. For a fixed simulation runlength of $N=3\times 10^{5}$, several step-size values $\alpha\in[10^{-3},10^{-1}]$ were tested in independent experiments with initial condition $u_{0}=u^{*}$. Each experiment used the same noise sequence $\{A_{n}\}$ and $w^{\bullet}_{n}=0$, so that $\varepsilon_{m}=0$.

For each fixed $\alpha$, the empirical MSE $\|u^{\alpha}_{N}-u^{*}\|_{2}^{2}$ was estimated via Monte Carlo as follows: letting $\{u^{\alpha}_{n}\}$ be a sequence of estimates of $u^{*}$ obtained from (3) with step-size $\alpha$, the MSE was approximated empirically through

[TABLE]

where $N^{\circ}=\lfloor 0.5N\rfloor$.

Fig. 1 (b) shows a plot of $\{\|u^{\alpha}_{N}-u^{*}\|_{2}^{2}\}$ as a function of $\alpha$ in a logarithm scale. Also plotted for comparison with the expected bounds are the functions $\tau_{1}(\alpha)=\alpha$ and $\tau_{2}(\alpha)=\alpha^{2}$. We see that the empirical MSE scales with $\alpha$, as expected from the bounds in Thm. 2.3 for a time-invariant optimization problem with no observation noise (i.e., $\varepsilon_{m}={\bar{\gamma}}=0$).

3.2 Real-Time Optimal Power Flow

The next example aims to illustrate an application of the stochastic OFO framework within the context of power systems optimization. The goal is to optimize the operation of collections of distributed energy resources (DERs) in a power distribution network in real time. Similarly to what was done in [3, 5, 8], we frame this task as a time-varying optimal power flow (OPF) problem.

Setup:  We consider the IEEE 33-node test feeder [1], in which the controllable nodes (agents) were populated with photovoltaic (PV) systems (see Fig. 2 for details).

For each $n$, the desired input associated with the $i^{\text{th}}$ agent in the network is $u^{i}_{n}=(P^{i}_{n},Q^{i}_{n})^{\intercal}\in\mathbb{R}^{2}$, in which $P^{i}_{n}$ and $Q^{i}_{n}$ represent desired net active and reactive power injections to the distribution network, respectively. The actual net active and reactive power injections applied to the system by the $i^{\text{th}}$ agent is denoted $x^{i}_{n}=(P^{x,i}_{n},Q^{x,i}_{n})^{\intercal}\in\mathbb{R}^{2}$ so that $u_{n},x_{n}\in\mathbb{R}^{2{\cal A}}$. The compliance model was chosen so that $\{A_{n}\}$ could be recovered exactly from $\{\hat{x}_{n},u_{n}\}$ (see the discussion at the end of Section 2.1):

[TABLE]

where $\Phi_{n+1}\in\mathbb{R}^{\cal A}$. The system’s outputs are the voltage magnitudes measured at each of the agents $y_{n}\in\mathbb{R}^{\cal A}$.

As mentioned after (5), the OFO architecture proposed in this paper requires knowledge of the Jacobian of the system map. We approximate the mapping from $x_{n}$ to $y_{n}$ by a linear relationship based upon [4]: for each $n$,

[TABLE]

in which $\{Dr_{n}\}$ models the voltage contributions from the uncontrollable loads in the network.

While the above linearized model is used to implement the OFO algorithm, we note that the output measurements $\{{\hat{y}}_{n}\}$ are obtained by solving the exact nonlinear power flow equations for each $n$.

The regularized optimization problem (7) was considered with

[TABLE]

where $\kappa_{P},\kappa_{Q},\kappa_{y}$ are positive constants and $\{{\bar{P}}_{n}^{i}\}$ represents the real power available at the $i^{\text{th}}$ agent for each time instant $n$. The feasible set for the $i^{\text{th}}$ agent represents the PV inverter constraints:

[TABLE]

in which $S_{i,\max}$ is the inverter rating for agent $i$.

Results:  Load profiles for the uncontrollable nodes of the network consisted of real power usage data obtained from Smart* UMass apartment dataset. The sequences representing the real power available at each time instant at the PV generators of the controllable nodes $\{{\bar{P}}_{n}^{i}:1\leq i\leq{\cal A}\}$ were obtained from the Smart* UMass solar panel dataset [2].

Two OFO architectures were implemented for comparison: the stochastic algorithm proposed by this paper (8) (S-OFO) and the deterministic algorithm in [5] (D-OFO), in which $\nabla\widehat{\Upgamma}^{(n)}=C^{\intercal}\nabla g_{y}^{(n)}({\hat{y}}_{n})+\nabla g_{x}^{(n)}(\hat{x}_{n})+\eta u_{n}$. Both choices were implemented with $\alpha=5\times 10^{-2}$ and $\eta=10^{-3}$. The constants $\kappa_{P}=4$, $\kappa_{Q}=1$ and $\kappa_{y}=8$ were selected for the objective function.

The two algorithms were compared over $M=100$ independent experiments with equal load profiles and real available power curves, but different initial conditions and volatility: $\{\leftindex^{j}{u}_{0},\leftindex^{j}{\Phi}_{n}:1\leq j\leq M\}$. Performance was based upon two metrics: Approximate expected average power curtailment (PC) and approximate expected average voltage deviation (VD):

[TABLE]

Tab. 1 displays results from several experiments with different choices of $\{\Phi_{n}\}$. It is possible to see more power curtailing resulting from the D-OFO algorithm when there is more uncertainty at the controllable nodes. That is, when there is a larger presence of unexpected loads contributing to the power injections at these nodes.

Irrespective of the choice of $\{\Phi_{n}\}$, however, we consistently see much less voltage violations resulting from S-OFO when compared to its deterministic counterpart. A similar pattern was also observed for a different choice of objective that aimed to penalized power curtailment over voltage deviations, in which $k_{P}=8$, $k_{Q}=1$ and $k_{y}=4$.

This advantage of S-OFO over D-OFO is further illustrated in Fig. 3, which shows voltage magnitude profiles as functions of $n$ for the $4^{\text{th}}$ agent (node $29$) when $\Phi_{n}\sim\text{Unif[-1/2,1]}$. In this case, the stochastic algorithm leads to less voltage violations. Also plotted in Fig. 3 are lines at $1.05$ p.u. and $0.95$ p.u., which define a desired operating region for the system. While we do not enforce this through inequality constraints in the optimization problem as in [5], we include it in the plot to enable comparison.

We see from Fig. 3 that in the case of zero power curtailing (i.e., $P^{x,4}_{n}={\bar{P}}^{4}_{n}$), voltages repeatedly surpass the desired upper bound of $1.05$ p.u.. This is improved when the D-OFO algorithm is applied, but sporadic voltage violations are still apparent. Voltage magnitudes are more stabilized when the S-OFO algorithm is used.

4 Conclusions

This paper has extended the OFO framework to a stochastic setting, in which agents may not be compliant to the commands issued by the central controller or multiple local controllers. There are many paths for further research:

$\bullet$  The results in this paper restrict to the special case in which $\{\Phi_{n}\}$ is an i.i.d. sequence. Could we relax (A3) and allows for the the case in which $\{\Phi_{n}\}$ is a Markov chain?

$\bullet$  Prior work in SA/SGD has established that passing $\{u_{n}\}$ through a low-pass filter leads to better MSE bounds when $\{\Phi_{n}\}$ is a mixture of sinusoids [12]. Could we apply such techniques to the stochastic setting of this paper?

$\bullet$  Several references in the OFO literature consider explicit inequality constraints in the optimization problem and corresponding primal-dual methods. We conjecture that the analysis in this paper can be easily extended to address expectation inequality constraints. Could we extend this framework even further and allow for other types of constraints?

Appendix A Technical Proofs

We begin this section by establishing Lemma 2.1.

Proof of Lemma 2.1 (i)  For each $n$, denote $x^{*}_{n}=\varphi(u^{*}_{n},\Phi_{n+1})$ and $y^{*}_{n}=h^{(n)}(x^{*}_{n})$. By the chain rule we obtain,

[TABLE]

in which $\nabla\tilde{g}^{(n)}_{y}(y_{n})\mathbin{:=}\nabla g^{(n)}_{y}(y_{n})-\nabla g^{(n)}_{y}(y^{*}_{n})$. The Lipschitz continuity conditions in (A1) imply that $\|\nabla\tilde{g}^{(n)}_{y}(y_{n})\|\leq L_{g}\|CA_{n+1}\|\|\tilde{u}_{n}\|$, yielding \|\nabla_{u}\tilde{f}^{(n)}_{y}(u_{n},\Phi_{n+1})\|\leq b^{\text{\rm\tiny\ref{t:facts_f}}}\|A_{n+1}\|^{2}\|\tilde{u}_{n}\|, where b^{\text{\rm\tiny\ref{t:facts_f}}} is a constant. Applying similar steps as the ones outlined above to $\nabla_{u}f_{x}$ yields an analogous bound: $\|\nabla_{u}\tilde{f}^{(n)}_{x}(u_{n},\Phi_{n+1})\|\leq L_{g}\|A_{n+1}\|^{2}\|\tilde{u}_{n}\|$.

Now, by the triangle inequality we have

[TABLE]

Squaring both sides of the above equation and taking conditional expectations completes the proof with L_{f}=(b^{\text{\rm\tiny\ref{t:facts_f}}}+L_{g})\sigma_{\Delta}^{4} in view of the moment bounds in (A3).   \sqcap$$\sqcup

Proof of Lemma 2.1 (ii)  Let $x^{u}_{n}=\varphi(u,\Phi_{n+1})$, $y^{u}_{n}=h^{(n)}(x^{u}_{n})$, $x^{u^{\prime}}_{n}=\varphi(u^{\prime},\Phi_{n+1})$ and $y^{u^{\prime}}_{n}=h^{(n)}(x^{u^{\prime}}_{n})$. Moreover, let $\varrho^{t}(u,u^{\prime})=tu+(1-t)u^{\prime}$. Since $g^{(n)}_{x}$ and $g^{(n)}_{y}$ are assumed to be strongly convex, we have that for each $n,t$ and any $u,u^{\prime}\in{\cal U}^{(n)}$ that are ${\cal F}_{n}$-measurable,

[TABLE]

Now, since $u,u^{\prime}$ are ${\cal F}_{n}$-measurable and $\{\Phi_{n}\}$ is assumed i.i.d., we apply Jensen’s inequality to obtain the lower bound:

[TABLE]

Similarly, we have $\|{\bar{A}}(u-u^{\prime})\|^{2}\leq{\sf E}[\|A_{n+1}(u-u^{\prime})\|^{2}|{\cal F}_{n}]$.

Under (A5) it follows that the eigenvalues of the matrices ${\bar{A}}^{\intercal}C^{\intercal}C{\bar{A}}$ and ${\bar{A}}^{\intercal}{\bar{A}}$ are positive. Taking conditional expectations of both sides of (10) and using the above lower bound, we obtain strong convexity of $f$ in conditional mean, which implies part (ii) of Lemma 2.1 with

[TABLE]

\sqcap$$\sqcup

Proof of Lemma 2.1 (iii)  Let $x^{u}_{n}=\varphi(u,\Phi_{n+1})$ and $y^{u}_{n}=h(x^{u}_{n},r_{n})$. Using the triangle inequality, the assumed Lipschitz continuity of $\nabla g^{(n)}_{y}$ and $\nabla g^{(n)}_{x}$ in (A1) and Hölder’s inequality, we have the following for a constant $b^{\bullet}$,

[TABLE]

In view of the assumed bound $\|A_{n+1}\|_{4,n}\leq\sigma_{\Delta}$ in (A3) and the fact that $u^{*}_{n}$ is ${\cal F}_{n}$-measurable, it follows that there exists a potentially larger constant $b^{\circ}$ such that $\|x^{u}_{n}\|_{4,n}\leq b^{\circ}$ and $\|y^{u}_{n}\|_{4,n}\leq b^{\circ}$. Together with the above equation, this also leads to: \|\nabla_{u}f^{(n)}(u,\Phi_{n+1})\|_{2,n}\leq b^{\text{\rm\tiny\ref{t:facts_f}}}, in which b^{\text{\rm\tiny\ref{t:facts_f}}} is a constant.

Applications of Jensen’s inequality and the triangle inequality complete the proof for \sigma_{f}=2b^{\text{\rm\tiny\ref{t:facts_f}}}.   \sqcap$$\sqcup

The proof of Thm. 2.3 is largely based upon obtaining a contraction for the MSE $\|u_{n}-u^{*}_{n}\|_{2}^{2}$. The next result follows from parts (i) and (ii) of Lemma 2.1 and is important in defining a contractive factor for the MSE.

Corollary A.1**.**

Under the assumptions of Lemma 2.1 (ii), the following holds:

[TABLE]

where $\Upsilon_{\alpha}=1-2\alpha{\overline{\mu}}_{f}+\alpha^{2}L_{f}^{2}$.   \sqcap$$\sqcup

Next, we obtain an uniform bound on the error resulting from estimating the sequence of true gradients $\{\nabla_{u}f^{(n)}\}$ via the measurement-constructed sequence $\{\nabla\Upgamma^{(n)}\}$.

Lemma A.2**.**

Suppose that (A1)–(A4) hold. Then for a constant b^{\text{\rm\tiny\ref{t:estimation_haclL}}},

[TABLE]

Proof.

From the triangle inequality, we have

[TABLE]

in which

[TABLE]

We proceed to bound each term. Assumptions (A3), (A4) and the triangle inequality imply the upper bound: $\|A_{n+1}^{\circ}\|_{4,n}\leq\varepsilon_{m}+\sigma_{\Delta}$. Together with the assumed Lipschitz continuity of $\nabla g^{(n)}_{y}$ in (A1), the bounds in (A4) and Hölder’s inequality, we have the following: for a constant $b$,

[TABLE]

in which the uniform bound on $\|\nabla g^{(n)}_{y}(y_{n})\|_{4,n}$ follows from similar steps as in the proof of Lemma 2.1 (iii).

Repeating the above process for $\nabla g_{x}^{(x)}$ yields analogous bounds: ${\cal G}^{c}_{n}\leq b(\varepsilon_{m}+\sigma_{\Delta})\varepsilon_{m}$, ${\cal G}^{d}_{n}\leq b\varepsilon_{m}$. This completes the proof with b^{\text{\rm\tiny\ref{t:estimation_haclL}}}=2b(\varepsilon_{m}+\sigma_{\Delta}+1).   \sqcap$$\sqcup

The following shorthand notation will be used throughout the remainder of the appendix: let ${\cal E}_{n}\mathbin{:=}{\cal E}^{a}_{n}-\alpha{\cal E}^{b}_{n}$ and ${\cal M}_{n}\mathbin{:=}{\cal M}^{a}_{n}+{\cal M}^{b}_{n}$, in which

[TABLE]

Moreover, let $\mathcal{D}_{n}\mathbin{:=}(u^{*}_{n-1}-u_{n}^{*})^{\intercal}(u_{n}-u_{n-1}^{*})$.

We proceed to bounding the terms $E[{\cal E}_{n}{\cal M}_{n}\mid{\cal F}_{n-1}]$ and $E[\mathcal{D}_{n}\mid{\cal F}_{n-1}]$ in the next two lemmas. Together with the identity in Corollary A.1, bounds on the first term are crucial in establishing Lemma A.5, while bounds on the latter are used to obtain Thm. 2.2.

Lemma A.3**.**

Under (A1) and (A3), it follows that

[TABLE]

Consequently, ${\sf E}[{{\cal M}^{a}_{n}}^{\intercal}{\cal E}^{a}_{n}\mid{\cal F}_{n-1}]=0$.

Proof.

The result in (11) follows directly from the dominated convergence theorem (i.e., the order of expectation and differentiation can be exchanged in the second term of ${\cal M}^{a}_{n}$). To prove the remaining identity, we use the facts that ${\cal E}^{a}_{n}$ is ${\cal F}_{n}$-measurable and $\{\Phi_{n}\}$ is i.i.d.. Then, we have that

[TABLE]

\sqcap$$\sqcup

Lemma A.4**.**

Under (A1)–(A4), the following bounds hold:

(i)*  $\displaystyle|{\sf E}[{{\cal M}^{a}_{n}}^{\intercal}{\cal E}^{b}_{n}\mid{\cal F}_{n-1}]|\leq 2\sigma_{f}L_{f}b_{\cal U}$*

(ii)*  \displaystyle|{\sf E}[{{\cal M}^{b}_{n}}^{\intercal}{\cal E}^{a}_{n}\mid{\cal F}_{n-1}]|\leq 2b^{\text{\rm\tiny\ref{t:estimation_haclL}}}\varepsilon_{m}b_{\cal U}*

(iii)*  \displaystyle|{\sf E}[{{\cal M}^{b}_{n}}^{\intercal}{\cal E}^{b}_{n}\mid{\cal F}_{n-1}]|\leq 2b^{\text{\rm\tiny\ref{t:estimation_haclL}}}\varepsilon_{m}L_{f}b_{\cal U}*

Proof.

The Cauchy-Schwarz inequality yields

[TABLE]

Under (A2), we have that $\|{\cal E}^{a}_{n}\|_{2,n-1}\leq 2b_{\cal U}$, which also implies the following, via Lemma 2.1 (i): $\|{\cal E}^{b}_{n}\|_{2,n-1}\leq 2L_{f}b_{\cal U}$. Moreover, applications of (11), Lemma 2.1 (iii) and Lemma A.2 lead to the bounds: $\|{\cal M}^{a}_{n}\|_{2,n-1}\leq\sigma_{f}$ and \|{\cal M}^{b}_{n}\|_{2,n-1}\leq b^{\text{\rm\tiny\ref{t:estimation_haclL}}}\varepsilon_{m}, completing the proof.   \sqcap$$\sqcup

Lemma A.5**.**

Under the assumptions of Thm. 2.2, the following holds

[TABLE]

with q_{\alpha}=b^{\text{\rm\tiny\ref{t:one_step_couple_haobj}}}[\alpha^{2}(\xi+\sqrt{\xi})+\alpha\varepsilon_{m}], in which $\xi=\sigma_{f}+\varepsilon_{m}$ and b^{\text{\rm\tiny\ref{t:one_step_couple_haobj}}} is a constant depending upon $b_{\cal U}$ and $L_{f}$.

Proof.

Let $\upbeta_{n-1}=u^{*}_{n-1}-\alpha\nabla_{u}{\sf E}[f^{(n-1)}(u^{*}_{n-1},\Phi_{n})]$. By linearity of the projection operator and the fact that $u^{*}_{n-1}$ satisfies a fixed point equation, we have that

[TABLE]

in which the second inequality follows from the non-expansiveness property of the projection operator and ${\cal H}_{n}=2\alpha^{2}[{{\cal M}^{a}_{n}}^{\intercal}{\cal E}^{b}_{n}+{{\cal M}^{b}_{n}}^{\intercal}{\cal E}^{b}_{n}]-2\alpha{{\cal M}^{b}_{n}}^{\intercal}{\cal E}^{a}_{n}$.

Upon taking conditional expectations of both sides and applying Corollary A.1 and Lemma A.3, we obtain the upper bound

[TABLE]

It remains to bound the last two terms in the right hand side. Using the triangle inequality and the bounds in Lemma A.4, we have

[TABLE]

in which b^{\text{\rm\tiny\ref{t:one_step_couple_haobj}}} is a constant depending upon $b_{\cal U}$ and $L_{f}$. The remaining term is bounded similarly: applying (11), Lemma 2.1 (iii) and Lemma A.2 we have \|{\cal M}_{n}\|_{2,n-1}\leq\sigma_{f}+b^{\text{\rm\tiny\ref{t:estimation_haclL}}}\varepsilon_{m}, which completes the proof.   \sqcap$$\sqcup

Lemma A.6**.**

Under the assumptions of Thm. 2.2, the following bound holds:

[TABLE]

in which

[TABLE]

where $q_{\alpha}$ is given by Lemma A.5.

Proof.

The Cauchy-Schwarz inequality yields

[TABLE]

where the last inequality was obtained from an application of Lemma A.5 and the triangle inequality. Repeating this process recursively yields the desired result.   \sqcap$$\sqcup

Equipped with Lemma A.5 and Lemma A.6, we are ready to establish the finite-time bounds in Thm. 2.2.

Proof of Thm. 2.2  Expanding the square in $\|\tilde{u}_{n}-u^{*}_{n-1}+u^{*}_{n-1}\|^{2}_{2,n-1}$ yields

[TABLE]

where the last bound follows from applications of Lemma A.5 and Lemma A.6. Taking expectations of both sides and repeating this process recursively yields the final result.   \sqcap$$\sqcup

Finally, Thm. 2.3 follows as a corollary to Thm. 2.2.

Proof of Thm. 2.3  Upon choosing $\alpha<\frac{{\overline{\mu}}_{f}}{2L^{2}_{f}}$, it follows that $\Upsilon_{\alpha}\leq\sqrt{\Upsilon_{\alpha}}<1$ and $\Upupsilon_{\alpha}\leq 1-{\overline{\mu}}_{f}\alpha$. Thus, we obtain the following, via the geometric series formula: for a fixed $N$ and each $k$,

[TABLE]

Substituting the above identity into (6), using the geometric series formula once more and taking the limit supremum of both sides yields

[TABLE]

The proof is complete upon using the fact that $\sqrt{\Upupsilon_{\alpha}}<1$ to bound the last term and substituting $q_{\alpha}$ from Thm. 2.2:

[TABLE]

\sqcap$$\sqcup