#  An Online Feedback Optimization Approach to Voltage Regulation in Inverter-
Based Power Distribution Networks

Alejandro D. Domínguez-García, Madi Zholbaryssov, Temitope Amuda, and
Olaoluwapo Ajala Domínguez-García, Amuda, and Ajala are with the Department of
Electrical and Computer Engineering of the University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA. Email: {aledan, tamuda2,
<EMAIL_ADDRESS>is with Typhoon HIL Inc., Somerville, MA
02143, USA. Email: Email<EMAIL_ADDRESS>This work was
supported in part by the U.S. Department of Energy’s Office of Energy
Efficiency and Renewable Energy (EERE) under Solar Energy Technologies Office
(SETO) Agreement Number EE0009025; and the C3.ai Digital Transformation
Institute.

###### Abstract

We address the problem of controlling the reactive power setpoints of a set of
distributed energy resources (DERs) in a power distribution network so as to
mitigate the impact of variability in uncontrolled power injections associated
with, e.g., renewable-based generation. We formulate the control design
problem as a stochastic optimization problem, which we solve online using a
modified version of a projected stochastic gradient descent (PSGD) algorithm.
The proposed PSGD-based algorithm utilizes sensitivities of changes in bus
voltage magnitudes to changes in DER reactive power setpoints; such
sensitivities are learned online via a recursive least squares estimator
(rLSE). To ensure proper operation of the rLSE, the sequence of incremental
changes in DER reactive power setpoints needs to be persistently exciting,
which is guaranteed by a mechanism built into the controller. We analyze the
stability of the closed-loop system and showcase controller performance via
numerical simulations on the IEEE 123-bus distribution test feeder.

##  I Introduction

Driven by a global effort to decarbonize the electricity sector, electric
power distribution networks are undergoing major transformations in both power
production and demand. These transformations include the massive integration
of variable generation, e.g., photovoltaic installations, new types of loads,
e.g., plug-in electric vehicles, and storage devices; such energy assets are
commonly referred to as distributed energy resources (DERs). Coupled with more
frequent (and possibly random) network topology changes, and the increasing
uncertainty in load demand, the rapidly increasing deployment of DERs poses
numerous operational challenges for which existing control schemes in power
distribution networks are not well equipped to handle. In light of this, the
objective of this paper is to address one such challenge—ensuring effective
regulation of voltage magnitude across all buses of the power distribution
network.

Currently, voltage regulation in power distribution networks is accomplished,
for the most part, through tap-changing under-load transformers (TCULs) and
fixed/switched capacitor banks (see, e.g., [1, p. 16]). However, while these
devices are effective in managing slow changes in voltage (minutes to hours),
they are not suitable for managing fast voltage fluctuations (seconds to
minutes) arising, e.g., from rapid changes in renewable-based power
generation. This problem can be effectively addressed by controlling the
reactive power injected into the distribution network by power-electronic
inverter-interfaced DERs—a solution that has been actively pursued in the last
decade (see, e.g, [2, 3] and the references therein).

In this paper, we also pursue the idea of utilizing inverter-interfaced DERs
for voltage regulation. Building on our earlier work on data-driven control
algorithms in [3, 4, 5], we design a voltage regulation scheme that does not
rely on an a priori known model of the system to be controlled; instead, the
scheme utilizes data to estimate such a model online, while simultaneosly
executing a feedback control algorithm. We demonstrate that the proposed
approach is adaptive to constantly varying system conditions and disturbances,
and is capable of utilizing the reactive power support capabilities of
inverter-based resources to effectively provide voltage regulation at a fast
time-scale (i.e., milliseconds to seconds).

The control design problem is cast as a stochastic optimization problem whose
goal is to determine the reactive power setpoints of DERs so as to minimize
the expectation of bus voltage deviations from their nominal values. In order
to solve this problem, we utilize a projected stochastic gradient descent
(PGSD) algorithm (see e.g., [6]). Proper initialization and execution of the
algorithm for a single step essentially results in a feedback controller that
utilizes measured voltage deviations to adjust reactive power setpoints of
DERs. In its basic form, the algorithm relies on knowing the sensitivities of
changes in bus voltages with respect to changes in reactive power injections.
Instead of obtaining these sensitivities offline via a model, we design a
recursive least squares estimator (rLSE) that learns them online by using
real-time voltage measurement data and the sequence of reactive power
setpoints generated by the controller. The rLSE is executed in parallel with
the controller, and in order to ensure proper operation of the estimator, the
sequence of incremental changes in DER reactive power setpoints needs to be
persistently exciting. To this end, we modify the basic control algorithm to
include a mechanism that ensures this.

Most of the existing literature on utilization of inverter-interfaced DERs for
voltage regulation rely on the use of exact models of the network (see, e.g.,
[2] and the references therein). However, these methodologies have practical
limitations in that, due to the limited number of sensors in power
distribution networks, the models are hard to obtain in practice. More
recently, there have been several papers that propose the use of data-driven
techniques for addressing the voltage regulation problem in power distribution
networks, as well as in other related control and estimation problems. For
example, the authors of [4] present a data-driven framework for coordinating
the active and reactive power injections of DERs to provide voltage regulation
in radial networks. The proposed scheme utilizes estimates of network voltage
sensitivities obtained online by fitting measurements to the so-called
LinDistFlow model; by contrast, the method we propose here is not restricted
to radial networks, and it does not assume any particular model structure. The
authors of [3] propose a data-driven voltage regulation approach based on the
estimated topology and line parameters of radial power distribution systems.
The authors of [5] propose a model-free control scheme for regulating voltage,
frequency, and line flows that assumes no prior information on the system,
and, hence, is highly adaptable to intermittent operating conditions.
Additionally, a number of papers (see, e.g., [7, 8, 9]) have proposed data-
driven approaches that use measurement data to directly perform controller
synthesis without system identification. The use of linear estimation
techniques to learn sensitivities of various power system state variables to
control inputs has also been exploited in numerous applications in bulk power
system monitoring and control [10, 11, 12].

The framework proposed in this paper is closely related to that presented in
[13], which proposes a model-free, real-time optimal power flow solver through
feedback optimization. Similar to the approach we adopt here, [13] proposes a
data-driven online approach to learn a model of input–output sensitivities.
The objective function of the online feedback optimization formulated is the
operational cost of the inputs, and the inequality constraints are the
operational limitations on the inputs. The sensitivity estimate at each
iteration is computed using both the previous estimate and the current
measurements only, and a projected gradient solution method for the
optimization problem, with persistent excitation of the control inputs
performed at each iteration. By contrast, in our work, we formulate the
control problem as a stochastic optimization problem, where the objective
function is the expected value of the norm of the deviation of the outputs
from their nominal values. In addition, the sensitivity estimate in our work
is computed using all previous measurements with less weights assigned to
older measurements. Finally, while the proposed scheme in [13] does not
account for the constraints on control inputs when attempting to generate
persistently exciting inputs, our work does take this important issue into
account.

The remainder of this paper is organized as follows. We begin by formulating
the voltage regulation problem in Section II. In order to solve such problem,
Section III describes a model-based feedback optimization controller that
relies on sensitivities of voltage deviations to reactive power injections;
such sensitivities are learnt online via a rLSE estimator described in Section
IV. Then, our proposed controller, which is based on combining the ideas in
Section III and Section IV, is presented in Section V. In Section VI, we
present the main results of our convergence analysis, whereas in Section VII,
we provide numerical simulation results that demonstrate the effectiveness of
our control scheme. Concluding remarks are given in Section VIII.

##  II Problem Formulation

Consider a power distribution network with $n+1$ buses indexed by the elements
in $\mathcal{N}=\\{0,1,2,\dots,n\\}$, where the $0$ element corresponds to the
bus at which the network is connected to an external system, e.g., a sub-
transmission grid. Assume the network has $m$ reactive-power-capable DERs
indexed by the elements in $\mathcal{C}=\\{1,2\dots,m\\}$. Let $v_{i}(t)$
denote the magnitude of the phasor associated with the voltage at bus
$i,~{}i\in\mathcal{N}$, at time $t$. Assume that $v_{0}(t)=V_{0}$, where
$V_{0}$ is a positive constant, for all $t\geq 0$, and define
$v(t):=\big{[}v_{1}(t),v_{2}(t),\dots,v_{n}(t)\big{]}^{\top}$. Also, let
$q_{i}(t)$ denote the reactive power injected into the network by reactive-
power capable DER $i,~{}i\in\mathcal{C}=\\{1,2\dots,m\\}$, at time $t$, and
define $q(t):=\big{[}q_{1}(t),q_{2}(t),\dots,q_{m}(t)\big{]}^{\top}$.

Assume that, initially, $q(t_{0})=q_{0}$, where $q_{0}\in\mathds{R}^{m}$ is
given. Then, the objective is to adjust $q(t_{k})=:q_{k},~{}k\geq 1,$ via a
feedback controller so as to regulate $v(t)$ to some
$v^{*}=\big{[}v_{1}^{*},v_{2}^{*},\dots,v_{n}^{*}\big{]}^{\top},$ where
$v_{i}^{*}>0,~{}i=1,2,\dots,n,$ denotes some nominal value. Let
$\\{t_{k}^{+}\\}_{k\geq 1}$, $t_{k}<t_{k}^{+}<t_{k+1}$, denote the sequence of
time instants at which $v(t)$ is measured, and define $v_{k}:=v(t_{k}^{+})$.
Then, we have that

$\displaystyle v_{k}=h\big{(}q_{k},w_{k}\big{)},$ (1)

where $w_{k}=\big{[}w_{1,k},w_{2,k},\dots,w_{d,k}\big{]}^{\top}$, with
$w_{j,k}:=w_{j}(t_{k}),~{}j=1,2,\dots,d,$ representing exogenous disturbances
associated with, e.g., DER active power generation, load active and reactive
power demand, and network parameters.

###### Remark 1

In general, it is difficult to analytically characterize the function
$h(\cdot,\cdot)$ as this essentially entails obtaining an analytical solution
to the power flow problem, which is defined by a set of nonlinear equations
describing the active and reactive power balance at each bus (see, e.g., [14,
pp. 323-330]). To be more specific, these equations map the magnitudes and
phase angles of the phasors associated with bus voltages to the active and
reactive power injections at all buses; hence one would need to invert such
mapping to analytically characterize $h(\cdot,\cdot)$. Therefore, in the
remainder we will assume that $h(\cdot,\cdot)$ is not known.

Assume that values taken by $w_{k},~{}k\geq 1$, are of the form

$\displaystyle w_{k}=w_{k}^{\circ}+\xi_{k},$ (2)

with the value taken by $w_{k}^{\circ}$ known and slowly changing as $k$
evolves, and where $\xi_{k}$ is not a priori known but can be described by a
random vector $\Xi_{k}\in\mathds{R}^{d}$. In the remainder, we assume that the
random vectors $\Xi_{k},~{}k\geq 0,$ are independent and identically
distributed (i.i.d.) with zero mean, i.e.,

$\displaystyle\text{E}\big{[}\Xi_{k}\big{]}=0,\quad k\geq 1,$ (3)

and whose covariance matrix is such that

$\displaystyle\text{E}\big{[}\|\Xi_{k}\|_{2}^{2}\big{]}\leq\sigma,$ (4)

where $\sigma$ denotes some positive constant. Now, since the values taken by
$w_{k}^{\circ},~{}k\geq 1$, are known, we can equivalently describe the
relation in (1) as follows:

$\displaystyle v_{k}=h_{k}(q_{k},\xi_{k}),$ (5)

where $h_{k}(q_{k},\xi_{k})\coloneqq h(q_{k},w_{k}^{\circ}+\xi_{k})$. In the
remainder, we will assume $q_{0}$ is such that

$\displaystyle v_{0}$ $\displaystyle=h(q_{0},w_{0}^{\circ})$
$\displaystyle=v^{*},$ (6)

with the value taken by $w_{0}^{\circ}$ known, i.e., initially there is no
uncertainty in the value that the exogenous disturbances takes, and the
reactive powers injected by the DERs are such that the magnitudes of all bus
voltages are equal to their nominal values.

##  III Model-Based Feedback Optimization

In order to determine the value of $q_{k},~{}k\geq 1$, consider the following
optimization problem:

$\displaystyle\underset{\varphi_{k}}{\mbox{minimize}}$
$\displaystyle\quad\frac{1}{2}\text{E}_{\Xi_{k}}\bigg{[}\big{\|}v^{*}-h_{k}\big{(}\varphi_{k},\Xi_{k}\big{)}\big{\|}_{2}^{2}\bigg{]}$
(7a) subject to
$\displaystyle\quad\underline{q}_{k}\leq\varphi\leq\overline{q}_{k},$ (7b)

where $\text{E}_{\Xi_{k}}[\cdot]$ denotes expectation over the distribution of
$\Xi_{k}$, and $\underline{q}_{k}$ and $\overline{q}_{k}$ denote the lower and
upper bounds, respectively, on the reactive power that can be injected by DERs
into the network. Then, a local minimum of (7), which we denote by
$\varphi^{*}_{k}$, can be approximately obtained as
$\varphi^{*}_{k}\approx\lim_{r\to\infty}\varphi^{r}_{k}$, with the evolution
of $\varphi^{r}_{k},$ governed by

$\displaystyle\varphi^{r+1}_{k}=\Bigg{[}\varphi^{r}_{k}+$
$\displaystyle\gamma_{k}\Bigg{(}\frac{\partial
h_{k}(\varphi,\xi)}{\partial\varphi}\Bigg{|}_{\varphi=\varphi^{r}_{k},~{}\xi=\xi_{k}^{r}}\Bigg{)}^{\top}$
$\displaystyle\times\Big{(}v^{*}-h_{k}\big{(}\varphi^{r}_{k},\xi_{k}^{r}\big{)}\Big{)}\Bigg{]}_{\underline{q}_{k}}^{\overline{q}_{k}},\quad
r\geq 0,$ (8)

where $\gamma_{k}$ is a positive constant, and $\xi_{k}^{r},~{}r\geq 0,$
denote samples from the distribution of $\Xi_{k}$. The algorithm in (8) is
essentially a PSGD algorithm for the problem in (7) (see, e.g., [6] and the
references therein).

We initialize the PSGD algorithm in (8) as follows

$\displaystyle\varphi^{0}_{k}=q_{k-1},$ $\displaystyle\xi_{k}^{0}=\xi_{k-1},$
(9)

where $\xi_{k-1}$ denotes the realized value of $\Xi_{k-1}$, and instead of
running it to completion, we will only execute one iteration and set
$q_{k}=\varphi^{1}_{k}$. Then, it follows that the evolution of $q_{k}$ is
governed by

$\displaystyle q_{k}=\Bigg{[}q_{k-1}+$
$\displaystyle\gamma_{k}\Bigg{(}\frac{\partial
h_{k}(\varphi,\xi)}{\partial\varphi}\Bigg{|}_{\varphi=q_{k-1},~{}\xi=\xi_{k-1}}\Bigg{)}^{\top}$
$\displaystyle\times\Big{(}v^{*}-h_{k}\big{(}q_{k-1},\xi_{k-1}\big{)}\Big{)}\Bigg{]}_{\underline{q}_{k}}^{\overline{q}_{k}}.$
(10)

The one-step PSGD-based algorithm proposed in (10) for solving the problem in
(7) is a special case of the setting in [15], which proposes a framework to
sequentially solve stochastic minimization problems with slowly varying cost
functions that are convex.

###### Remark 2

We assume for subsequent developments that the problem (7) changes slowly with
$k$, while the one-step PSGD-based algorithm is executed, which entails
assuming that the $h_{k}(\cdot,\cdot)$’s, $\underline{q}_{k}$ and
$\overline{q}_{k}$ vary slowly with $k$. However, the algorithm will also work
in the presence of large but infrequent variations due to, e.g., weather
changes, as long as there are sufficiently long time intervals during which
the algorithm is able to converge.

Note that by using $\xi_{k-1}$ to initialize the PSGD-based algorithm, we are
using a sample from $\Xi_{k-1}$; however, recall that we have assumed that the
$\Xi_{l}$’s are i.i.d., so we are effectively sampling from the distribution
of $\Xi_{k}$. Then, since we have assumed that $h_{k}(\cdot,\cdot)$ changes
slowly with $k$, we can make the following approximations:

$\displaystyle h_{k}\big{(}q_{k-1},\xi_{k-1}\big{)}$ $\displaystyle\approx
h_{k-1}\big{(}q_{k-1},\xi_{k-1}\big{)}$ $\displaystyle=v_{k-1},$ (11)
$\displaystyle\frac{\partial
h_{k}(\varphi,\xi)}{\partial\varphi}\Bigg{|}_{\varphi=q_{k-1},~{}\xi=\xi_{k-1}}$
$\displaystyle\approx\frac{\partial
h_{k-1}(\varphi,\xi)}{\partial\varphi}\Bigg{|}_{\varphi=q_{k-1},~{}\xi=\xi_{k-1}}.$
(12)

Thus, if we assume measurements of $v_{k},~{}k\geq 1,$ are available, we can
update the value of $q_{k}$ as follows:

$\displaystyle
q_{k}=\Big{[}q_{k-1}+\gamma_{k}S_{k-1}^{\top}\big{(}v^{*}-v_{k-1}\big{)}\Big{]}_{\underline{q}_{k}}^{\overline{q}_{k}},$
(13)

where

$\displaystyle S_{k-1}=\frac{\partial
h_{k-1}(\varphi,\xi)}{\partial\varphi}\Bigg{|}_{\varphi=q_{k-1},~{}\xi=\xi_{k-1}}\in\mathds{R}^{n\times
m}.$ (14)

Note that in order to execute (13), we would need to compute $\partial
h_{k-1}(\varphi,\xi)/\partial\varphi$ for $\varphi=q_{k-1}$ and
$\xi=\xi_{k-1}$. This computation can be done by manipulating the power flow
Jacobian without necessarily solving the power flow equations; however, it
requires knowing the value of $\xi_{k-1}$, which we have assumed it is not a
priori known. To circumvent this issue, instead of using $S_{k-1}$ in (13), we
will use an estimate obtained using measurements of
$\big{\\{}q_{l},v_{l}\big{\\}}_{l=0}^{k-1}$; the construction of such an
estimate is detailed next.

##  IV Online Matrix Sensitivity Estimator

Define $\Delta v_{k}:=v_{k}-v_{k-1}$, $\Delta q_{k}:=q_{k}-q_{k-1}$ and
$\Delta\xi_{k}:=\xi_{k}-\xi_{k-1}$; then, by using (5), we have that

$\displaystyle v_{k-1}+\Delta v_{k}=h_{k}\big{(}q_{k-1}+\Delta
q_{k},\xi_{k-1}+\Delta\xi_{k}\big{)},$ (15)

Now, by expanding the right-hand side of (15) about
$\big{(}q_{k-1},\xi_{k-1}\big{)}$ using the Taylor series expansion, and using
the approximation in (12), it follows that

$\displaystyle\Delta v_{k}=S_{k-1}\Delta q_{k}+\epsilon_{k},$ (16)

where $\epsilon_{k}$ represents higher order terms in $\Delta q_{k}$, all the
terms in $\Delta\xi_{k}$, and the error associated with the use of (12). Then,
by assuming that $\epsilon_{k}$ is much smaller that $S_{k-1}\Delta q_{k}$, we
have that

$\displaystyle\Delta v_{k}\approx S_{k-1}\Delta q_{k}.$ (17)

Now let $\widehat{S}^{*}_{k}$ denote an estimate of $S_{k}$, which we can
obtain using the method of least squares as follows:

$\displaystyle\widehat{S}_{k}^{*}=\underset{S}{\operatorname{argmin}}~{}\ell(S),$
(18)

where

$\ell(S):=\sum_{l=1}^{k}\lambda^{k-l}\big{\|}\Delta v_{l}-S\Delta
q_{l}\big{\|}_{2}^{2},$

with $\lambda\in(0,1)$ denoting the forgetting factor that allows the
estimator to assign exponentially less weight to older measurements and adapt
to changes in operating conditions. The solution to the least-squares problem
(18) can be derived as follows. We first compute the gradient of the loss
function $\ell(S)$ in (18), which yields the following expression:

$\displaystyle\nabla\ell(S)=-2\sum_{l=1}^{k}\lambda^{k-l}\big{(}\Delta
v_{l}-S\Delta q_{l}\big{)}\Delta q_{l}^{\top}.$ (19)

Then, we equate the gradient (19) to zero and obtain

$\displaystyle\sum_{l=1}^{k}\lambda^{k-l}\Delta v_{l}\Delta
q_{l}^{\top}=S\left(\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top}\right).$ (20)

Now in order to solve for $S$ in (20), we need to invert the matrix that
multiplies $S$. Since this matrix is the sum of $k$ $(m\times m)$-dimensional
rank-one matrices, sufficient conditions for ensuring its invertibility are
that: i) $k\geq m$, and ii) the sequence $\\{\Delta q_{l}\\}_{l=1}^{k}$ is
persistently exciting.111A discrete-time signal $x[t]$, $t=1,2,\dots,$ is
persistently exciting if, for every $k$, there exist an integer $l$ and
constants $\varrho_{1},\varrho_{2}>0$ such that the matrices
$\varrho_{1}I-\sum_{t=k}^{k+l}x[t]x[t]^{\top}$ and
$\sum_{t=k}^{k+l}x[t]x[t]^{\top}-\varrho_{2}I$ are positive definite [16].
Assuming these two conditions are satisfied, we can now solve for $S$ in (20)
and obtain $\widehat{S}_{k}^{*},~{}k\geq m$:

$\displaystyle\widehat{S}_{k}^{*}$
$\displaystyle=\left(\sum_{l=1}^{k}\lambda^{k-l}\Delta v_{l}\Delta
q_{l}^{\top}\right)\Bigg{(}\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top}\Bigg{)}^{-1}.$ (21)

The issue of ensuring that the sequence $\\{\Delta q_{l}\\}_{l=1}^{k}$ is
persistently exciting is addressed in the next section. However, even if this
can be addressed satisfactorily, the estimate $\widehat{S}_{k}^{*}$ obtained
using (21) is only valid for $k\geq m$; this means that in practice we would
have to wait for $k$ steps before we can obtain our first estimate. To address
this, we use an algorithm that will recursively generate a sequence
$\big{\\{}\widehat{S}_{k}\big{\\}}_{k\geq 1}$ that can be shown to converge to
$\widehat{S}_{k}^{*}$ for $k$ large enough. The update rules for this
recursive algorithm are as follows:

$\displaystyle\widehat{S}_{k}=$
$\displaystyle~{}\widehat{S}_{k-1}+\Big{(}\Delta v_{k}-\widehat{S}_{k-1}\Delta
q_{k}\Big{)}\Delta q_{k}^{\top}F_{k},$ $\displaystyle F_{k}=$
$\displaystyle~{}\lambda^{-1}F_{k-1}$
$\displaystyle-\frac{\lambda^{-2}}{1+\lambda^{-1}\Delta
q_{k}^{\top}F_{k-1}\Delta q_{k}}F_{k-1}\Delta q_{k}\Delta
q^{\top}_{k}F_{k-1},$ (22)

with

$\displaystyle\widehat{S}_{0}$ $\displaystyle=\frac{\partial
h(\varphi,w)}{\partial\varphi}\Bigg{|}_{\varphi=q_{0},~{}w=w_{0}^{\circ}},$
$\displaystyle F_{0}$ $\displaystyle=I_{m},$ (23)

where $I_{m}$ denotes the $(m\times m)$-dimensional identity matrix. The
algorithm is obtained by using the matrix inversion lemma to recursively
invert the matrix $\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top}$; see Appendix VIII for the derivation.

###### Remark 3

Note that in order to initialize the algorithm in (22), we need to compute
$\partial h(\varphi,w)/\partial\varphi$ for $\varphi=q_{0}$ and
$w=w_{0}^{\circ}$. Unlikely the computation of $S_{k-1}$ in (14), in this
case, we can indeed compute $\widehat{S}_{0}$ by manipulating the power flow
Jacobian without necessarily solving the power flow equations because $q_{0}$
and $w_{0}^{\circ}$ are known. Furthermore, the initialization of
$\widehat{S}_{0}$ as in (23) is not crucial as we will see in the numerical
simulation results presented in Section VII.

###### Remark 4

Note that if the algorithm in (22) were to be executed for $k>m$, with
$\widehat{S}_{m}=\widehat{S}_{m}^{*}$ obtained using (21) for $k=m$, and

$F_{m}=\Bigg{(}\sum_{l=1}^{m}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top}\Bigg{)}^{-1},$

the sequence generated by the algorithm would be
$\big{\\{}\widehat{S}_{l}^{*}\big{\\}}_{l\geq m}$, i.e., for each $k>m$, the
algorithm would generate the exact solution to (20) as given in (21).

##  V Online Feedback Optimization

Now, on the one hand, by replacing $S_{k-1}$ by $\widehat{S}_{k-1}$ in (13),
we obtain that

$\displaystyle
q_{k}=\Big{[}q_{k-1}+\gamma_{k}\widehat{S}_{k-1}^{\top}\big{(}v^{*}-v_{k-1}\big{)}\Big{]}_{\underline{q}_{k}}^{\overline{q}_{k}},\quad
k\geq 1.$ (24)

On the other hand, by inspecting (21), one can see that an estimate of
$\widehat{S}_{k}$ can be computed if the matrix

$\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta q_{l}^{\top}$

is invertible; this can be ensured as mentioned earlier if the sequence
$\big{\\{}\Delta q_{l}\big{\\}}_{l=1}^{k}$ is persistently exciting. However,
it is not clear a priori that (24) will generate such persistently exciting
sequence. To address this issue, we modify (24) to add a mechanism to ensure
that indeed the sequence $\big{\\{}\Delta q_{l}\big{\\}}_{l=1}^{k}$ is
persistently exciting. First, by assuming that
$q_{k-1}\in[\underline{q}_{k-1},\overline{q}_{k-1}]$, we can rewrite (24) as
follows:

$\displaystyle
q_{k}=q_{k-1}+\Big{[}\gamma_{k}\widehat{S}_{k-1}^{\top}\big{(}v^{*}-v_{k-1}\big{)}\Big{]}_{\underline{\Delta
q}_{k}}^{\overline{\Delta q}_{k}},\quad k\geq 1,$ (25)

where $\underline{\Delta q}_{k}=\underline{q}_{k}-q_{k-1}$ and
$\overline{\Delta q}_{k}=\overline{q}_{k}-q_{k-1}$. Then, following the ideas
in [5], we will make three modifications to (25) as follows:

  1. M1.

The term inside the projection term on the right hand side of (25) is modified
as follows:

$\displaystyle\gamma_{k}\widehat{S}_{k-1}^{\top}\big{(}\Delta
v^{*}_{k}+cz_{k}\big{)},$ (26)

where $c$ is some constant, $\Delta v^{*}_{k}=v^{*}-v_{k-1}$, and

$z_{k}=\begin{cases}0,~{}\text{if}~{}\\{\Delta
v^{*}_{l}\\}_{l=1}^{k}~{}\text{is persistently exciting},\\\ \text{sampled
from}~{}(-\Delta v^{*}_{k})U(0,a_{1}),~{}\text{otherwise},\end{cases}$ (27)

with $a_{1}\in(0,1)$ and $U(x,y)$ denoting the continuous uniform distribution
over the interval $(x,y)$.

  2. M2.

The incremental lower and upper capacity limit in (25) are respectively set to

$\displaystyle\underline{q}_{k}+\underline{\eta}_{k},$
$\displaystyle\overline{q}_{k}-\overline{\eta}_{k},$ (28)

where $\underline{\eta}_{k}$ is sampled from $U(0,a_{2}|\underline{\Delta
q}_{k}|)$ and $\overline{\eta}_{k}$ is sampled from
$U(0,a_{2}|\overline{\Delta q}_{k}|)$, with $a_{2}\in(0,1)$.

  3. M3.

By using the two modifications above, define the following quantity:

$\displaystyle\widetilde{\Delta
q}_{k}=\Big{[}\gamma_{k}\widehat{S}_{k-1}^{\top}\big{(}\Delta
v^{*}_{k}+cz_{k}\big{)}\Big{]}_{\underline{\Delta
q}_{k}+\underline{\eta}_{k}}^{\overline{\Delta q}_{k}-\overline{\eta}_{k}}.$
(29)

Let $\mbox{Null}\big{(}\widehat{S}_{k}\big{)}$ denote the null space of the
matrix $\widehat{S}_{k}$. Then, if the sequence $\Big{\\{}\big{\\{}\Delta
q_{l}\big{\\}}_{l=1}^{k-1},\widetilde{\Delta q}_{k}\Big{\\}}$ is persistently
exciting, we update the value of $q_{k}$ as follows:

$\displaystyle q_{k}=q_{k-1}+\widetilde{\Delta q}_{k},\quad k\geq 1,$ (30)

otherwise we update its value as follows:

$\displaystyle q_{k}=q_{k-1}+\widetilde{\Delta
q}_{k}+\gamma_{k}\alpha_{k}\nu_{k},\quad k\geq 1,$ (31)

with $\alpha_{k}$ sampled from $U(-b_{k},b_{k})$, where $b_{k}\geq 0$ is
arbitrarily chosen so that

$\displaystyle\underline{\Delta q}_{k}$ $\displaystyle\leq\widetilde{\Delta
q}_{k}+\gamma_{k}b_{k}\nu_{k}\leq\overline{\Delta q}_{k},$ (32)
$\displaystyle\underline{\Delta q}_{k}$ $\displaystyle\leq\widetilde{\Delta
q}_{k}-\gamma_{k}\alpha_{k}\nu_{k}\leq\overline{\Delta q}_{k},$ (33)

for some arbitrarily chosen
$\nu_{k}\in\mbox{Null}\big{(}\widehat{S}_{k}\big{)}$; this ensures that
$\underline{q}_{k}\leq q_{k}\leq\overline{q}_{k}$.

The rationale behind the modifications above is to introduce excitation across
the entire space in which the $\Delta q_{k}$’s take values. To see this, first
consider the case when the capacity constraints are not active at instant $k$,
i.e., $\widetilde{\Delta
q}_{k}=\gamma_{k}\widehat{S}_{k-1}^{\top}\big{(}\Delta
v^{*}_{k}+cz_{k}\big{)}.$ Let $\mbox{Row}\big{(}\widehat{S}_{k}\big{)}$ denote
the row space of $\widehat{S}_{k}\in\mathds{R}^{n\times m}$. Then, the idea
behind adding the term $cz_{k}$ to $\Delta v_{k}^{*}$ in Modification M1 is to
introduce excitation in the subspace spanned by the columns of
$\widehat{S}^{\top}_{k},$ i.e., $\mbox{Row}\big{(}\widehat{S}_{k}\big{)}$,
whereas the idea behind adding the term $\gamma_{k}b_{k}\nu_{k}$ in
Modification M3 is to introduce excitation in
$\mbox{Null}\big{(}\widehat{S}_{k}\big{)}$, which is the orthogonal complement
of $\mbox{Row}\big{(}\widehat{S}_{k}\big{)}$. Then, since for any
$x\in\mathds{R}^{m}$ we have that $x=u+v,$ with
$u\in\mbox{Row}\big{(}\widehat{S}_{k}\big{)}$ and
$v\in\mbox{Null}\big{(}\widehat{S}_{k}\big{)}$, this mechanism allows us to
introduce excitation in separate parts of the whole space $\mathds{R}^{m}$ as
needed. The idea behind modifying the incremental capacity limits as described
in Modification M2 is as follows. If the capacity constraints are active at
instant $k$, by modifying the incremental capacity limits as in (28), we are
leaving headroom to add the term $\gamma_{k}b_{k}\nu_{k}$ if needed, while
still ensuring that the constraints in (32) – (33) imposed by the actual
incremental capacity limits are satisfied.

##  VI Convergence Analysis

Next, we provide the conditions under which the proposed online feedback
optimization controller in (29) – (33) converges. To this end, we first need
to rewrite the model in (1) as follows. Define $\delta v_{k}=v_{k}-v^{*}$,
$\delta q_{k}=q_{k}-q_{0}$, and $\delta w_{k}=w_{k}-w_{0}^{\circ}$ (note that
$w_{k}-w_{0}^{\circ}=w_{k}^{\circ}-w_{0}^{\circ}+\xi_{k}$); then, we have that

$\displaystyle v^{*}+\delta v_{k}$ $\displaystyle=h(q_{0}+\delta
q_{k},w_{0}^{\circ}+\delta w_{k}),\quad k\geq 0.$ (34)

Now by using the Taylor series expansion to expand the right hand side of (34)
about $(q_{0},w_{0}^{\circ})$, we obtain that

$\displaystyle v_{k}=S_{\varphi}q_{k}+S_{w}\xi_{k}+\mu_{k}+\eta_{k},$ (35)

where

$\displaystyle S_{\varphi}$ $\displaystyle=\frac{\partial
h(\varphi,w)}{\partial\varphi}\Bigg{|}_{\varphi=q_{0},~{}w=w_{0}^{\circ}}$
$\displaystyle S_{w}$ $\displaystyle=\frac{\partial h(\varphi,w)}{\partial
w}\Bigg{|}_{\varphi=q_{0},~{}w=w_{0}^{\circ}}$ $\displaystyle\mu_{k}$
$\displaystyle=v^{*}-S_{\varphi}q_{0}+S_{w}\big{(}w_{k}^{\circ}-w_{0}^{\circ}\big{)},$
(36)

with $\eta_{k}$ representing higher-order terms in $\delta q_{k}$ and $\delta
w_{k}$.

The next result, which is established using standard analysis techniques (see,
e.g., [17]), shows that the sequence $\\{q_{k}\\}_{k\geq 1}$ generated by (27)
– (33) converges almost surely to a solution of the optimization problem in
(7) provided that the sensitivity estimates,
$\big{\\{}\widehat{S}_{k}\big{\\}}_{k\geq 1}$, are unbiased. In doing so, we
set $c=0$ in (29). Also, we assume that (i) $\mu_{k}=\mu,~{}k\geq 0,$ where
$\mu$ is some constant, which is consistent with the assumption we made
earlier that $w_{k}^{\circ}$ slowly changes with $k$; and (ii)
$\eta_{k}=0,~{}k\geq 0$, which is reasonable since this term captures the
higher-order terms of the Taylor series expansion of $h(q,w)$ about
$(q_{0},w_{0}^{\circ})$. We plan to address these issues in future work.

###### Proposition 1

Consider the online feedback optimization controller in (27) – (33) with
$c=0$, and where the sequence $\big{\\{}\widehat{S}_{k}\big{\\}}_{k\geq 0}$
are the unbiased estimates of the sensitivity matrix $S_{\varphi}$, namely,
$\text{E}\Big{[}\widehat{S}_{k}-S_{\varphi}\,|\,\mathcal{F}_{k}\Big{]}=0$,
with $\mathcal{F}_{k}$ denoting the accumulated collection of states up to
instant $k$, $\big{\\{}(v_{l},q_{l})\big{\\}}_{l=0}^{k}$, and
$\text{E}\left[\big{\|}\widehat{S}_{k}-S_{\varphi}\big{\|}^{2}\right]\leq\sigma_{\varphi}$,
with $\sigma_{\varphi}$ denoting some positive constant. Let
$\mathcal{X}_{k}^{*}$ denote the set of optimal solutions of (7). We assume
that the following conditions hold:

(a) $\underline{q}_{k}$ and $\overline{q}_{k}$ are constant, $\mu_{k}$ is
constant, $\text{E}\big{[}\eta_{k}\big{]}=0$ and
$\text{E}\big{[}\|\eta_{k}\|^{2}_{2}\big{]}\leq\sigma$.

(b) $\gamma_{k}$ is a diminishing step size, namely,

$\sum_{k=1}^{\infty}\gamma_{k}=\infty,\quad\sum_{k=1}^{\infty}\gamma_{k}^{2}<\infty.$
(37)

Then, the sequence $\\{q_{k}\\}_{k\geq 1}$ converges almost surely to some
point in $\mathcal{X}_{k}^{*}$.

###### Proof:

Define

$\displaystyle f(\varphi)$
$\displaystyle:=\frac{1}{2}\text{E}_{\Xi_{k}}\left[\big{\|}v^{*}-h_{k}\big{(}\varphi,\Xi_{k}\big{)}\big{\|}_{2}^{2}\right],$
$\displaystyle g_{k}$ $\displaystyle:=\frac{\partial
f(\varphi)}{\partial\varphi}\Bigg{|}_{\varphi=q_{k}}=-S_{\varphi}^{\top}(v^{*}-S_{\varphi}q_{k}-\mu_{k}),$

$\delta_{k}:=(S_{\varphi}-\widehat{S}_{k})^{\top}\Delta v^{*}_{k+1}$; then, it
follows from (29) – (31) that

$\displaystyle q_{k+1}$
$\displaystyle=\Big{[}q_{k}+\gamma_{k+1}\widehat{S}_{k}^{\top}\Delta
v^{*}_{k+1}+\gamma_{k+1}\alpha_{k+1}\nu_{k+1}\Big{]}_{\underline{q}_{k+1}+\underline{\eta}_{k+1}}^{\overline{q}_{k+1}-\overline{\eta}_{k+1}}$
$\displaystyle=\Big{[}q_{k}+\gamma_{k+1}\big{(}S_{\varphi}+(\widehat{S}_{k}-S_{\varphi})\big{)}^{\top}\Delta
v^{*}_{k+1}$
$\displaystyle\quad+\gamma_{k+1}\alpha_{k+1}\nu_{k+1}\Big{]}_{\underline{q}_{k+1}+\underline{\eta}_{k+1}}^{\overline{q}_{k+1}-\overline{\eta}_{k+1}}$
$\displaystyle=\Big{[}q_{k}-\gamma_{k+1}({g_{k}}+S_{\varphi}^{\top}S_{w}{\xi_{k}}+S_{\varphi}^{\top}\eta_{k}$
$\displaystyle\quad+\delta_{k}-\alpha_{k+1}\nu_{k+1})\Big{]}_{\underline{q}_{k+1}+\underline{\eta}_{k+1}}^{\overline{q}_{k+1}-\overline{\eta}_{k+1}}.$
(38)

Define
$\zeta_{k}:=S_{\varphi}^{\top}S_{w}\xi_{k}+S_{\varphi}^{\top}\eta_{k}+\delta_{k}-\alpha_{k+1}\nu_{k+1}$;
then, by using the non-expansiveness property of the Euclidean projection
operator, we have that

$\displaystyle\|q_{k+1}-q^{*}\|^{2}$
$\displaystyle\leq\big{\|}q_{k}-q^{*}-\gamma_{k+1}\big{(}g_{k}+\zeta_{k}\big{)}\big{\|}^{2}$
$\displaystyle=\|q_{k}-q^{*}\|^{2}-2\gamma_{k+1}\big{(}g_{k}+\zeta_{k}\big{)}^{\top}(q_{k}-q^{*})$
$\displaystyle\quad+\gamma_{k+1}^{2}\|g_{k}+\zeta_{k}\|^{2}.$ (39)

It follows from the convexity property that

$\displaystyle-g_{k}^{\top}(q_{k}-q^{*})\leq f(q^{*})-f(q_{k}).$ (40)

By applying (40) to (39), we obtain that

$\displaystyle\|q_{k+1}-q^{*}\|^{2}$
$\displaystyle\leq\|q_{k}-q^{*}\|^{2}-2\gamma_{k+1}(f(q_{k})-f(q^{*}))$
$\displaystyle\quad-2\gamma_{k+1}\zeta_{k}^{\top}(q_{k}-q^{*})+\gamma_{k+1}^{2}\|g_{k}+\zeta_{k}\|^{2}.$
(41)

Define

$g^{*}:=\frac{\partial f(\varphi)}{\partial\varphi}\Bigg{|}_{\varphi=q^{*}};$

then, the following Lipschitz condition holds trivially for some $L>0$:

$\displaystyle\|g_{k}-g^{*}\|=\|S_{\varphi}^{\top}S_{\varphi}(q_{k}-q^{*})\|\leq
L\|(q_{k}-q^{*})\|.$ (42)

By using (42) and the fact that $2x^{\top}y\leq x^{2}+y^{2}$, for any
$x,y\in\mathds{R}^{n}$, we have that

$\displaystyle\|g_{k}+\zeta_{k}\|^{2}$
$\displaystyle=\big{\|}g_{k}-g^{*}+g^{*}+\zeta_{k}\big{\|}^{2}$
$\displaystyle\leq 2\|g_{k}-g^{*}\|^{2}+2\|g^{*}+\zeta_{k}\|^{2}$
$\displaystyle\leq 2L^{2}\|q_{k}-q^{*}\|^{2}+2\|g^{*}+\zeta_{k}\|^{2}.$ (43)

Since $\text{E}\big{[}\delta_{k}\,|\,\mathcal{F}_{k}\big{]}=0$,
$\text{E}\big{[}S_{\varphi}^{\top}S_{w}\xi_{k}\,|\,\mathcal{F}_{k}\big{]}=0$,
$\text{E}\big{[}S_{\varphi}^{\top}\eta_{k}\,|\,\mathcal{F}_{k}\big{]}=0$, and
$\text{E}\big{[}\alpha_{k+1}\nu_{k+1}\,|\,\mathcal{F}_{k}\big{]}=0$, we have
that

$\displaystyle\text{E}\big{[}\zeta_{k}\,|\,\mathcal{F}_{k}\big{]}$
$\displaystyle=\text{E}\big{[}S_{\varphi}^{\top}S_{w}\xi_{k}+S_{\varphi}^{\top}\eta_{k}+\delta_{k}-\alpha_{k+1}\nu_{k+1}\,|\,\mathcal{F}_{k}\big{]}$
$\displaystyle=0.$ (44)

Then, by taking an expectation of (41) and applying (43) and (44), we obtain
that

$\displaystyle\text{E}\left[\|q_{k+1}-q^{*}\|^{2}\,|\,\mathcal{F}_{k}\right]$
$\displaystyle\leq(1+2L^{2}\gamma_{k+1}^{2})\|q_{k}-q^{*}\|^{2}$
$\displaystyle\quad-2\gamma_{k+1}(f(q_{k})-f(q^{*}))$
$\displaystyle\quad-2\gamma_{k+1}\text{E}[\zeta_{k}\,|\,\mathcal{F}_{k}]^{\top}(q_{k}-q^{*})$
$\displaystyle\quad+2\gamma_{k+1}^{2}\text{E}\left[\|g^{*}+\zeta_{k}\|^{2}\,|\,\mathcal{F}_{k}\right])$
$\displaystyle=(1+2L^{2}\gamma_{k+1}^{2})\|q_{k}-q^{*}\|^{2}$
$\displaystyle\quad-2\gamma_{k+1}(f(q_{k})-f(q^{*}))$
$\displaystyle\quad+2\gamma_{k+1}^{2}\text{E}\left[\|g^{*}+\zeta_{k}\|^{2}\,|\,\mathcal{F}_{k}\right]).$
(45)

Further, it can be easily shown that

$\displaystyle\text{E}\left[\|g^{*}+\zeta_{k}\|^{2}\,|\,\mathcal{F}_{k}\right]<\infty.$
(46)

By applying the Robbins-Siegmund Theorem (see, e.g., [18, Lemma 11]) to (45),
we conclude that $q_{k}$ converges almost surely to some point in
$\mathcal{X}_{k}^{*}$. ∎

The next result provides the cost error bound at the time instant when the
sequence $\\{\Delta q_{k}\\}_{k\geq 1}$ is no longer persistently exciting.
This indeed will be the case if the rLSE in (22) is used to generate the
sequence $\\{\widehat{S}_{k}\\}_{k\geq 1}$ as $\widehat{S}_{k}$ will cease to
be an unbiased estimate of the sensitivity matrix, $S_{\varphi}$, as
$k\to\infty$.

###### Proposition 2

Suppose we have for some $\alpha$ and $T>0$ that

$\displaystyle\sum_{l=1}^{T}\lambda^{T-l}\Delta q_{l}\Delta
q_{l}^{\top}<\alpha I.$ (47)

Suppose that $q_{T}>\underline{q}_{T}+\underline{\eta}_{T}$ and
$q_{T}<\overline{q}_{T}-\overline{\eta}_{T}$. Then, the following relation
holds for any $q^{*}\in\mathcal{X}_{T}^{*}$:

$\displaystyle
f(q_{T})-f(q^{*})<(\overline{q}_{T}-\underline{q}_{T})\frac{\sqrt{m\alpha}+B}{\gamma_{T}},$
(48)

where $B$ is an upper bound for $\|\zeta_{l}\|$, $l\geq 1$.

###### Proof:

It follows from (47) that $\alpha I-\sum_{l=1}^{T}\lambda^{T-l}\Delta
q_{l}\Delta q_{l}^{\top}$ is a positive definite matrix. Hence, by [19,
Corollary 7.1.5], its trace is positive, namely,

$\displaystyle\operatorname{Tr}\left(\alpha
I-\sum_{l=1}^{T}\lambda^{T-l}\Delta q_{l}\Delta q_{l}^{\top}\right)>0,$ (49)

where $\operatorname{Tr}(A)=\sum_{i}a_{i,i}$ denotes the trace of matrix
$A=[a_{i,j}]$. Then, by using the assumption that
$q_{T}>\underline{q}_{T}+\underline{\eta}_{T}$ and
$q_{T}<\overline{q}_{T}-\overline{\eta}_{T}$, and applying the triangle
inequality, we have that

$\displaystyle\sqrt{m\alpha}$
$\displaystyle>\sqrt{\operatorname{Tr}\left(\sum_{l=1}^{T}\lambda^{T-l}\Delta
q_{l}\Delta q_{l}^{\top}\right)}$
$\displaystyle\geq\|\gamma_{T}(g_{T}+\zeta_{T})\|\geq\gamma_{T}(\|g_{T}\|-\|\zeta_{T}\|)$
$\displaystyle\geq\gamma_{T}\|g_{T}\|-B,$ (50)

Hence,

$\displaystyle\|g_{T}\|<\frac{\sqrt{m\alpha}+B}{\gamma_{T}}.$ (51)

By using (40) and the fact that $q_{T}\in[\underline{q}_{T},\overline{q}_{T}]$
and $q^{*}\in[\underline{q}_{T},\overline{q}_{T}]$, we obtain that

$\displaystyle\|g_{T}\|(\overline{q}_{T}-\underline{q}_{T})\geq
f(q_{T})-f(q^{*}).$ (52)

By combining (52) with (51) the result in (48) follows. ∎

Figure 1: IEEE 123-bus distribution test feeder. Figure 2: Trajectories of
active and reactive power at bus 1.

##  VII Numerical Results

Here, we present numerical simulation results illustrating the effectiveness
of our proposed online feedback optimization controller for voltage regulation
in power distribution systems. To this end, we employ the test system depicted
in Fig. 1, which is a modified version of the three-phase balanced IEEE
123-bus distribution test feeder presented in [4, 20], with reactive-power-
capable DERs added to the network at the following buses: $19$, $26$, $38$,
$49$, $56$, $64$, $78$, $89$, and $99$. At other buses in the distribution
network, random perturbations in the active and reactive power demand are
introduced every $100$ milliseconds (for example, see Fig. 2 for an
illustration of the active and reactive power demand at bus 1). In this case
we have that $n=122$ and $m=9$. Then, we set $F_{0}=I_{9}$ and
$\widehat{S}_{0}$ to the matrix that results from removing the last $113$
columns of the matrix $I_{122}$. Additionally, all the components of
$\underline{q}_{k}$ and $\overline{q}_{k}$ are set to $-0.5$ pu and $0.5$ pu,
respectively, for all $k$. We also set $a_{1}=0.8$, $a_{2}=0.8$,
$\lambda=0.995$ and $\gamma_{k}=0.95$, $k\in\\{1,\dots,100\\}$,
$\gamma_{k}=0.1$, $k\in\\{101,\dots,400\\}$.

Figure 3: Trajectory of voltage magnitude at bus 1 for c = 0. Figure 4:
Trajectory of voltage magnitude at bus 1 for c = 0.5. Figure 5: Trajectory of
MSE of voltage magnitude for c = 0.

Figures 3 and 4 show the trajectory followed by bus 1 voltage magnitude under
i) no control action, i.e., $q_{k}=q_{0}$ for all $k\geq 1$, and ii) the
action of the proposed controller for $c=0$ and $c=0.5$. We also display the
mean square error (MSE) of voltage magnitudes at all buses, i.e
$\frac{\sum_{i=1}^{n}(v_{i}^{*}-v_{i})^{2}}{n}$ in Fig. 5. All our numerical
simulations demonstrate that our proposed controller is extremely effective at
maintaining bus voltage magnitudes closed to their nominal value.

##  VIII Concluding Remarks

In this paper we have proposed a controller for voltage regulation in power
distribution networks using reactive-power capable DERs. The proposed
controller is based on a PSGD algorithm for solving online a sequence of
optimization problems, each of which capturing the objectives and constraints
of the voltage regulation problem at a particular time instant. By assuming
the cost functions of each of these problems vary slowly with time, we can
show that the PSGD-based algorithm acts as a feedback controller. In order to
execute the controller, it is necessary to know the sensitivities of changes
in bus voltage magnitudes with respect to changes in reactive power
injections; we assume this are not a priori known and use a rLSE to estimate
them.

We showcased the performance of the controller via numerical simulations on
the IEEE 123-bus system. In addition, under certain simplifying assumptions,
we showed that the the sequence of DER setpoints generated by the controller
converges almost surely to a solution of the aforementioned optimization
problem when the estimates used by the controller are unbiased. In this
regard, while the estimates generated by the rLSE will initially be unbiased,
in the limit this will no longer be the case; however, the simulation results
show that the controller still performs satisfactorily. We plan to investigate
this issue in future work and also investigate the effect on the convergence
analysis of relaxing the other assumptions we made in establishing the results
in this paper.

### Derivation of Sensitivity Estimator

In order to derive (22), consider the expression in (21):

$\displaystyle\widehat{S}_{k}$
$\displaystyle=\left(\sum_{l=1}^{k}\lambda^{k-l}\Delta v_{l}\Delta
q_{l}^{\top}\right)\Bigg{(}\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top}\Bigg{)}^{-1},$ (53)

Define

$\displaystyle F_{k}^{-1}=\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top};$ (54)

then, we have that

$\displaystyle F_{k}^{-1}$
$\displaystyle=\lambda\sum_{l=1}^{k-1}\lambda^{k-1-l}\Delta q_{l}\Delta
q_{l}^{\top}+\Delta q_{k}\Delta q_{k}^{\top}$ $\displaystyle=\lambda
F_{k-1}^{-1}+\Delta q_{k}\Delta q_{k}^{\top}.$ (55)

Then, by using the matrix inversion lemma (see, e.g., [19]), it follows that

$\displaystyle F_{k}$ $\displaystyle=\Big{(}\lambda F_{k-1}^{-1}+\Delta
q_{k}\Delta q_{k}^{\top}\Big{)}^{-1}$ $\displaystyle=\lambda^{-1}F_{k-1}$
$\displaystyle-\frac{\lambda^{-2}}{1+\lambda^{-1}\Delta
q_{k}^{\top}F_{k-1}\Delta q_{k}}F_{k-1}\Delta q_{k}\Delta
q^{\top}_{k}F_{k-1}.$ (56)

Now, we manipulate the expression in (53) as follows:

$\displaystyle\widehat{S}_{k}$
$\displaystyle=\left(\sum_{l=1}^{k}\lambda^{k-l}\Delta v_{l}\Delta
q_{l}^{\top}\right)\Bigg{(}\sum_{l=1}^{k}\lambda^{k-l}\Delta q_{l}\Delta
q_{l}^{\top}\Bigg{)}^{-1}$
$\displaystyle=\left(\lambda\underbrace{\sum_{l=1}^{k-1}\lambda^{k-1-l}\Delta
v_{l}\Delta q_{l}^{\top}}_{\widehat{S}_{k-1}F_{k-1}^{-1}}+\Delta v_{k}\Delta
q_{k}^{\top}\right)F_{k}$
$\displaystyle=\left(\lambda\widehat{S}_{k-1}F_{k-1}^{-1}+\Delta v_{k}\Delta
q_{k}^{\top}\right)F_{k}$
$\displaystyle=\left(\widehat{S}_{k-1}\Big{(}F_{k}^{-1}-\Delta q_{k}\Delta
q_{k}^{\top}\Big{)}+\Delta v_{k}\Delta q_{k}^{\top}\right)F_{k}$
$\displaystyle=\widehat{S}_{k-1}+\Big{(}\Delta v_{k}-\widehat{S}_{k-1}\Delta
q_{k}\Big{)}\Delta q_{k}^{\top}F_{k},$ (57)

where the next to the last equality follows from (55).

## References

  * [1] “IEEE guide for identifying and improving voltage quality in power systems,” _IEEE Std 1250-2018 (Revision of IEEE Std 1250-2011)_ , pp. 1–63, 2018.
  * [2] B. A. Robbins, C. N. Hadjicostis, and A. D. Domínguez-García, “A two-stage distributed architecture for voltage control in power distribution systems,” _IEEE Trans. Power Syst._ , vol. 28, no. 2, pp. 1470–1482, 2013\.
  * [3] H. Xu, A. D. Domínguez-García, V. V. Veeravalli, and P. W. Sauer, “Data-driven voltage regulation in radial power distribution systems,” _IEEE Trans. Power Syst._ , vol. 35, no. 3, pp. 2133–2143, 2020\.
  * [4] H. Xu, A. D. Domínguez-García, and P. W. Sauer, “Data-driven coordination of distributed energy resources for active power provision,” _IEEE Trans. Power Syst._ , vol. 34, no. 4, pp. 3047–3058, 2019.
  * [5] M. Zholbaryssov and A. D. Domínguez-García, “Safe data-driven secondary control of distributed energy resources,” _IEEE Trans. Power Syst._ , vol. 36, no. 6, pp. 5933–5943, 2021.
  * [6] H. Kushner and G. Yin, _Stochastic Approximation and Recursive Algorithms and Applications_ , ser. Stochastic Modelling and Applied Probability. Springer New York, 2003.
  * [7] P. Hidalgo-Gonzalez, R. Henriquez-Auba, D. S. Callaway, and C. J. Tomlin, “Frequency regulation using data-driven controllers in power grids with variable inertia due to renewable energy,” in _Proc. IEEE Power Energy Society General Meeting_ , 2019, pp. 1–5.
  * [8] S. S. Madani, C. Kammer, and A. Karimi, “Data-driven distributed combined primary and secondary control in microgrids,” _IEEE Trans. on Control Syst. Technol._ , pp. 1–8, 2020.
  * [9] X. Chen, J. I. Poveda, and N. Li, “Safe model-free optimal voltage control via continuous-time zeroth-order methods,” in _Proc. of the IEEE Conference on Decision and Control_ , 2021, pp. 4064–4070.
  * [10] Y. C. Chen, A. D. Domínguez-García, and P. W. Sauer, “Measurement-based estimation of linear sensitivity distribution factors and applications,” _IEEE Trans. Power Syst._ , vol. 29, no. 3, pp. 1372–1382, 2014.
  * [11] K. E. Van Horn, A. D. Domínguez-García, and P. W. Sauer, “Measurement-based real-time security-constrained economic dispatch,” _IEEE Trans. Power Syst._ , vol. 31, no. 5, pp. 3548–3560, 2016.
  * [12] X. Zhang, V. Purba, M. Hong, and S. Dhople, “A sum-of-squares optimization method for learning and controlling photovoltaic systems,” in _Proc. American Control Conference_ , 2020, pp. 2376–2381.
  * [13] M. Picallo, L. Ortmann, S. Bolognani, and F. Dörfler, “Adaptive real-time grid operation via online feedback optimization with sensitivity estimation,” _Electric Power Systems Research_ , vol. 212, 2022.
  * [14] A. Bergen and V. Vittal, _Power Systems Analysis_. Prentice Hall, 2000.
  * [15] C. Wilson, V. V. Veeravalli, and A. Nedich, “Adaptive sequential stochastic optimization,” _IEEE Trans. on Automat. Contr._ , vol. 64, no. 2, pp. 496–509, 2019.
  * [16] K. J. Åström and B. Wittenmark, _Adaptive Control: Second Edition_. Mineola, New York: Dover, 2008\.
  * [17] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro, “Robust stochastic approximation approach to stochastic programming,” _SIAM Journal on Optimization_ , vol. 19, no. 4, pp. 1574–1609, 2009.
  * [18] B. T. Polyak, _Introduction to Optimization_. Optimization Software, 1987.
  * [19] R. A. Horn and C. R. Johnson, _Matrix Analysis_ , 2nd ed. Cambridge University Press, 2013.
  * [20] _IEEE PES Test Feeder_ , IEEE, https://cmte.ieee.org/pes-testfeeders/resources/.


# Charge–Trapping–Induced Compensation of the Ferroelectric Polarization in
FTJs: optimal conditions for a synaptic device operation

R. Fontanini  M. Segatto  K. S. Nair  M. Holzer  F. Driussi  Member, IEEE  I.
Häusler  C. T. Koch  C. Dubourdieu  Senior Member, IEEE  V. Deshpande and D.
Esseni  Fellow, IEEE This work was supported by the European Union through the
BeFerroSynaptic project, Grant 871737.R. Fontanini
(fontanini.riccardo@spes.uniud.it), M. Segatto, F. Driussi and D. Esseni are
with the Polytechnic Department of Engineering and Architecture (DPIA),
University of Udine, via delle Scienze 206, Udine, Italy. K. S. Nair, M.
Holzer and C. Dubourdieu are with Helmholtz–Zentrum Berlin, Hahn–Meitner–Platz
1, Berlin, Germany and with Freie Universität Berlin, Physical Chemistry,
Arnimallee 22, 14195, Berlin, Germany.V. Deshpande is with Helmholtz–Zentrum
Berlin, Hahn–Meitner–Platz 1, Berlin, Germany.I. Häusler and C. T. Koch are
with Institut für Physik, Humboldt-Universität zu Berlin, Berlin, Germany.R.
Fontanini, M. Segatto, K. S. Nair and M. Holzer contributed equally to this
work.

###### Abstract

In this work, we present a clear evidence, based on numerical simulations and
experiments, that the polarization compensation due to trapped charge strongly
influences the ON/OFF ratio in Hf0.5Zr0.5O2 (HZO) based Ferroelectric Tunnel
Junctions (FTJs). Furthermore, we identify and explain compensation conditions
that enable an optimal operation of FTJs. Our results provide both key
physical insight and design guidelines for the operation of FTJs as
multi–level synaptic devices.

Charge trapping, Ferroelectric Tunneling Junction, Optimization, Synaptic
Device.

##  1 Introduction

The demand for energy efficient hardware platforms for artificial intelligence
has raised a growing interest for memristive devices capable of implementing
non volatile, multi–level conductance [1]. Prominent applications of such
synaptic devices are crossbar arrays for artificial deep neural networks [2],
and hybrid memristive–CMOS circuits for spike–based neuromorphic computing
[3].

The field driven polarization switching and the tunneling read current qualify
Ferroelectric Tunnel Junctions (FTJs) as high–impedance and low–energy
synaptic devices. In this respect, in order to overcome the challenge of
ultra–thin ferroelectric layers (1–3 nm) needed in Metal–Ferroelectric–Metal
(MFM) stacks, the FTJ design has been recently steered towards
Metal–Ferroelectric–Insulator–Metal (MFIM) devices (see Fig. 1), that are
particularly suited for the integration in the BEOL of CMOS circuits [5, 4].

Figure 1: a), b) Sketches of the Stack Anneal–Bottom (SAB) FTJ and of the Post
Anneal–Deposition (PAD) device, respectively. The main difference between SAB
and PAD is the annealing process undergone by the Al2O3 oxide. All the devices
have a top electrode area of 95$\times$95 $\mu$m2. c) Cross section TEM image
of one the FTJ stacks fabricated and characterized in this work (SAB).

However, the operation of MFIM FTJs inherently relies on delicate trade–offs.
In fact, the depolarization field enables a polarization dependent tunnelling
read current, but it also tends to destabilize the polarization in retention
and read mode. Moreover, in ferroelectric–dielectric (FE-DE) stacks the charge
injection and trapping substantially influence the stability of polarization
and its switching dynamics [6], as it has also been recognized and debated for
ferroelectric–based FETs [7, 8, 9]. Clearly a simulation driven optimization
of FTJs is necessary to meet the sometimes conflicting design targets. To this
purpose, it is indispensable that a credible modelling framework for the FTJs
is developed. This can only be achieved through an appropriate comparison and
validation with experiments.

In this work we present an unprecedented joint effort based on an in–house
developed numerical model and experiments to understand and to quantitatively
investigate the trade–offs implied in the operation of FTJs as synaptic
devices. The calibration of the simulator against experiments clarifies some
crucial aspects of the device operation related to the charge trapping inside
the MFIM stack. In particular, our results show that either a small or a very
large trapping induce, respectively, negligible or complete compensation of
the ferroelectric polarization, that hinder the FTJ operation in both cases.
In this respect, we report optimal compensation conditions that can be
exploited to optimize FTJs as synaptic devices with multiple conductance
levels.

##  2 Device fabrication and experiments

Two types of MFIM FTJ stacks were fabricated with the same nominal thickness
for FE and DE films, but through different process sequences. The sketch of
the devices is reported in Fig. 1. The first stack named ‘Stack Anneal–Bottom’
[SAB, Fig. 1a)] was fabricated as follows: a 30 nm TiN layer was sputtered on
$p$$+$$+$ Si substrate as the bottom electrode, then a 3 nm Al2O3 film
followed by a 10 nm Hf0.5Zr0.5O2 (HZO) layer were deposited by atomic layer
deposition at 250∘C. Finally, a 30 nm W layer was deposited as top electrode
by sputtering at room temperature. The entire stack underwent a
crystallization anneal by RTP at 400∘C for 120 s in N2 ambient. The top
electrodes were patterned by Ti/Pt lift–off followed by W etch in H2O2. A TEM
image of a SAB FTJ is shown in Fig. 1c).

The second FTJ named ‘Post–Anneal Deposition’ [PAD, Fig. 1b)] was fabricated
as follows. After the deposition of the 30 nm bottom TiN electrode onto the
$p$$+$$+$ Si substrate, a 10 nm HZO layer followed by a 30 nm TiN film were
first deposited. The obtained TiN/HZO/TiN stack was RTP annealed at 400∘C for
120 s in N2 ambient. After the anneal, the top TiN layer was completely etched
away and the 3 nm Al2O3 film was deposited at 250∘C. The deposition and
patterning of top W pads was then performed with the process described earlier
for the SAB device. Note that, in the PAD FTJ case, the HZO was not in contact
with the Al2O3 layer during annealing, while in the SAB stack the HZO and
Al2O3 layers were annealed together. This is expected to lead to different
HZO/Al2O3 interfaces in the two devices and, thus, to different chemical and
electrical boundary conditions for the switching and stabilization of the HZO
polarization. In particular, one may expect a significantly larger density of
traps at the annealed HZO/Al2O3 interface of the SAB stack compared to the PAD
case. TiN/HZO/TiN capacitors with a 10 nm HZO layer were also fabricated as a
reference.

After fabrication, SAB and PAD FTJs are woken–up through 2000 cycles of a
triangular waveform at 1 kHz with 4.5 V amplitude. The number of wake–up
cycles corresponds to what is needed to wake–up the reference TiN/HZO/TiN
capacitors. Then, polarization–voltage (P–V) and current–volage (I–V) loops
are measured by applying 1 ms long triangular pulses to the top electrode,
with the substrate at ground potential. To avoid the breakdown of devices, the
maximum voltage amplitudes used for SAB and PAD FTJs are 5 and 4.5 V,
respectively.

##  3 Modelling framework and calibration

For the FTJ simulations, we used the in–house numerical model whose details
can be found in [10, 11, 12]. In particular, the FE domain dynamics is solved
through multi–domain Landau, Ginzburg, Devonshire (LGD) equations, accounting
for the 3D electrostatics [10], whose only free parameters are the LGD
constants $\alpha_{i}$, $\beta_{i}$, $\gamma_{i}$, the resistivity $\rho$ of
the FE dynamics (setting the time scale), and the domain wall coupling
constant $k$. In all simulations, we used $\rho=100\ \Omega\cdot$m,
consistently with the values recently reported for HZO based capacitors [13,
14, 15], and set $k\simeq 0$ in accordance with recent first principle
calculations [16] and a detailed comparison between our model and transient
operation of ferroelectric based devices [15].

Figure 2: Examples for the band diagram across the FE–DE stack in read mode.
The polarization is taken as positive when it points toward the DE–FE
interface. $V_{T}$ is the external applied voltage at the metal contacting the
ferroelectric. $\Phi_{M}$ is the work function of the electrodes.
$\chi_{F}=2.4$ eV, $\chi_{D}=1.6$ eV are the electron affinity of the HZO
ferroelectric and Al2O3 dielectric, $E_{f,MD}$, $E_{f,MF}$ are the Fermi
levels of the MD and MF electrodes. (a) The voltage drop $V_{D}$ in the DE is
such that the HZO conduction band minimum, $E_{C,I}$, at the FE–DE interface
is smaller than the Fermi level at the MD electrode, thus enabling a
tunnelling injection, $I_{tunn}$, limited by the dielectric; (b) A smaller
$V_{D}$ compared to results in (a) is expected in case of a Poole–Frenkel
current, $I_{PF}$, in shallow HZO traps [17].

Concerning the read current, , it is assumed to be dominated by the tunnelling
across the Al2O3 layer, even if additional mechanisms assisted by defects are
also possible [17]. In this respect, Fig. 2 qualitatively shows that a similar
band bending across the Al2O3 layer is necessary for either a current
dominated by tunnelling ($I_{tunn}$, a) or by a Poole–Frenkel mechanism in
shallow HZO traps ($I_{PF}$, b).

The model also includes charge injection and trapping, which is described in
terms of areal densities at the FE–DE interface [12]. However, it is
acknowledged that the trap concentrations of this work should be considered as
effective areal densities including also a possible charge trapping inside the
FE and DE films. The model solves the trap occupation $f_{T}$ through first
order dynamic equations for each trap energy level $E_{T}$ at the FE–DE
surface, assuming a Fermi occupation functions in the metals MF and MD,
respectively at the ferroelectric and dielectric side (see Fig. 2) [12]. A
detailed balance condition ensures that the $f_{T}$ is given by the Fermi
function at equilibrium (i.e. in steady state at $V_{T}=0$ V). As discussed
above, the capture rates $c_{MD0}$, $c_{MF0}$ from the MD and the MF
electrodes, respectively, are attributed to tunnelling and depend on the area
($\sigma_{T}$) and energy ($\sigma_{E}$) cross sections of the traps. In each
simulation domain the charges $Q_{acc}$, $Q_{don}$ in acceptor and donor
traps, respectively, depend on $f_{T}$ and are proportional to the trap
densities $N_{acc}$, $N_{don}$. The trap dynamics are solved in all domains
self–consistently with the LGD equations for the ferroelectric behaviour [12].

Figure 3: Simulated and experimental charge versus electric field
characteristic for a TiN/HZO/TiN device. The mean values of the gaussian
distributed LGD constants are reported in the figure.

We calibrated the model against experiments and, in particular, Fig. 3
compares the experimental charge versus the electric field, $E_{F}$, curves
for the TiN$/$HZO$/$TiN structure and the simulations with the LGD model
described above. Such comparison allowed us to extract the LGD constants. In
particular, in order to reproduce the experiments well, we considered
domain–to–domain variations of the LGD parameters by assuming a normal
distribution of $\alpha$, $\beta$ and $\gamma$. The mean values of such
distributions are $\alpha=-5.8\cdot 10^{8}$ m/F, $\beta=2.7\cdot 10^{8}$
m5/F/C2 and $\gamma=1.3\cdot 10^{11}$ m9/F/C4, while the standard deviation of
each constant are $\sigma_{\alpha}=25\%$, $\sigma_{\beta}=5\%$ and
$\sigma_{\gamma}=8\%$. Throughout this work, we will use these nominal values,
if not otherwise stated.

In all simulations, the work functions of W and TiN electrodes were taken as
4.5 eV, while for the tunneling mass and the permittivity we assumed
$m_{FE}=0.38\cdot m_{0}$ and $\varepsilon_{FE}=34\cdot\varepsilon_{0}$ for the
HZO layer, $m_{DE}=0.15\cdot m_{0}$ and
$\varepsilon_{DE}=10\cdot\varepsilon_{0}$ for the Al2O3 dielectric, with
$m_{0}$ being the free electron mass and $\varepsilon_{0}$ the permittivity of
vacuum.

##  4 Evidence of trap contribution in FTJ operation

Figure 4 compares simulations and experiments concerning the
polarization–voltage ($P$–$V$) curves of the SAB devices. As it can be seen,
simulations neglecting charge injection and trapping (green line) exhibit a
much more stretched $P$–$V$ curves compared to experiments. These simulation
results are actually qualitatively consistent with experiments reported for a
HZO capacitor serially connected to a discrete ceramic capacitor ensuring zero
charge injection [6, 18]. Similarly, in Fig. 5, also current–voltage ($I$–$V$)
simulations hardly show any switching current when trapping is neglected.

A fairly good agreement with $P$–$V$ and $I$–$V$ experiments can instead be
achieved by including, in the simulations, an adequate charge density at the
FE–DE interface (see blue lines in Figs. 4 and 5). This is a clear evidence of
the contribution of the charge trapping in the operation of FTJs, that is in
line with recent results reported for FeFET devices [7, 8, 9]. It is worth
mentioning that in Figs. 4, 5 and 6 the LGD constants had to be changed
compared to Fig. 3, in order to increase the coercive field $E_{c}$, thus
improving the agreement with experiments. This increase of the apparent
$E_{c}$ in MFIM compared to MFM structures has been previously reported and it
has been ascribed to the division of external voltage between DE and FE layer.
The capacitance of the ferroelectric increases as the field approaches
$E_{c}$, causing an increased drop across the DE layer and thus higher
apparent $E_{c}$ than in MFM [19].

Figure 4: Polarization versus voltage characteristics for SAB FTJs (blue open
circles) measured by using the triangular waveform shown in the inset; the
Al2O3 thickness is $t_{D}$=$3$ nm. Corresponding simulations for no trapped
charge (green line) or by accounting for acceptor and donor type traps (blue
line) are shown.

Figure 5: Current versus voltage curves corresponding to the $P$–$V$ plots in Fig. 4. Corresponding trap densities $N_{acc}$ and $N_{don}$ are the same as in Fig. 4 and are summarized in Tab. 1. Table 1: Trapping cross sections $\sigma_{T,acc}$, $\sigma_{T,don}$, $\sigma_{E}$ used in simulations throughout the paper and trap densities used in the simulations of Fig. 4 (SAB) and Fig. 6 (PAD). | $N_{acc}$ | $N_{don}$ | $\sigma_{T,acc}$ | $\sigma_{T,don}$ | $\sigma_{E}$
---|---|---|---|---|---
| $[cm^{-2}eV^{-1}]$ | $[cm^{-2}eV^{-1}]$ | $[cm^{2}]$ | $[cm^{2}]$ | $[eV]$
SAB | $5.25\cdot 10^{13}$ | $2.5\cdot 10^{13}$ | $3.5\cdot 10^{-14}$ | $8\cdot 10^{-16}$ | $7\cdot 10^{-3}$
PAD | $1.12\cdot 10^{13}$ | $1\cdot 10^{13}$ | $7\cdot 10^{-14}$ | $7\cdot 10^{-14}$ | $7\cdot 10^{-3}$

Figure 6: Comparison similar that in Fig. 4 between simulations and
experiments for the $P$–$V$ curves for PAD devices. Corresponding trap
densities are reported in Tab. 1. Green line shows the no traps case, while
red line is the simulation with traps included.

Figure 6 reports an analysis similar to that in Fig. 4 for the PAD devices.
Even for the PAD FTJs, the simulations neglecting trapping (green line)
exhibit large discrepancies with experiments (symbols). However, the trap
density matching simulations with experiments is smaller than the
concentration used for the SAB devices (see Tab. 1). As already mentioned in
Sec. 2, this difference is expected based on the different annealing
conditions for the two devices. The difference in the simulated interfacial
charge in SAB and PAD samples is best illustrated by Fig. 7, reporting the
average polarization and trapped charge $Q_{int}=(Q_{acc}+Q_{don})$ along a
triangular $V_{T}$ waveform. Figure 7 shows that the polarization is partly
compensated by $Q_{int}$ in both devices, but the effect is more prominent in
SAB FTJs. Again, this is ascribed to the different annealing sequence between
the two devices and Fig. 7 provides an insight about the role of the trap
density at the HZO/Al2O3 interface in the FTJ operation.

Figure 7: Simulated polarization ($P$, dashed lines) and trapped interface
charge ($Q_{int}$, symbols) averaged across the device area for the SAB and
PAD simulations in Figs. 4 and 6; the $V_{T}$ waveform is shown as a yellow
solid line. The used interfacial trap densities are reported in Tab. 1.

Figure 8: Measured current density for a SAB device at a read voltage
$V_{R}=2$ V (black diamonds) versus the corresponding voltage drop, $V_{D}$,
across the 3 nm Al2O3 layer estimated by using simulations with different
density, $N_{acc}$ [cm-2eV-1], of acceptor traps. In these simulations the
density of donors does not affect the final result. Measured tunnelling
current in thin SiO2 layers (red symbols) and corresponding simulations with
the tunnelling model of this work are reported as a reference.

In read mode, the compensation of the positive polarization is expected to
reduce the band bending and thus the read current $I_{R}$ (see Fig. 2). In
this regard, Fig. 8 shows the measured ON current density of about 0.12 A/m2
(black diamonds) for SAB FTJs, together with the tunnelling current density in
metal–SiO2–metal (MIM) systems [20]; incidentally, our tunnelling model is in
good agreement with experiments in the MIM systems. For the FTJ device, the
voltage drop, $V_{D}$, across the Al2O3 layer is estimated by simulations and
for different $N_{acc}$. Figure 8 shows that the estimated $V_{D}$ for
$N_{acc}\simeq 5\cdot 10^{13}$ cm-2eV-1 (see SAB parameters in Tab. 1) is
clearly smaller than the value $V_{D}=(\Phi_{M}-\chi_{F})/q=2.1$ V required to
have a read current $I_{R}$ limited by the $I_{tunn}$ through the Al2O3 film
alone [see Fig. 2a)]. Consistently with this picture, the measured $I_{R}$ is
much smaller than the $I_{tunn}$ in MIM systems, and it may be limited by
Poole–Frenkel conduction [see $I_{PF}$ in Fig. 2b)].

##  5 Polarization–compensation–aware design of the FTJ

Figure 9: Minimum, $I_{R,min}$, and maximum, $I_{R,max}$, simulated read
current (right $y$–axis) versus the trap density at read voltage $V_{R}=2$ V.
$I_{R,min}$ and $I_{R,max}$ correspond to a SET pulse voltage of
$V_{SET,min}=2.5$ V and $V_{SET,max}=4.5$ V, respectively. The corresponding
tunnel electroresistance, TER $=(I_{R,max}-I_{R,min})/I_{R,min}$ is reported
(left $y$–axis).

In recent designs of neuromorphic processors, a reasonable target for the
minimum read current is set to about 100 pA [21, 22], which requires a read
current density of about $I_{R}/A\simeq 100$ A/m2 for an FTJ area of $A\approx
1$ $\mu$m2. Hereafter, we use our calibrated simulations to examine the
optimal design of FTJs in terms of $I_{R}$ and tunneling electroresistance TER
$=(I_{R,max}-I_{R,min})/I_{R,min}$. For the simulations in this section, the
LGD constants are those in Fig. 3, calibrated on the MFM device, in order to
precisely account for the characteristics of the integrated HZO layer.

Figure 10: Average polarization, $P$, interface charge, $Q_{int}$, and net or
compensated polarization $(P+Q_{int})$ versus the trap density in read
condition ($V_{R}=2$ V) and for different SET voltages. At low $N_{acc}$, the
$(P+Q_{int})$ degrades due to depolarization field and back–switching, while
at large $N_{acc}$, the $(P+Q_{int})$ decreases due to a large $P$
compensation.

Figure 11: Simulation results: (a) HZO conduction band minimum, $E_{C,I}$, in
read mode versus the trap density, where the Fermi level at the MD electrode
is $E_{F,MD}=0$ eV (see inset). (b) Ferroelectric domain back–switching, BSW,
defined as the percentage of domains having a negative $P$ in retention after
being set to a positive $P$ during the SET pulse.

Figure 9 reports the simulated minimum and maximum read current densities and
TER for a FTJ having an Al2O3 thickness scaled down to $t_{D}=2.5$ nm and for
different assumed trap densities $N_{acc}=N_{don}$. Interestingly the TER
exhibits an optimal value for $N_{acc}$ around $4\cdot 10^{13}$ cm-2eV-1,
which stems from the behavior of the net or compensated polarization
$(P+Q_{int})$ illustrated in Fig. 10. In fact, at low $N_{acc}$ values, the
compensation of the positive $P$ is too weak, so that the depolarization field
$E_{F}$ increases and destabilizes the polarization. In this respect, we
calculated the difference between the number of domains with positive
spontaneous polarization $P$ at the maximum value of the SET voltage and at
zero external bias (i.e. during retention) and then we divided this value by
the overall number of domains. This is indeed the fraction of the domains set
to positive $P$ during the SET pulse that back–switch to negative $P$ during
retention. We defined such value as back–-switching in Fig. 11b), which shows
that at low $N_{acc}$ a significant back–switching of the ferroelectric
domains occurs. The large back–-switching at low $N_{acc}$ deteriorates the
TER in Fig. 9.

At large $N_{acc}$, instead, there is an excessive compensation that reduces
$(P+Q_{int})$. This implies an increase of the HZO conduction band minimum,
$E_{C,I}$, at the FE–DE interface during the read mode, which is reported in
Fig. 11a). Such an $E_{C,I}$ raise leads to the large $I_{R}$ drop at high
$N_{acc}$ shown in Fig. 9. Moreover, a large compensation tends to merge the
states obtained by different $V_{SET}$ values, which seriously hinders a
multi–level read current and thus the operation as a synaptic device.

Figure 12: Simulation results for the FTJ of Fig. 9 under optimal compensation
conditions corresponding to $N_{acc}=4\cdot 10^{13}$ cm-2eV-1. The box plots
reveal overlapping distributions for the current levels for small area FTJs
with $n_{D}=100$ ferroelectric domains. Simulations have been performed for a
small device area of 2500 nm2 (corresponding to a domain size $d=5$ nm).

The results of Fig. 9 and the inspection of the internal quantities performed
in Figs. 10 and 11, clearly show that an optimal compensation condition exists
for the operation of FTJs. In this respect, as discussed above, both
negligible and excessive charge–trapping–induced compensation of the
ferroelectric polarization are detrimental for the operation of the FTJs under
study. Indeed, there exists a quite delicate balance between the suppression
of back–switching and TER values enabling a multi–level operation in read
mode. In Fig. 12, for a given FTJ structure, we leverage the optimal design
condition in Fig. 9 and show that eight current levels can be placed in the
memory window of the FTJ for the optimal compensation condition corresponding
to $N_{acc}=4\cdot 10^{13}$ cm-2eV-1, thus enabling a 3–bit synaptic weight
resolution.

##  6 Conclusion

We presented a joint effort between numerical modelling and experiments to
demonstrate the crucial role of polarization compensation due to trapped
charges in HZO based FTJs. We identified and physically explained the optimal
compensation condition for a robust operation of FTJs.

Even if a tailoring of trap densities may be challenging from a technological
standpoint, we experimentally demonstrated that appropriate processing steps
and annealing conditions can lead to different trap densities at the FE–DE
interface of fabricated FTJs. Furthermore, as reported in [23], the insertion
of an ultra–thin metal layer at the FE–DE interface may be an interesting
design option to control the compensation condition in FTJs.

## Acknowledgment

The authors acknowledge support of Dr. Adnan Hammud, FHI–Berlin, for TEM
preparation.

## References

  * [1] S. Yu, ”Neuro-inspired computing with emerging nonvolatile memorys,” Proceedings of the IEEE, vol. 106, no. 2, pp. 260-285, 2018, doi: 10.1109/JPROC.2018.2790840
  * [2] S. Ambrogio, P. Narayanan, H. Tsai, R. M. Shelby, I. Boybat, C. di Nolfo, S. Sidler, M. Giordano, M. Bodini, N. C. P. Farinha, B. Killeen, C. Cheng, Y. Jaoudi, G. W. Burr, ”Equivalent-accuracy accelerated neural-network training using analogue memory” Nature 558, 60–67, 2018, doi: 10.1038/s41586-018-0180-5
  * [3] E. Chicca and G. Indiveri , ”A recipe for creating ideal hybrid memristive-CMOS neuromorphic processing systems”, Applied Physics Letter 116, 120501, 2020, doi: 10.1063/1.5142089
  * [4] R.Fontanini, J.Barbot, M.Segatto, S.Lancaster, Q.Duong, F.Driussi, L.Grenouillet, F.Triozon, J.Coignus, T.Mikolajick, S.Slesazeck, and D. Esseni, “Polarization switching and interface charges in BEOL compatible ferroelectric tunnel junctions,” in 2021 IEEE 51st European Solid-State Device Research Conference (ESSDERC), 2021.
  * [5] S. Slesazeck, T. Mikolajick, ”Nanoscale resistive switching memory devices: a review” Nanotechnology, vol. 30, p. 352003, 2019, doi: 10.1088/1361-6528/ab2084
  * [6] H. W. Park, S. D. Hyun, I. S. Lee, S. H. Lee, Y. B. Lee, M. Oh, B. Y. Kim, S. G. Ryoo, and C. S. Hwang, “Polarizing and depolarizing charge injection through a thin dielectric layer in a ferroelectric–dielectric bilayer”, Nanoscale, vol. 13, pp. 2556–2572, 2021, doi: 10.1039/D0NR07597C
  * [7] K. Toprasertpong, M. Takenaka and S. Takagi, ”Direct Observation of Interface Charge Behaviors in FeFET by Quasi-Static Split C-V and Hall Techniques: Revealing FeFET Operation,” 2019 IEEE International Electron Devices Meeting (IEDM), 2019, pp. 23.7.1-23.7.4, doi: 10.1109/IEDM19573.2019.8993664.
  * [8] J. Li, Y. Qu, M. Si, X. Lyu, and P. D. Ye, “Multi-Probe characterization of ferroelectric/dielectric interface by C-V, P-V and conductance methods”, 2020 IEEE Symposium on VLSI Technology, pp. 1–2, 2020, doi: 10.1109/VLSITechnology18217.2020.9265069
  * [9] S. Deng, Z. Jiang, S. Dutta, H. Ye, W. Chakraborty, S. Kurinec, S. Datta, K. Ni, ”Examination of the Interplay Between Polarization Switching and Charge Trapping in Ferroelectric FET”, 2020 IEEE International Electron Devices Meeting (IEDM), 2020, pp. 4.4.1-4.4.4, doi: 10.1109/IEDM13553.2020.9371999.
  * [10] T. Rollo, F. Blanchini, G. Giordano, R. Specogna, and D. Esseni, “Stabilization of negative capacitance in ferroelectric capacitors with and without a metal interlayer” Nanoscale, vol. 12,, pp. 6121–6129, 2020, doi: 10.1039/C9NR09470A
  * [11] R. Fontanini, M. Massarotto, R. Specogna, F. Driussi, M. Loghi and D. Esseni, ”Modelling and design of FTJs as high reading–impedance synaptic devices,” 2021 5th IEEE Electron Devices Technology & Manufacturing Conference (EDTM), 2021, pp. 1–3, doi: 10.1109/EDTM50988.2021.9420863.
  * [12] R. Fontanini, M. Segatto, M. Massarotto, R. Specogna, F. Driussi, M. Loghi, and D. Esseni, “Modelling and design of FTJs as multi-level low energy memristors for neuromorphic computing”, IEEE Journal of the Electron Devices Society, pp. 1–1, 2021, doi: 10.1109/JEDS.2021.3120200.
  * [13] M. Kobayashi, N. Ueyama, K. Jang, and T. Hiramoto, “Experimental Study on Polarization-Limited Operation Speed of Negative Capacitance FET with Ferroelectric HfO2,” in 2016 IEEE International Electron Devices Meeting (IEDM), pp. 314–317, 2016, doi: 10.1109/IEDM.2016.7838402.
  * [14] T. Kim, J. A. del Alamo and D. A. Antoniadis, ”Dynamics of HfZrO2 Ferroelectric Structures: Experiments and Models,” 2020 IEEE International Electron Devices Meeting (IEDM), 2020, pp. 21.4.1-21.4.4, doi: 10.1109/IEDM13553.2020.9372013.
  * [15] D. Esseni, R. Fontanini, “Macroscopic and microscopic picture of negative capacitance operation in ferroelectric capacitors”, Nanoscale, vol. 13, NO. 21, pp. 9641–9650, 2021, doi: 10.1039/D0NR06886A.
  * [16] H. J. Lee, M. Lee, K. Lee, J. Jo, H. Yang, Y. Kim, S. C. Chae, U. Waghmare, J. H. Lee, ”Scale–free ferroelectricity induced by flat phonon bands in HfO2”, Science, volume 369, n. 6509, pp. 1343-1347, 2020, doi:10.1126/science.aba0067.
  * [17] H. Schroeder, “Poole–Frenkel–effect as dominating current mechanism in thin oxide films – An illusion?!,” Journal of Applied Physics, vol. 117, no. 21, p. 215103, 2015, doi: 10.1063/1.4921949.
  * [18] M. Hoffmann, M. Gui, S. Slesazeck, R. Fontanini, M. Segatto, D. Esseni, and T. Mikolajick, “Intrinsic nature of negative capacitance in multidomain Hf0.5Zr0.5O2–based ferroelectric/dielectric heterostructures,” Advanced Functional Materials, p. 2108494, 2021, doi: 10.1002/adfm.202108494.
  * [19] Jiang, L., Hu, W., Wei, Z., Xu, W. and Meng, H, ”High–Performance Organic Single-Crystal Transistors and Digital Inverters of an Anthracene Derivative”, Advanced Materials, vol. 21, pp. 3649–3653, 2009, doi: 10.1002/adma.200900503.
  * [20] N. Matsuo, T. Miura, A. Urakami, T. Miyoshi, ”Analysis of Direct Tunneling for Thin SiO2 Film”, Japanese Journal of Applied Physics, vol. 38, p. 3967, 1999, doi: 10.1143/JJAP.38.3967.
  * [21] N. Qiao, G. Indiveri, ”Scaling Mixed-Signal Neuromorphic Processors to 28 nm FD-SOI Technologies” IEEE Biomedical Circuits and Systems Conference (BioCAS),2016, doi: 10.1109/BioCAS.2016.7833854.
  * [22] M. Sharifshazileh, K. Burelo, T. Fedele, J. Sarnthein and G. Indiveri, ”A Neuromorphic Device for Detecting High-Frequency Oscillations in Human iEEG,” 2019 26th IEEE International Conference on Electronics, Circuits and Systems (ICECS), pp. 69–72, 2019, doi: 10.1109/ICECS46596.2019.8965192.
  * [23] D. J. Frank, P. M. Solomon, C. Dubourdieu, M. M. Frank, V. Narayanan and T. N. Theis, ”The Quantum Metal Ferroelectric Field-Effect Transistor,” IEEE Transactions on Electron Devices, vol. 61, no. 6, pp. 2145–2153, 2014, doi: 10.1109/TED.2014.2314652.