# Incremental Volt/Var Control for Distribution Networks via Chance-
Constrained Optimization

Antonin Colot∗, Elisabetta Perotti∗, Mevludin Glavic, and Emiliano Dall’Anese
∗ Equal contribution of the authors.This work was supported in part by the
National Science Foundation (NSF) awards 1941896. The work of E. Perotti was
supported by Schmidt Science Fellows, in partnership with the Rhodes Trust.
The work of A. Colot was supported by the Research Fellow Fellowship of the
F.R.S-FNRS.A. Colot and M. Glavic are with the Montefiore Institute,
University of Liège (email<EMAIL_ADDRESS>mevludin.glavic@uliege.be).E. Perotti and E. Dall’Anese are with the
Department of ECEE, University of Colorado Boulder (email:
<EMAIL_ADDRESS>emiliano.dallanese@colorado.edu).

###### Abstract

This paper considers an incremental Volt/Var control scheme for distribution
systems with high integration of inverter-interfaced distributed generation
(such as photovoltaic systems). The incremental Volt/Var controller is
implemented with the objective of minimizing reactive power usage while
maintaining voltages within safe limits sufficiently often. To this end, the
parameters of the incremental Volt/Var controller are obtained by solving a
chance-constrained optimization problem, where constraints are designed to
ensure that voltage violations do not occur more often than a pre-specified
probability. This approach leads to cost savings in a controlled, predictable
way, while still avoiding significant over- or under-voltage issues. The
proposed chance-constrained problem is solved using a successive convex
approximation method. Once the gains are broadcast to the inverters, no
additional communication is required since the controller is implemented
locally at the inverters. The proposed method is successfully tested on a low-
voltage 42-nodes network.

###### Index Terms:

Incremental Volt/Var, distributed energy resources, distribution networks,
chance-constrained optimization.

##  I Introduction

The increasing integration of distributed energy resources (DERs) is driving a
paradigm shift in electrical power networks, moving away from centralized
power plants and embracing decentralized energy systems. The growing share of
inverter-interfaced renewable energy resources (RESs) in the electricity
production mix, aimed to meet climate targets, poses a number of challenges on
all aspects of electrical distribution networks (DNs), ranging from the
planning to the real-time control [1, 2]. The goal of this paper is to address
the voltage regulation problem in DNs by leveraging inverters’ capabilities.

Literature Review. Traditionally, voltage regulation in DNs was achieved using
load tap changers or switchgears. However, the increasing variability
introduced by DERs can shorten their lifespan, and they may become
insufficient to resolve voltage issues [3]. On the other hand, continuous
improvements in power electronics converters create new possibilities for the
control of DERs [4], enabling new tools for voltage regulation. The real-time
control of DERs’ powers can be categorized based on the control architecture.
In a centralized approach, the system operator determines the DERs’ power
injections by solving a given instance of the optimal power flow (OPF)
problem. The central controller needs a precise knowledge of the network,
which may be hard to achieve in practice [5]. Even though recent feedback-
based methods [6, 7] do not require perfect knowledge of non-controllable
power injections, they rely on a communication infrastructure that is not
present in existing DNs. In this paper, we focus on decentralized strategies,
which offer simplicity and low implementation cost since they rely only on
local measurements to perform control actions [8, 9].

Volt/Var controllers determine reactive power injections from, e.g., a static
function of the local voltage measurements. The slopes of individual Volt/Var
curves can be tuned to achieve various objectives by solving appropriate
optimization problems [10]. Since these static feedback laws can lead to
oscillatory behaviors [11], incremental strategies, based on voltage
measurements and the past reactive power setpoint [12], are generally favored.
Chance-constrained approaches to design optimal rules for non-incremental
Volt/Var control are proposed in [13], [14]. These works consider a separate
set of gains for each DER. This may be challenging in practice, as it requires
knowledge of each DER location and an advanced communication infrastructure to
properly dispatch the gains.

We also mention some representative works in the context of data-driven
methods to learn a Volt/Var controller [15]. For learning-based strategies, it
is often difficult to analyze the closed-loop stability [16]. Some exceptions
are, e.g. [17, 18]; however, the controllers in these works do not minimize
reactive power usage, potentially causing additional losses in the DN. In
[19], closed-loop stability for a general class of Volt/Var controllers is
guaranteed but historical data are needed for training the learning-based
controllers. In the context of learning-based Volt/Var, however, it is
difficult to account for network topology changes. Indeed, changes in topology
may require collecting new data and re-training learning-based controllers,
which is time-consuming. In [20], authors propose a linear Volt/Var controller
and a methodology for adapting the controller to varying network topologies.
However, this requires solving multiple instances of the OPF.

(a)

(b)

Figure 1: (a) Proposed voltage regulation strategy. Based on the forecast
$\rho$ and the probability $\epsilon$ to violate voltage limits, the
controller gains ${\bf x}$ are computed offline and then dispatched to the
controllable DERs. (b) Illustrative explanation of the impact of the parameter
$\epsilon$ on the total amount of voltage violations. A smaller $\epsilon$
results in a more constrained optimization problem and therefore in fewer
voltage violations.

Statement of Contributions: We design a local incremental Volt/Var controller
scheme that allows to specify the maximum amount of voltage violation that is
tolerated. Fig. 1 illustrates the key steps of our approach. We formulate a
chance-constrained optimization problem to ensure that voltage violations
cannot occur more often than according to a given, prefixed probability. To
the authors’ knowledge, this is the first work that relies on chance-
constrained optimization to determine the gains of an incremental Volt/Var
controller. In particular, we consider a single set of gains suitable for all
DERs, thus keeping the offline communication infrastructure requirements low.

With respect to a standard Volt/Var controller, this approach aims to minimize
the generator’s reactive power injections. The controller gains are determined
in advance based on forecast data for power generation and loads. There is no
need of historical data of any type. Therefore, not only can our approach
easily handle planned topology changes but, even in case of unplanned changes,
it is still reliable since network conditions are taken into account as
feedback. To account for a variety of possible distributions of the
forecasting errors, a conservative approximation of the chance constraints is
derived [21]. Our approach uses an algorithm based on Successive Convex
Approximation (SCA) methods [22] to derive a convex, and therefore tractable
formulation of our chance-constrained problem.

Outline: In Section II, we define the DN model and formulate the problem.
Section III introduces our incremental local Volt/Var controller and analyses
its stability. Section IV presents a strategy to determine the controller
gains that optimize the performance of our controller. Section V provides
details regarding the implementation of the feedback controllers. Numerical
simulations in Section VI compare the proposed approach to benchmark methods.
Section VII concludes the paper.

##  II Problem formulation

###  II-A Power system model

Consider a balanced three-phase power distribution network111 _Notation.
Upper-case (lower-case boldface) letters are used for matrices (column
vectors); $(.)^{\top}$ denotes the transposition, $(.)^{*}$ denotes the
complex conjugate and $(.)^{-1}$ denotes the inverse matrix; $j$ the imaginary
unit and $|.|$ the absolute value of a number. If we consider a given vector
${\bf x}\in^{N}$, $\mathsf{diag}(\cdot)$ returns a $N\times N$ matrix with the
element of ${\bf x}$ in its diagonal. For vectors ${\bf x}\in\mathbb{R}^{n}$
and ${\bf u}\in\mathbb{R}^{m}$, $\|{\bf x}\|$ denotes the $\ell_{2}$-norm and
$({\bf x},{\bf u})\in\mathbb{R}^{n+m}$ denotes their vector concatenation. We
denote as $\mathbf{0}$ a vector with all zeros (the dimensions will be clear
from the context). We denote $\mathbb{R}_{>0},\mathbb{R}_{\geq 0}$ and
$\mathbb{R}_{\leq 0}$ the set of strictly positive, positive and negative real
numbers. $\mathds{1}$ denotes the $N\times 1$ column vector and $\mathds{I}$
the $N\times N$ identity matrix._ with $N+1$ nodes and hosting $G$ DERs. The
node 0 is taken to be the point of common coupling, while
$\mathcal{N}:=\\{1,...,N\\}$ is the set of remaining nodes. We consider a
phasor model of the single-phase equivalent DN where we define
$u_{k}=v_{k}e^{j\delta_{k}}\in\mathbb{C}$ the voltage at node $k$. The voltage
at node 0 is set to $u_{0}=V_{0}$. Using Ohm’s and Kirchhoff’s laws, one has
the usual relationship:

${\bf s}=\mathsf{diag}{({\bf u})}({\bf y}^{*}V_{0}+Y^{*}{\bf u}^{*}),$ (1)

where ${\bf u}:=\\{u_{k}\\}_{k\in{\cal N}}$ collect the voltages at every
node, and ${\bf s}={\bf p}_{\mathrm{av}}-{\bf p}_{\mathrm{l}}+j({\cal
A}\widetilde{\bf q}-{\bf q}_{\mathrm{l}})\in\mathbb{C}^{N}$ is the net power
injection at nodes $n\in{\cal N}$. The network bus admittance matrix is
partitioned such that ${\bf y}\in\mathbb{C}^{N},Y\in\mathbb{C}^{N\times N}$.
We define ${\bf p}_{\mathrm{av}},{\bf p}_{\mathrm{l}},{\bf
q}_{\mathrm{l}}\in\mathbb{R}^{N}$ as the vectors collecting the non-
controllable active power injection, the active and reactive power consumption
at nodes $n\in{\cal N}$, respectively. We define $\widetilde{\bf
q}\in\mathbb{R}^{G}$ as the vector collecting the controllable reactive power
of the $G$ DERs. The matrix ${\cal A}\in\mathbb{R}^{N\times G}$ maps the index
of a DER to the node where it is located.

Equation (1) describes the power flow equations. For a given vector of net
power injection ${\bf s}$, one can solve this non-linear system of equations
using numerical methods to find the vector of voltage phasors ${\bf u}$.
Notice that the system of equations (1) may have zero, one or many solutions.
For the rest of this paper, we make the following assumption.

######  Assumption 1 (Existence and uniqueness of a practical solution of the
power flow equations)

There exists at least one solution to the power flow equations (1). If
multiple solutions exist, we only consider the _practical_ solution, i.e., in
the neighborhood of the nominal voltage profile, we pick the high voltage, and
small line currents solution. $\Box$

Let us define ${\bf z}:=({\bf p}_{\text{av}},{\bf p}_{\text{l}},{\bf
q}_{\text{l}})$ as the concatenation of the non-controllable powers, the
algebraic map $H:=\mathbb{R}^{G+3N}\rightarrow\mathbb{R}^{N}$ and ${\bf
v}:=\\{v_{k}\\}_{k\in{\cal N}}$ the vector collecting the voltage magnitudes
at every node. For convenience, we denote ${\bf v}=H(\widetilde{\bf q},{\bf
z})$ where $H$ relates the net power injections to the _practical_ solution of
the power flow equations defined in (1). Although the analytical form of the
map $H$ is not known, its existence and uniqueness have been discussed for
balanced [23, 24] and multi-phase DNs [25].

###  II-B Problem setup

The deployment of inverter-interfaced generation in DNs might induce voltage
quality issues. However, DERs can also be used to provide voltage regulation
services via, e.g., reactive power compensation. Ideally, one wants to
minimize the total reactive power usage while maintaining voltages inside a
given feasible set. We formulate an optimal reactive power flow problem as
follows:

$\displaystyle(\text{P0}(t))\quad\min_{\widetilde{\bf q}\in\mathcal{Q}}$
$\displaystyle\quad f(\widetilde{\bf q})$ (2) s.t. $\displaystyle\quad
H(\widetilde{\bf q},{\bf z}(t))\in\mathcal{V},$

where $\mathcal{Q}\subset\mathbb{R}^{G}$ is the set of feasible values for the
reactive powers, $\mathcal{V}\subset\mathbb{R}^{N}$ is the set of feasible
values for voltage magnitudes, and $f:\mathbb{R}^{G}\rightarrow\mathbb{R}$ is
a differentiable cost function. Notice that (P0($t$)) is a time-varying
optimization problem. Indeed, the non-controllable power injections ${\bf
z}(t)$ vary with time as they depend on users’ habits and weather conditions.
Therefore, also the optimal reactive power injections are time dependent;
their time-scale is determined by the variability of the non-controllable
power injections and is usually within seconds [26]. However, collecting data,
solving problem (P0($t$)) and broadcasting the setpoints to the inverters
every second is challenging because of the non-linear nature of the power flow
equations (1) represented by $H$, and the communication burden associated with
large DNs. This paper aims to solve the following problem.

###### Problem 1

Design feedback controllers to approximate the solution of (P0($t$)) with
limited computational resources and in a decentralized fashion, i.e., each
controller uses only local voltage measurements to compute its reactive power
output and the operation does not necessitate continuous communications with
neighbours or a centralized entity. $\Box$

To outline the proposed framework, we consider the case where there is one DER
per node in the network and the map $H$ does not change over time.
Additionally, we make the following assumption.

######  Assumption 2 (Feasibility)

For any ${\bf q}\in\mathcal{Q}$, there exists a solution such that ${\bf
v}(t)\in\mathcal{V}$. $\Box$

This assumption ensures that there is enough reactive power reserve to
regulate the voltages. We will review this assumption and our setup choices
later, and discuss how they could be relaxed.

##  III Design of feedback controllers

Within the next section, we identify three different time scales. In order
from the shortest to the longest, they correspond to the time scales that
characterize the following events: i) controller law updates, ii) forecast
updates and iii) controller gains updates.

#### Incremental Volt/Var control

To begin with, we consider only the shortest time scale, which is related to
the controller law updates. Let us discretize the temporal domain as
$t=k\tau$, where $k\in\mathbb{N}_{+}$ and $\tau\in\mathbb{R}_{>0}$ is a given
time interval, small enough to resolve variations in the time-dependent
disturbances, i.e., less than a second. We consider the following feedback
controller:

${\bf q}_{k+1}={\bf
q}_{k}+\eta(\mathds{1}-{\mbox{\boldmath$\nu$}}_{k})-(1-\eta)\alpha{\bf
q}_{k},$ (3)

where $\alpha\in\mathbb{R}_{\geq 0}$, $\eta\in[0,1]$, and
${\mbox{\boldmath$\nu$}}_{k}=X{\bf q}_{k}+{\mbox{\boldmath$\rho$}}_{k}$ with
$X\in\mathbb{R}^{N\times N}$ is a linear approximation of the power flow
equations and ${\mbox{\boldmath$\rho$}}_{k}=H(\mathbf{0},{\bf z}_{k})$ denotes
the voltage profile obtained by setting the controllable reactive powers to 0.
The linearized power flow equations can be derived from the bus injection
model [26], or from the branch flow model [27]. For the rest of this paper, we
consider the linearized power flow equations based on the branch flow model
[28], since it guarantees $X$ being positive definite [29]. This approximation
has been used in [27, 16], and its quality has been discussed in [29]. Notice
that ${\mbox{\boldmath$\rho$}}_{k}=H(\mathbf{0},{\bf z}_{k})$ is derived from
the true power flow equations. Substituting this approximation in (3), the
controller can be written as:

${\bf q}_{k+1}=A(\eta,\alpha){\bf q}_{k}+B(\eta,{\mbox{\boldmath$\rho$}}),$
(4)

where $A(\eta,\alpha):=(1-(1-\eta)\alpha)\mathds{I}-\eta X$ and
$B(\eta,{\mbox{\boldmath$\rho$}}):=\eta(\mathds{1}-{\mbox{\boldmath$\rho$}})$.

#### Existence and uniqueness of the equilibrium

Denoting ${\mbox{\boldmath$\rho$}}={\mbox{\boldmath$\rho$}}_{k}$ for a given
$k$, the equilibrium for (4) is defined as:

$\displaystyle{\bf q}^{*}$ $\displaystyle=\left[\eta
X+(1-\eta)\alpha\mathds{I}\right]^{-1}\eta(\mathds{1}-{\mbox{\boldmath$\rho$}})$
(5) $\displaystyle{\mbox{\boldmath$\nu$}}^{*}$ $\displaystyle=X\left[\eta
X+(1-\eta)\alpha\mathds{I}\right]^{-1}\eta(\mathds{1}-{\mbox{\boldmath$\rho$}})+{\mbox{\boldmath$\rho$}}.$

Since $X$ is positive definite and $\alpha\in\mathbb{R}_{\geq 0}$,
$\eta\in[0,1]$, the matrix $\eta X+(1-\eta)\alpha\mathds{I}$ in (5) is always
invertible and the equilibrium is unique. One can check that for $\eta=0$ and
$\alpha>0$, ${\bf q}^{*}=\mathbf{0}$ and
${\mbox{\boldmath$\nu$}}^{*}={\mbox{\boldmath$\rho$}}$, while for $\eta>0$ and
$\alpha=0$ or $\eta=1$, ${\bf
q}^{*}=X^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}})$ and
${\mbox{\boldmath$\nu$}}^{*}=\mathds{1}$.

Increasing the gain $\alpha$ decreases the use of reactive power, while
increasing the gain $\eta$ steers the voltage magnitudes to the nominal
voltage profile.

#### Stability analysis

The controller (4) is asymptotically stable if and only if
$\rho(A(\eta,\alpha))<1$, where $\rho(\cdot)$ denotes the spectral radius.
This condition is verified if:

$\mathbf{0}<(1-\eta)\alpha\mathds{1}+\eta{\mbox{\boldmath$\lambda$}}_{X}<2\mathds{1},$
(6)

where ${\mbox{\boldmath$\lambda$}}_{X}\in\mathbb{R}^{N}$ is the vector
containing the eigenvalues of the matrix $X$. Moreover the matrix $X$ is
positive definite by construction, and since $\eta,\alpha$ satisfy
$\alpha\in\mathbb{R}_{\geq 0}$, $\eta\in[0,1]$, we always have
$(1-\eta)\alpha\mathds{1}+\eta{\mbox{\boldmath$\lambda$}}_{X}\geq\mathbf{0}$.
The equality holds only if $\eta=\alpha=0$, which guarantees a stable
controller since ${\bf q}_{k+1}={\bf q}_{k}$.

#### Multi-phase unbalanced distribution networks

DNs are often highly unbalanced, and may have different connection
configurations (wye-connected or delta-connected). In order to study the
stability properties of the controller (3) for unbalanced DNs, one would need
to use a different approximation of the power flow equations, e.g., [25, 30].
In our formulation, we leverage the positive definite property of the matrix
$X$ to show the uniqueness and existence of the equilibrium, as well as the
stability of the controller. To adapt this framework to unbalanced DNs, one
would need to study the properties of the new matrix $X$ to characterize the
equilibrium and stability properties of the controller (3). However, as
mentioned in [30], due to the structure of the distribution lines, $X$ is
often positive definite for unbalanced networks. In such cases, our
methodology can be readily applied to multi-phase networks.

##  IV Design of the controller gains

The performance of the controller (4) depends on the choice of $\eta$ and
$\alpha$. In the following section, we introduce an optimization-based method
to design the gains.

###  IV-A Time-varying formulation

Given a matrix $X$, a time-varying vector
${\mbox{\boldmath$\rho$}}_{k}=H(\mathbf{0},{\bf z}_{k})$, and feasible sets
$\mathcal{Q}$ and $\mathcal{V}$, we formulate the following problem at time
$k\tau$:

$\displaystyle(\text{P1}_{k})\quad\min_{\alpha\in\mathbb{R}_{\geq
0},\eta\in[0,1]}$ $\displaystyle\quad\|{\bf q}_{k}(\alpha,\eta)\|^{2}$ (7a)
s.t. $\displaystyle\quad{\bf q}_{k}(\alpha,\eta)\in\mathcal{Q}$ (7b)
$\displaystyle\quad X{\bf
q}_{k}(\alpha,\eta)+{\mbox{\boldmath$\rho$}}_{k}\in{\cal V}$ (7c)
$\displaystyle\quad(1-\eta)\alpha\mathds{1}+\eta{\mbox{\boldmath$\lambda$}}_{X}<2\mathds{1}$
(7d)

where ${\bf q}_{k}(\alpha,\eta)=[\eta
X+(1-\eta)\alpha\mathbb{I}]^{-1}\eta(\mathds{1}-{\mbox{\boldmath$\rho$}}_{k})$
is a non-linear function of $\eta$ and $\alpha$. For a given time $k\tau$, the
goal is to minimize the reactive power usage $\|{\bf
q}_{k}(\alpha,\eta)\|^{2}$, while satisfying operational constraints, by
appropriately selecting the gains $\eta$ and $\alpha$.

An alternative formulation of problem (P1k) can be derived using an
appropriate change of variable. We introduce the optimization variable ${\bf
x}=[\frac{\alpha}{\eta},-\alpha]^{\top}\in\mathbb{R}_{\geq
0}\times\mathbb{R}_{\leq 0}$, rewrite ${\bf q}_{k}({\bf
x})=\left[X+\mathbf{1}^{\top}{\bf
x}\mathds{I}\right]^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}}_{k})$, specify
$\mathcal{Q}$ and ${\cal V}$ in terms of box constraints, and reformulate the
problem ($\text{P1}_{k}$) as

$\displaystyle(\text{P2}_{k})\quad\min_{{\bf x}}$ $\displaystyle\quad
h_{0}^{k}({\bf x})$ (8) s.t. $\displaystyle\quad h_{i}^{k}({\bf x})\leq
0\quad\forall i\in\\{1,...,8\\}$

where

$\displaystyle h_{0}^{k}({\bf x})=\|{\bf q}_{k}({\bf x})\|^{2},\quad
h_{1}^{k}({\bf x})={\bf q}_{k}({\bf x})-{\bf q}_{\mathrm{max}},$ (9)
$\displaystyle h_{2}^{k}({\bf x})=-{\bf q}_{k}({\bf x})+{\bf
q}_{\mathrm{min}},$ $\displaystyle h_{3}^{k}({\bf x})=X{\bf q}_{k}({\bf
x})+{\mbox{\boldmath$\rho$}}_{k}-{\mbox{\boldmath$\nu$}}_{\mathrm{max}},$
$\displaystyle h_{4}^{k}({\bf x})=-X{\bf q}_{k}({\bf
x})-{\mbox{\boldmath$\rho$}}_{k}+{\mbox{\boldmath$\nu$}}_{\mathrm{min}},$
$\displaystyle h_{5}({\bf x})=\left(\mathbf{1}^{\top}{\bf
x}-2\right)\mathds{1}+{\mbox{\boldmath$\lambda$}}_{X},$ $\displaystyle
h_{6}({\bf x})=-\mathbf{1}^{\top}{\bf x},\quad h_{7}({\bf x})=-x_{1},\quad
h_{8}({\bf x})=x_{2},$

with $x_{1},\,x_{2}$ scalar components of ${\bf x}$ and
$\mathbf{1}=[1,1]^{\top}\in\mathbb{R}^{2}$. The objective function $h_{0}$ is
a scalar function, as well as $h_{6},\,h_{7}$ and $h_{8}$. All other
constraint functions $h_{i}$ with $i=1,...,5$ are vector-valued. Function
$h_{5}({\bf x})$ gives a tighter bound on (7d) with the equality reached for
$\eta=1$, such that the controller gains derived from $(\text{P2}_{k})$ ensure
asymptotic stability of the controller defined in (3) as long as $\eta<1$. The
conditions $\eta<1$ is always verified as $\eta=1$ implies
${\mbox{\boldmath$\nu$}}^{*}=\mathds{1}$ which leads to a suboptimal solution
(unless
${\mbox{\boldmath$\nu$}}_{\mathrm{min}}={\mbox{\boldmath$\nu$}}_{\mathrm{max}}=\mathds{1}$).

As explained in Section II-B, collecting measurements at every nodes, solving
$(\text{P2}_{k})$ and then dispatching the controller gains to the
controllable DERs in real-time is unfeasible because of the communication
burden and the computational time required to solve $(\text{P2}_{k})$. One
could envision solving $(\text{P2}_{k})$ offline for every time $k\tau$ using
forecasts of ${\mbox{\boldmath$\rho$}}_{k}$. However, it is not realistic to
have such frequent forecast updates, since $\tau$ should be sufficiently small
to cope with the DN dynamics. Furthermore, ${\mbox{\boldmath$\rho$}}_{k}$ is
affected by large uncertainties as it inherits them from ${\bf z}_{k}$ through
${\mbox{\boldmath$\rho$}}_{k}=H(\mathbf{0},{\bf z}_{k})$. Moreover, we would
like to find optimal controller gains $\eta$ and $\alpha$ over a longer time
period, to avoid broadcasting new values at every time $k\tau$, or to avoid
storing a large number of gains in each controller. In the next section we
address this issue by reformulating our problem in a chance-constrained
fashion.

$k\tau$$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$\Delta\tau$$\Delta\tau$$\Delta\tau$$T\tau$

Figure 2: Comparison between different time scales of the problem, assuming
$\tau=1s$.

###  IV-B Chance-constrained formulation

Given that our forecast ${\mbox{\boldmath$\rho$}}_{k}=H(\mathbf{0},{\bf
z}_{k})$ is subject to uncertainty in the vector ${\bf z}_{k}$, we implement
probabilistic constraints in our optimization problem in order to enforce
voltage regulation with prescribed probability. The available DERs’ powers at
time $k\tau$ are modeled by ${\bf p}_{\text{av},k}=\bar{{\bf
p}}_{\text{av},k}+{\mbox{\boldmath$\delta$}}_{\text{av},k}$, while the active
and reactive loads are expressed as ${\bf p}_{\text{l},k}=\bar{{\bf
p}}_{\text{l},k}+{\mbox{\boldmath$\delta$}}_{\text{pl},k}$ and ${\bf
q}_{\text{l},k}=\bar{{\bf
q}}_{\text{l},k}+{\mbox{\boldmath$\delta$}}_{\text{ql},k}$, respectively.
Writing
${\mbox{\boldmath$\delta$}}_{k}:=({\mbox{\boldmath$\delta$}}_{\text{av},k},{\mbox{\boldmath$\delta$}}_{\text{pl},k},{\mbox{\boldmath$\delta$}}_{\text{ql},k})$,
and $\bar{\bf z}_{k}:=(\bar{{\bf p}}_{\text{av},k},\bar{{\bf
p}}_{\text{l},k},\bar{{\bf q}}_{\text{l},k})$, we have
${\mbox{\boldmath$\rho$}}_{k}=H(\mathbf{0},\bar{\bf
z}_{k}+{\mbox{\boldmath$\delta$}}_{k})$ where ${\mbox{\boldmath$\delta$}}_{k}$
follows a given distribution function.

It is reasonable to assume that a new forecast for
${\mbox{\boldmath$\rho$}}_{k}$ will be available after a certain time interval
$\Delta\tau$ with $\Delta\in\mathbb{N}_{+}$ (for example, $\Delta$ could be
such that $\Delta\tau=15$ minutes). Concretely, this means that within a time
window of duration $\Delta\tau$, ${\mbox{\boldmath$\rho$}}_{k}$ will be the
same, regardless of the value of the index $k$. Therefore, we will drop the
index $k$ in the following, and consider ${\mbox{\boldmath$\rho$}}^{m}$
instead to underline that the forecast ${\mbox{\boldmath$\rho$}}^{m}$ will be
updated at each time $t=m\Delta\tau$ with $m\in\mathbb{N}_{+}$. This
introduces a longer time scale, whose magnitude is related to how often
forecast updates occur.

Ultimately, our goal is to determine controller gains to be deployed over an
even longer time interval $T\tau$ where $T=b\Delta$ with $b\in\mathbb{N}_{+}$,
e.g., $b$ is such that $T\tau=1$ hour. Fig. 2 provides a visual representation
of the relationship between the three different time scales relevant to this
problem.

We consider an extension of $(\text{P2}_{k})$ as a multi-period optimization
problem:

$\displaystyle(\text{P3})\quad\min_{{\bf x}}$
$\displaystyle\quad\sum_{m=1}^{b}\mathbb{E}\\{h_{0}({\bf
x};{\mbox{\boldmath$\rho$}}^{m})\\}$ (10) s.t.
$\displaystyle\quad\mathrm{Pr}\\{h_{i,n}({\bf
x};{\mbox{\boldmath$\rho$}}^{m})\leq 0\\}\geq 1-\epsilon_{i}$
$\displaystyle\qquad\forall i\in\\{1,2,3,4\\},n\in{\cal N},\,m\in\\{1,..,b\\}$
$\displaystyle\quad h_{i}({\bf x})\leq 0\quad\forall i\in\\{5,6,7,8\\}$

where we sum over $b$ intervals of magnitude $\Delta$, corresponding in total
to a period of time $T\tau$. Notice that for $b=1$, we recover a single-
interval formulation. $\mathrm{Pr}\\{A\\}$ denotes probability of an event $A$
to happen, meaning that in problem (P3) the constraints $h_{i}({\bf
x};{\mbox{\boldmath$\rho$}}^{m})\leq 0$ for $i=1,2,3,4$ are satisfied with a
probability $1-\epsilon_{i}$, where $\epsilon_{i}\in(0,1)$. Solving problem
(P3) for the optimization variable ${\bf x}$, we can retrieve the values of
the gains $\alpha$ and $\eta$ to be deployed during a time interval of length
$T\tau$. When following this approach we lose optimality in exchange for
convenience: (P3) can be solved offline using a coarser forecast and
controller gains are designed to cover wider time windows, which is a great
advantage from a practical point of view. Finally, notice that the objective
function $h_{0}({\bf x};{\mbox{\boldmath$\rho$}}^{m})$ also depends on
${\mbox{\boldmath$\rho$}}^{m}$, thus we minimize its expected value.

We seek a tractable approximation for the chance constraints in (10) since we
do not know the probability distribution function of
${\mbox{\boldmath$\delta$}}_{k}$, neither the map $H$. The chance constraints
to be approximated are of the form $\mathrm{Pr}\\{h({\bf
x};{\mbox{\boldmath$\rho$}})\leq 0\\}\geq 1-\epsilon$, where the function
$h({\bf x},{\mbox{\boldmath$\rho$}})$ depends on the optimization variable
${\bf x}$ and the random vector $\rho$. Consider $\psi(x)=[1+x]_{+}$, where
$[x]_{+}:=\max\\{x,0\\}$, a so-called generating function
$\psi:\mathbb{R}\rightarrow\mathbb{R}$ nonnegative, nondecreasing, and convex
that satisfies the conditions $\psi(x)>\psi(0)\,\forall x>0$ and $\psi(0)=1$.
Given a positive scalar $z>0$, we have that the following bound holds for all
$z>0$ and ${\bf x}$ [21]:

$\inf_{z\in\mathbb{R}}\\{\mathbb{E}_{\mbox{\boldmath$\rho$}}\\{[h({\bf
x},{\mbox{\boldmath$\rho$}})+z]_{+}\\}-z\epsilon\\}\leq 0.$ (11)

Each probabilistic constraint in (10) will be replaced by the approximation
(11):

$\mathbb{E}_{{\mbox{\boldmath$\rho$}}^{m}}\\{[h_{i,n}({\bf
x};{\mbox{\boldmath$\rho$}}^{m})+u_{i,n}^{m}]_{+}\\}-u_{i,n}^{m}\epsilon\leq
0,$ (12)

where $u_{i,n}^{m}$ are real and positive auxiliary optimization variables.
Moreover, since the max operator $[.]_{+}$ is not differentiable, we replace
it with a smooth approximation and define:

$\displaystyle g_{i,n}({\bf x},u_{i,n};{\mbox{\boldmath$\rho$}}^{m})=$ (13)
$\displaystyle\quad\frac{1}{2}\left(h_{i,n}({\bf
x})+u_{i,n}+\sqrt{\xi^{2}+(h_{i,n}({\bf
x})+u_{i,n})^{2}}\right)-u_{i,n}\epsilon$ $\displaystyle\qquad\forall
i\in\\{1,2,3,4\\},n\in{\cal N},m\in\\{1,..,b\\}$

with $\xi$ small and non-zero. Differentiability of the functions $g_{i}$ will
be required later in Section IV-C where we will convexify problem (P4)
introduced below. The expected values in (12) can be estimated empirically via
sample averaging for a sufficiently large number of samples $N_{s}$, leading
to a new formulation of the optimization problem:

$\displaystyle(\text{P4})\min_{{\bf x},u_{i,n}^{m}}$
$\displaystyle\quad\sum_{m=1}^{b}\quad\frac{1}{N_{s}}\sum_{s=1}^{N_{s}}h_{0}({\bf
x};{\mbox{\boldmath$\rho$}}^{m}[s])$ (14a) s.t.
$\displaystyle\quad\frac{1}{N_{s}}\sum_{s=1}^{N_{s}}g_{i,n}({\bf
x},u_{i,n}^{m};{\mbox{\boldmath$\rho$}}^{m}[s]))\leq 0$
$\displaystyle\qquad\forall i\in\\{1,2,3,4\\},n\in{\cal N},\,m\in\\{1,..,b\\}$
(14b) $\displaystyle\quad h_{i}({\bf x})\leq 0\quad\forall i\in\\{5,6,7,8\\}$
(14c)

where we will draw $N_{s}$ samples
${{\mbox{\boldmath$\rho$}}^{m}[s]}_{s=1}^{N_{s}}$ of the random vector
${\mbox{\boldmath$\rho$}}^{m}$. Problem (P4) constitutes a conservative
approximation of the initial chance constrained problem (P3), meaning that an
optimal solution to (P4) is a feasible suboptimal solution to (P3).

###  IV-C Solution via successive convex approximation

At first one might try to solve problem (P4) with any software package for
nonlinear optimization. However, it is not straightforward to implement the
inverse matrix contained in ${\bf q}({\bf x})$ in a computationally efficient
way. Therefore, we seek a different strategy that may be computationally more
affordable. In particular, we will leverage the algorithm proposed in [22]
which follows the ideas of SCA methods. More specifically, the method solves a
sequence of strongly convex inner approximation of an initial non-convex
problem. In particular, each intermediate problem is strongly convex and can
be written as:

$\displaystyle(\text{P5}({\bf x}_{p}))\min_{{\bf x},u_{i,n}^{m}}$
$\displaystyle\quad\sum_{m=1}^{b}\quad\frac{1}{N_{s}}\sum_{s=1}^{N_{s}}\widetilde{h}_{0}({\bf
x};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf x}_{p})$ (15) s.t.
$\displaystyle\quad\frac{1}{N_{s}}\sum_{s=1}^{N_{s}}\widetilde{g}_{i,n}({\bf
x},u_{i,n}^{m};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf x}_{p})\leq 0$
$\displaystyle\qquad\forall i\in\\{1,2,3,4\\},n\in{\cal N},\,m\in\\{1,..,b\\}$
$\displaystyle\quad h_{i}({\bf x})\leq 0\quad\forall i\in\\{5,6,7,8\\}$

where $\widetilde{g}_{i,n}({\bf
x},u_{i,n}^{m};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf x}_{p})$ approximates
$g_{i,n}({\bf x},u_{i,n}^{m};{\mbox{\boldmath$\rho$}}^{m}[s])$ around ${\bf
x}={\bf x}_{p}$. For given samples ${\mbox{\boldmath$\rho$}}^{m}[s]$, the
problem $(\text{P5}({\bf x}_{p}))$ is solved for successive values of ${\bf
x}_{p}$ until convergence. The surrogate functions in (15) are defined as:

$\widetilde{h}_{0}({\bf x};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf x}_{p})=\|{\bf
q}({\bf x}_{p})+({\bf x}-{\bf x}_{p})^{\top}\nabla{\bf q}({\bf
x}_{p})\|^{2}\\\ +\frac{d}{2}\|{\bf x}-{\bf x}_{p}\|^{2}$ (16)

and

$\displaystyle\widetilde{g}_{i,n}({\bf
x},u_{i,n}^{m};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf x}_{p})=$ (17)
$\displaystyle\quad\frac{1}{2}\left(\widetilde{h}_{i,n}({\bf
x})+u_{i,n}^{m}+\sqrt{\xi^{2}+(\widetilde{h}_{i,n}({\bf
x})+u_{i,n}^{m})^{2}}\right)-u_{i,n}^{m}\epsilon$ $\displaystyle\qquad\forall
i\in\\{1,2,3,4\\},n\in{\cal N},m\in\\{1,..,b\\}$

with

$\displaystyle\widetilde{h}_{i,n}({\bf x};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf
x}_{p})=$ (18) $\displaystyle h_{i,n}({\bf x}_{p})+({\bf x}-{\bf
x}_{p})^{\top}\nabla h_{i,n}({\bf x}_{p})+({\bf x}-{\bf
x}_{p})^{\top}M_{i,n}({\bf x}-{\bf x}_{p})$ $\displaystyle\quad\forall
i\in\\{1,2,3,4\\},n\in{\cal N},m\in\\{1,..,b\\}$

where $M_{i,n}\in\mathbb{R}^{2\times 2}$ is derived according to Appendix -A
to ensure that $\widetilde{h}_{i,n}({\bf
x};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf x}_{p})$ is a global majorizer of
$h_{i,n}({\bf x}_{p};{\mbox{\boldmath$\rho$}})$. To lighten the notation, we
omit the ${\mbox{\boldmath$\rho$}}^{m}$ dependence on the RHS of equations
(16–18) but recall that ${\bf q}({\bf
x};{\mbox{\boldmath$\rho$}})=\left[X+\mathbf{1}^{\top}{\bf
x}\mathds{I}\right]^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}})$ with
${\mbox{\boldmath$\rho$}}=H(\mathbf{0},{\bf z})$, and the functions
$h_{i,n}({\bf x};{\mbox{\boldmath$\rho$}})$ are defined in equation (9). The
approximated functions $\widetilde{h}_{0}$, $\widetilde{g}_{i}$ and
$\widetilde{h}_{i}$ defined in equations (16–18) must satisfy some assumptions
described in [22], specified in Appendix -B and verified with our choice of
surrogate functions, in order to be suitable for the SCA method.

Next we present our algorithm to solve (P5(${\bf x}_{p}$)). First, let us
define the set ${\cal K}$ defined by equations (14c), i.e., the set defined by
the convex constraints of problem (P4). Let us also define the set ${\cal X}$
defined by equations (14c) and (14b), such that ${\cal X}\subset{\cal K}$.
Then, Algorithm 1 is proposed to solve (P5(${\bf x}_{p}$)).

Algorithm 1 _Optimal Gain Design via SCA_

Initialization: $\gamma_{p}\in(0,1]$, ${\bf x}_{0}\in{\cal X}$. Set $p=0$.
[1.] If $\|{\bf x}_{p}-{\bf x}_{p-1}\|<e$ with $e>0$, then STOP.
[2.] Compute the solution ${\bf x}^{*}({\bf x}_{p})$ of (P5(${\bf x}_{p}$)).
[3.] Set ${\bf x}_{p+1}={\bf x}_{p}+\gamma_{p}({\bf x}^{*}({\bf x}_{p})-{\bf
x}_{p})$
[4.] $p\leftarrow p+1$ and go to step 1.

The convergence properties of Algorithm 1 are based on [22, Theorem 1], and
they are presented in Appendix -B.

##  V Implementation of the controllers

We assume that each controller is equipped with sensing capabilities, i.e., it
is capable of measuring the voltage magnitudes at the node where it is
located. For any given DER $g\in\mathcal{G}$ connected to node $n\in{\cal N}$,
the following incremental Volt/Var control is implemented:

$\displaystyle q_{g,k+1}$
$\displaystyle=q_{g,k}+\eta(1-v_{n,k})-(1-\eta)\alpha q_{g,k}$ (19a)
$\displaystyle p_{g,k+1}$
$\displaystyle=\min{\left(p_{g,k},\sqrt{s_{g}^{2}-q_{g,k+1}^{2}}\right)},$
(19b)

with $s_{g}$ the nominal rated size of DER $g$. Equation (19a) represents the
reactive power update of DER $g$ connected at node $n$. The approximated
voltage $\nu_{n,k}$ written in the initial controller formulation (3) has been
replaced by the voltage measurement $v_{n,k}$, which makes this controller
fully decentralized (we no longer rely on the matrix $X$, and the impact of
other DERs’ reactive power is implicitly taken into account through the
network feedback). Equation (19b) indicates that we prioritize reactive power
over active power. By prioritizing reactive power over active power, we
further mitigate overvoltage issues as the active power injection is reduced
and reactive compensation is used. However, this induces active power
curtailment which is costly. We will investigate this issue in the next
section, where we will present our numerical results. The good behavior of the
controller with reactive power prioritization is verified throughout
simulations. We leave for future work the theoretical stability study of our
controller when considering active power curtailment.

Figure 3: Block diagram of the proposed framework.

Fig. 3 illustrates the different stages of the framework proposed in this
work. We compute the gains $\eta$ and $\alpha$ for a given time interval
$T\tau$ by solving the problem (P5$({\bf x}_{p})$) until convergence based on
forecasts of ${\bf z}_{k}$. The samples ${\mbox{\boldmath$\rho$}}^{m}[s]$ are
generated by solving multiple power flows for different values of ${\bf
z}_{k}=\bar{{\bf z}}_{k}+{\mbox{\boldmath$\delta$}}_{k}$, where
${\mbox{\boldmath$\delta$}}_{k}$ follows a given probability distribution
function. The gains can be computed the day before deployment, or hours ahead
depending on the forecast availability and the computational time required to
solve (P5$({\bf x}_{p})$) until convergence. They are then broadcast to the
controllable DERs. Notice that we do not need to differentiate between DERs as
the gains are the same for any DER connected to the network.

We are now ready to address the assumptions introduced in Section II-B. First,
our current framework enforces one and only one DER per node. This requirement
is rather restrictive, even though it has already been adopted in the
literature, e.g. [27]. We can easily relax this assumption by considering only
the entries of matrix $X$ that correspond to the nodes where a controllable
DER is located. The drawback is that we can only guarantee voltage
satisfaction for a subset ${\cal N}_{\mathrm{red}}\subset{\cal N}$ of nodes.
However, the effect of other controllable devices can be embedded in the map
$H$. For instance, on-load tap changers (OLTC) or switched capacitor banks can
drive the forecast voltage profiles
${\mbox{\boldmath$\rho$}}=H(\mathbf{0},{\bf z})$ inside the feasible set
${\cal V}$. Our methodology can be combined with other traditional regulation
methods, and an optimal combination of slow acting controllers, such as OLTC,
with our fast acting controllers is part of our future work. Furthermore, if
multiple DERs are connected to the same node, one can aggregate the DERs and
model them as one single device associated with one controller. The reactive
power setpoint produced by controller (3) is then appropriately dispatched to
the different DERs.

Second, in the present setup the network topology does not change with time.
However, in our current framework, topological changes that can be forecast
(because of maintenance or planned operation) can easily be integrated.
Indeed, those changes impact the matrix $X$. By appropriately choosing the
time interval $T\tau$ and recomputing $X$, one can derive gains that would be
well adapted to this new network topology. This assumption is much harder to
relax for learning-based methods such as the one proposed in [19] since it
requires building new datasets, and learning new equilibrium functions, which
can be time-consuming. When it comes to unplanned changes, e.g. sudden line
tripping or unplanned operations, our controllers take into account the
network conditions as a feedback, and do not worsen the situation.
Nevertheless, evaluation of the robustness with respect to unplanned changes
remains to be investigated in our future work.

Finally, Assumption 2 tackles the feasibility issue of the optimal reactive
power flow problem. It may happen that, under our controller architecture,
there is not enough reactive power reserve to satisfy the voltage constraints.
This problem is implicitly addressed through our chance-constrained
formulation. Indeed, increasing the value of $\epsilon$ enlarges the feasible
set. For $\epsilon=1$, the problem (P5$({\bf x}_{p})$) is always feasible.

##  VI Numerical experiments

(a)

(b)

(c)

Figure 4: (a) Low-voltage 42-nodes network, (b) aggregated non-controllable
active power injections, and active/reactive power consumption, (c) reactive
power setpoint update for controller at node 20 around hour 20:00.

We consider the low voltage network (0.4 kV) shown in Fig. 4(a). We used a
modified network from [31], in which photovoltaic power (PV) plants have been
placed at each node, with inverter-rated size picked randomly among
$\\{20,25,31\\}$ kVA. The DERs dynamics are not implemented, as they are
considered to be much faster than the controller dynamics. This is a
reasonable assumption because of the time-scale separation between the power
system phenomena and the different control loops as mentioned in [16]. As
such, when the controller update law produces a new reactive power setpoint,
it is instantaneously implemented by the controllable DER. In this paper, we
only consider PV plants as DERs, but any type of inverter-interfaced
generation for which the reactive power can be controlled could be considered.
Fig. 4(b) shows the aggregated loads and maximum available active power for PV
plants throughout the day. The data is from the Open Power System Data222Data
available at https://data.open-power-system-
data.org/household_data/2020-04-15, and have been modified to match the
initial loads and PV plants nominal values present in the network. The
reactive power demand is set such that the power factor is 0.95 (lagging).
This represents a typical summer day, with high PV production. We will show
that, under these conditions, the electrical DN undergoes both overvoltages
and undervoltages.

###  VI-A Simulation setup

In the following, we assume the controllable DERs to be equipped with the
following overvoltage protection.

_Overvoltage protection scheme:_ We consider an overvoltage protection for PV
plants, i.e., the plant is disconnected from the grid if the voltage goes
above 1.06 pu, or stays above 1.05 pu for 10 minutes. The DER reconnects if
the voltage remains below 1.05 pu for at least 1 minute. The disconnection
scheme is inspired by the CENELEC EN50549-2 standard [32], and has been
adjusted considering the voltage service limits used in this paper.

In the simulation, the voltage service limits are set to 1.05 and 0.95 p.u.,
respectively. The load and PV production profiles have a granularity of 1
second, i.e., active/reactive power consumption and maximum available active
power for PV plants change every 1 second. The time horizon $T\tau$ is set to
1 hour, and the forecast update $\Delta\tau$ to 30 minutes. The reactive power
setpoints update $\tau$ is set to 100 ms. We compare our proposed controller,
based on Algorithm 1 (OGD-CCO), with: (a) a static Volt/Var control (VoltVar);
and, (b) no control (ON/OFF).

#### Static Volt/Var control

It is inspired by the standard _IEEE Std 1547-2018_ , with maximum reactive
power consumed/absorbed set to 44% of the nominal power of the DER and reached
for voltages 1.05/0.95 pu, respectively. The deadband ranges between 0.99 and
1.01 pu. Active power prioritization is implemented. Therefore, the maximum
reactive power that can be produced or consumed corresponds to the minimum
between 44% of the nominal apparent power $s_{g}$ and the reactive power
reserve of the inverter $\sqrt{s^{2}_{g}-p^{2}_{g}}$, with $p_{g}$ the active
power injection of DER $g$.

#### No control

Set controller gains $\eta=\alpha=0$, and reactive powers to 0.

###  VI-B Results

In Fig. 5, we illustrate the dynamics of the controllers associated with DER
connected to two different nodes and the corresponding voltage magnitudes.

Figure 5: Controller dynamics and voltage magnitudes for two given nodes: node
34 and node 20.

We compare the voltage magnitudes resulting from the deployment of different
controllers against the case without controllers and without overvoltage
protection. We fixed $T\tau$ to 1 hour, i.e., the controller gains change
every hour. It is clear that the VoltVar control overuses reactive power, as
there is no need to control the voltages between hours 16:00 – 17:00, or after
22:00. The ON/OFF strategy sees large fluctuations in voltages due to constant
connection and disconnection of DERs. Fig. 4(c) illustrates the reactive power
setpoint update for the DER connected at node 20 at hour 20:00. In about 7
seconds, the reactive power setpoint varies with a magnitude of 5kvar, which
we believe is reasonable. However, if the reactive power setpoints are ramp
limited, this would not drastically affect the results since such large
variations can occur only during change of hours, and for certain specific
changes (in this scenario, at hours 8:00, 16:00, 18:00, 20:00, 21:00 and
22:00).

Figure 6: Cumulative distribution function for the vectors ${\bf
V}_{\mathrm{max}}=\\{\max{v}_{i,k}\\}_{i\in{\cal N},k\in\mathcal{T}_{h1-h2}}$,
and ${\bf V}_{\mathrm{min}}=\\{\min{v}_{i,k}\\}_{i\in{\cal
N},k\in\mathcal{T}_{h1-h2}}$, where $\mathcal{T}_{h1-h2}$ is the set of time
indices between hours $h_{1}$ and $h_{2}$.

The cumulative distribution functions (CDF) for maximum and minimum voltages
during the simulation are shown in Fig. 6. At each time $t$ between hours
$h_{1}$ and $h_{2}$, we pick the maximum and minimum voltages throughout the
network and store them into two distinct vectors. We then plot the cumulative
distribution functions of those vectors. For ${\cal T}_{0-24}$, i.e. for the
entire duration of the simulation, our proposed method crosses the 1.05 p.u.
line above 95% for the maximum voltage, and lower than 5% for the
undervoltage. This means that the voltage constraints are satisfied at least
95% of the time, which is consistent with our choice of $\epsilon=0.05$ and
$\epsilon=0.2$. We also show the CDF during two other intervals, between hours
7:00 – 8:00 for overvoltages and 18:00 – 19:00 for undervoltages. Notice that
a larger $\epsilon$ allows for more voltage constraint violations. The VoltVar
control sees overvoltages more than 10% of the time while performing similarly
to our method for undervoltages. The ON/OFF strategy sees overvoltages more
than 20% of the time.

Fig. 7 displays, from top to bottom, the energy lost in the lines, the energy
lost because of active power curtailment, and the cumulative reactive energy
usage for the different strategies considered in this work. The line losses
are computed based on the line currents squared times the line resistances.
Although the usage of reactive power is _practically free_ , it induces larger
line currents, hence larger power losses. We make a distinction between losses
caused by active power curtailment and line losses since the former are
covered by the network user while the latter are covered by the system
operator. However, the system operator can in turn increase network tariffs to
compensate for line losses due to overuse of reactive power compensation.

Our controller applies reactive power prioritization, which naturally induces
active power curtailment if the DER injects a large amount of active power
into the network. On the other hand, the VoltVar and ON/OFF strategies may
experience active power curtailment because of prolonged overvoltage
violations. Taking all these competing factor into consideration, the total
energy loss is much more important for the ON/OFF strategy, while it is
equivalent for our controller and the VoltVar control. One can see in Fig. 7
that the VoltVar control uses reactive power all day long, even when it is not
needed (e.g. between hours 20:00 – 21:00), while our proposed method better
manages the resources. One can also see that the parameter $\epsilon$ acts as
a lever on the total usage of reactive power. If the reactive power is very
expensive, and voltage violations are not extremely important, one could
increase the value of $\epsilon$ to reduce the total usage of reactive energy.

Figure 7: Energy lost in the lines and because of curtailment, and reactive
energy usage.

In Fig. 8, we show the maximum duration of voltage violations. The vector
$\Gamma_{{\cal V}}$ is built such that, if an overvoltage or undervoltage
occurs in the network at a given time step $k\tau$, we count and add up every
voltage violations for the subsequent time steps. When, at a subsequent time
step, no voltage violation occurs, the total number of voltage violation is
appended to the vector, and the counter is reset to 0. The vector
$\Gamma_{{\cal V}}$ indicates the duration of voltage violations. This is an
important number, since usually electrical devices can cope with short and
limited over or under voltages, but may be damaged during prolonged, large
excursions from nominal voltages. The parameter $\epsilon$ sets the maximum
time the voltage can exceed the voltage limits. For our given $\Delta\tau$ =
30 minutes, and with $\epsilon=5\%$, the maximum duration for voltage
violation is 90 seconds, while for $\epsilon=20\%$, it is 360 seconds.

Fig. 8 shows that a smaller $\epsilon$ leads to shorter voltage violations,
while both ON/OFF and VoltVar strategies lead to substantially longer voltage
violations.

Figure 8: Cumulative distribution function for voltage violations.

##  VII Conclusion

We propose an incremental Volt/Var control strategy for voltage regulation in
DNs. We show the stability of our controller, and introduce a methodology to
compute the gains of our controller based on a chance-constrained formulation
of an optimal reactive power flow problem. Our methodology only needs limited
offline communications, i.e., the same controller gains are broadcast to
individual controllers. Our chance-constrained formulation tackles
uncertainties in power injections. Moreover, the feasibility issue of local
Volt/Var control, i.e., if there is enough reactive power reserve to satisfy
the voltage constraints given the architecture of the controller, is
implicitly taken into account by allowing the system operator to tolerate a
prescribed probability of voltage violations. Our method shows better
performance compared to static Volt/Var curves with fixed parameters, and
limited and short voltage violations that are consistent with the prescribed
probability. Future works will investigate the combination of our fast-acting
controller with slower traditional regulation devices.

###  -A Derivations of $M$

In this section, we derive the matrices $M_{i,n}$ appearing in (18) to ensure
that $\widetilde{h}_{i,n}({\bf x};{\mbox{\boldmath$\rho$}}^{m}[s],{\bf
x}_{p})$ is a global majorizer of $h_{i,n}({\bf
x}_{p};{\mbox{\boldmath$\rho$}})$. First, let us build the Jacobian matrices
$\nabla h_{i}({\bf x})$ for $i=1,2,3,4$ from the Jacobian matrix $\nabla{\bf
q}({\bf x})$. Notice that the partial derivatives of the vector ${\bf q}({\bf
x})$ with respect to the first and second component of the vector ${\bf x}$
are identical:

$\displaystyle\frac{\partial{\bf q}({\bf x})}{\partial
x_{1}}=\frac{\partial{\bf q}({\bf x})}{\partial x_{2}}$
$\displaystyle=P\mathsf{diag}{(\gamma({\bf
x}))}P^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}})$ (20)

where $\gamma({\bf x})=[-\frac{1}{(\mathbf{1}^{\top}{\bf
x}+\lambda_{X,1})^{2}},...,-\frac{1}{(\mathbf{1}^{\top}{\bf
x}+\lambda_{X,N})^{2}}]^{\top}$. Since $X$ is positive definite, and
$\mathbf{1}^{\top}{\bf x}\geq 0$ by constraint, the vector $\gamma({\bf x})$
has bounded elements. We can now construct $\nabla{\bf q}({\bf
x})=\left[\frac{\partial{\bf q}({\bf x})}{\partial x_{1}},\frac{\partial{\bf
q}({\bf x})}{\partial x_{2}}\right]^{\top}\in\mathbb{R}^{2\times N}$ as

$\nabla{\bf q}({\bf x})=[P\mathsf{diag}{(\gamma({\bf
x}))}P^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}}),P\mathsf{diag}{(\gamma({\bf
x}))}P^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}})]^{\top}.$ (21)

Recalling the definitions (9) of the vector-valued functions $h_{i}$, we can
write:

$\displaystyle\nabla h_{1}({\bf x})=\nabla{\bf q}({\bf x}),\quad\nabla
h_{2}({\bf x})=-\nabla{\bf q}({\bf x})$ (22) $\displaystyle\nabla h_{3}({\bf
x})=X\nabla{\bf q}({\bf x}),\quad\nabla h_{4}({\bf x})=-X\nabla{\bf q}({\bf
x})$

In order to determine $M_{i,n}$ such that $M_{i,n}-\nabla^{2}h_{i,n}({\bf
x})\succeq 0$, we first derive the Hessian matrix of $h_{i,n}({\bf x})$:

$\nabla^{2}h_{i,n}({\bf x})=c_{i,n}({\bf x})\begin{pmatrix}1&1\\\
1&1\end{pmatrix}$ (23)

where $c_{i,n}({\bf x})\in\mathbb{R}$. We define $c_{i}({\bf x})=[c_{i,1}({\bf
x}),...,c_{i,N}({\bf x})]^{\top}$ and from eqs.(21) and (22) it is easy to
retrieve:

$\displaystyle c_{1}({\bf x})=P\mathsf{diag}{(d_{\gamma}({\bf
x}))}P^{-1}(\mathds{1}-{\mbox{\boldmath$\rho$}}),\quad c_{2}({\bf
x})=-c_{1}({\bf x})$ (24) $\displaystyle c_{3}({\bf x})=Xc_{1}({\bf x}),\quad
c_{4}({\bf x})=-Xc_{1}({\bf x})$

where $d_{\gamma}({\bf x})=[\frac{2}{(\mathbf{1}^{\top}{\bf
x}+\lambda_{X,1})^{3}},...,\frac{2}{(\mathbf{1}^{\top}{\bf
x}+\lambda_{X,N})^{3}}]^{\top}$. Given that the quantity
$\mathbf{1}^{\top}{\bf x}$ is always positive, setting $\mathbf{1}^{\top}{\bf
x}=0$ in $d_{\gamma}({\bf x})$ gives us an upper bound on the value of
$\nabla^{2}h_{i,n}({\bf x})$. Therefore, we write $\bar{c}_{i}=c_{i}({\bf
x})|_{\mathbf{1}^{\top}{\bf x}=0}$ for $i=1,2,3,4$ and find:

$M_{i,n}=2\begin{pmatrix}\max(0,\bar{c}_{i,n})&0\\\
0&\max(0,\bar{c}_{i,n})\end{pmatrix}\quad\forall n\in{\cal N}$ (25)

such that $\widetilde{h}_{i,n}({\bf x};{\bf x}_{p})$ is a global majorizer of
$h_{i,n}({\bf x})$ $\forall{\bf x}\in\mathbb{R}^{2}$.

###  -B Convergence guarantees for Algorithm 1

In this section, we show that the assumptions considered in [22] for the SCA
method to converge are applicable to our problem (P4). To this end, define
${\cal W}\in\mathbb{R}^{N}$ such that ${\mbox{\boldmath$\rho$}}\in{\cal W}$.

We start by observing the following properties and facts related to problem
(P4).

[A1] The set ${\cal K}$ is closed and convex. The set ${\cal K}$ is
constructed with linear constraints.

[A2] The functions ${\bf x}\mapsto h_{0}({\bf x},{\mbox{\boldmath$\rho$}})$
and ${\bf x}\mapsto g_{i,n}({\bf x},u_{i,n};{\mbox{\boldmath$\rho$}})$ are
continuously differentiable on ${\cal K}$, for any
${\mbox{\boldmath$\rho$}}\in{\cal W}$. This follows from the differentiability
of ${\bf q}({\bf x};{\mbox{\boldmath$\rho$}})$; since $X$ is positive definite
by construction, and $\mathbf{1}^{\top}{\bf x}\geq 0$ by constraint, it
follows that $\nabla{\bf q}({\bf x})$ always exists and is bounded.

[A3] The gradient $\nabla_{{\bf x}}h_{0}({\bf x};{\mbox{\boldmath$\rho$}})$ is
Lipschitz continuous on ${\cal K}$ with constant $L_{\nabla h_{0}}$, uniformly
in $\rho$. This follows from the fact that $\nabla^{2}{\bf q}({\bf x})$ is
bounded on ${\cal K}$.

[A4] For some ${\bf x}_{0}\in{\cal X}$, $\\{{\bf x}\in{\cal X}:h_{0}({\bf
x};{\mbox{\boldmath$\rho$}})\leq h_{0}({\bf
x}_{0};{\mbox{\boldmath$\rho$}})\\}$ is compact, for any given
${\mbox{\boldmath$\rho$}}\in{\cal W}$.

Next, we note that the surrogate function $\widetilde{h}_{0}:{\cal
K}\times{\cal W}\times{\cal X}\rightarrow\mathbb{R}$ is continuously
differentiable with respect to ${\bf x}$ and satisfies the following
conditions.

[A5] The map ${\bf x}\mapsto\widetilde{h}_{0}({\bf
x};{\mbox{\boldmath$\rho$}},{\bf x}_{p})$ is strongly convex on ${\cal K}$
with constant $c\geq d>0$, for any given ${\bf x}_{p}\in{\cal X}$ and
${\mbox{\boldmath$\rho$}}\in{\cal W}$.

[A6] We have that $\nabla\widetilde{h}_{0}({\bf
x}_{p};{\mbox{\boldmath$\rho$}},{\bf x}_{p})=\nabla h_{0}({\bf
x}_{p};{\mbox{\boldmath$\rho$}})$ for any given ${\bf x}_{p}\in{\cal X}$ and
${\mbox{\boldmath$\rho$}}\in{\cal W}$. to see this, define
$f(\cdot)=\|\cdot\|^{2}$ so and express $h_{0}({\bf x})=f({\bf q}({\bf x}))$
and $\widetilde{h}_{0}({\bf x};{\bf x}_{p})=f({\bf q}({\bf x}_{p})+({\bf
x}-{\bf x}_{p})^{\top}\nabla{\bf q}({\bf x}_{p}))+\frac{d}{2}\|{\bf x}-{\bf
x}_{p}\|^{2}$. We can write $\nabla h_{0}({\bf x})=\nabla{\bf q}({\bf
x})\nabla f({\bf q}({\bf x}))$, and $\nabla\widetilde{h}_{0}({\bf x};{\bf
x}_{p})=\nabla{\bf q}({\bf x})\nabla f({\bf q}({\bf x})+({\bf x}-{\bf
x}_{p})^{\top}\nabla{\bf q}({\bf x}))+d({\bf x}-{\bf x}_{p})$. It is clear
that at ${\bf x}={\bf x}_{p}$ we have $\nabla\widetilde{h}_{0}({\bf
x}_{p};{\bf x}_{p})=\nabla h_{0}({\bf x}_{p})$.

[A7] The map $({\bf x},{\bf x}_{p})\mapsto\nabla\widetilde{h}_{0}({\bf
x};{\mbox{\boldmath$\rho$}},{\bf x}_{p})$ is continuous on ${\cal
K}\times{\cal X}$ for any given ${\mbox{\boldmath$\rho$}}\in{\cal W}$.

Finally, the surrogate constraint functions $\widetilde{g}_{i,n}:{\cal
K}\times\mathbb{R}_{\geq 0}\times{\cal W}\times{\cal X}\rightarrow\mathbb{R}\
\forall i\in\\{1,2,3,4\\},n\in{\cal N}$ satisfy the following.

[A8] The map $({\bf x},u_{i,n})\mapsto\widetilde{g}_{i,n}({\bf
x},u_{i,n};{\mbox{\boldmath$\rho$}},{\bf x}_{p})$ is convex on ${\cal
K}\times\mathbb{R}_{\geq 0}$ for any given ${\mbox{\boldmath$\rho$}}\in{\cal
W}$ and ${\bf x}_{p}\in{\cal X}$. to see this, note that the composite
function $\widetilde{g}_{i,n}=f(\widetilde{h}_{i,n}+u_{i,n})-u_{i,n}\epsilon$
is given by the outer function
$f(y)=\frac{1}{2}\left(y+\sqrt{\xi^{2}+(y)^{2}}\right)$ which is strictly
increasing and convex and by the inner function $\widetilde{h}_{i,n}+u_{i,n}$
which is convex by construction; it then follows that $\widetilde{g}_{i,n}$ is
also convex.

[A9] $\widetilde{g}_{i,n}({\bf x}_{p},u_{i,n};{\mbox{\boldmath$\rho$}},{\bf
x}_{p})=g_{i,n}({\bf x}_{p},u_{i,n};{\mbox{\boldmath$\rho$}})$ for any given
${\bf x}_{p}\in{\cal X}$, $u_{i,n}\in\mathbb{R}_{\geq 0}$ and
${\mbox{\boldmath$\rho$}}\in{\cal W}$.

[A10] $g_{i,n}({\bf
x},u_{i,n};{\mbox{\boldmath$\rho$}})\leq\widetilde{g}_{i,n}({\bf
x},u_{i,n};{\mbox{\boldmath$\rho$}},{\bf x}_{p})$ for all ${\bf x}\in{\cal
K}$, $u_{i,n}\in\mathbb{R}_{\geq 0}$ and for any given ${\bf x}_{p}\in{\cal
X}$, ${\mbox{\boldmath$\rho$}}\in{\cal W}$.

[A11] The map $({\bf x},u_{i,n},{\bf x}_{p})\mapsto\widetilde{g}_{i,n}({\bf
x},u_{i,n};{\mbox{\boldmath$\rho$}},{\bf x}_{p})$ is continuous on ${\cal
K}\times\mathbb{R}_{\geq 0}\times{\cal X}$ for any given
${\mbox{\boldmath$\rho$}}\in{\cal W}$.

[A12] $\nabla\widetilde{g}_{i,n}({\bf
x}_{p},u_{i,n,p}^{m};{\mbox{\boldmath$\rho$}},{\bf x}_{p})=\nabla g_{i,n}({\bf
x}_{p},u_{i,n,p}^{m};{\mbox{\boldmath$\rho$}})$ for all ${\bf x}_{p}\in{\cal
X}$, $u_{i,n,p}^{m}\in\mathbb{R}_{\geq 0}$ and for any given ${\bf
x}_{p}\in{\cal X}$, ${\mbox{\boldmath$\rho$}}\in{\cal W}$. To see this, using
the chain rule we obtain

$\displaystyle\frac{\partial g_{i,n}({\bf x},u)}{\partial{\bf x}}$
$\displaystyle=\frac{1}{2}\left(1+\frac{h_{i,n}({\bf
x})+u}{\sqrt{\xi^{2}+(h_{i,n}({\bf x})+u)^{2}}}\right)\frac{\partial
h_{i,n}}{\partial{\bf x}},$ (26) $\displaystyle\frac{\partial g_{i,n}({\bf
x},u)}{\partial u}$ $\displaystyle=\frac{1}{2}\left(1+\frac{h_{i,n}({\bf
x})+u}{\sqrt{\xi^{2}+(h_{i,n}({\bf x})+u)^{2}}}\right)-\epsilon.$

The claim is verified since $\frac{\partial\widetilde{h}_{i,n}}{\partial{\bf
x}}({\bf x}_{p};{\mbox{\boldmath$\rho$}},{\bf x}_{p})=\frac{\partial
h_{i,n}}{\partial{\bf x}}({\bf x}_{p};{\mbox{\boldmath$\rho$}})$ and
$\widetilde{h}_{i,n}({\bf x}_{p};{\mbox{\boldmath$\rho$}},{\bf
x}_{p})=h_{i,n}({\bf x}_{p};{\mbox{\boldmath$\rho$}})$ for all ${\bf
x}_{p}\in{\cal X}$ and any given ${\mbox{\boldmath$\rho$}}\in{\cal W}$.

[A13] The map $({\bf x},u_{i,n},{\bf
x}_{p})\mapsto\nabla\widetilde{g}_{i,n}({\bf
x},u_{i,n};{\mbox{\boldmath$\rho$}},{\bf x}_{p})$ is continuous on ${\cal
K}\times\mathbb{R}_{\geq 0}\times{\cal X}$ for any given
${\mbox{\boldmath$\rho$}}\in{\cal W}$.

Given the above properties [A1]–[A13], the following holds when using
Algorithm 1 to solve (P4).

###### Proposition 1

Consider the problem (P4), and use Algorithm 1 with a step-size $\gamma_{p}$.
Assume that one of the following two conditions is satisfied:

  * •

$0<\inf_{p}\gamma_{p}\leq\sup_{p}\gamma_{p}\leq\gamma^{\mathrm{max}}\leq 1$
and $2c\geq\gamma^{\mathrm{max}}L_{\nabla h_{0}}$,

  * •

${\cal X}$ is compact, ${\bf x}^{*}({\bf x}_{p})\in{\cal X}$ is regular for
every ${\bf x}_{p}\in{\cal X}$ and $\gamma^{p}\in(0,1],\ \gamma_{p}\rightarrow
0,\ \sum_{p}\gamma_{p}=+\infty$.

Then, every regular limit point of $\\{{\bf x}_{p}\\}$ is a stationary
solution of (P4). Furthermore, none of such points is a local maximum of
$h_{0}({\bf x})$. $\Box$

The result of Proposition 1 follows from [22, Theorem 1]. It asserts that the
Algorithm 1 converges to a stationary solution of (P4), and this solution is
guaranteed not to be a local maximum.

## References

  * [1] R. Walling, R. Saint, R. C. Dugan, J. Burke, and L. A. Kojovic, “Summary of distributed resources impact on power delivery systems,” _IEEE Trans. on Power Delivery_ , vol. 23, no. 3, pp. 1636–1644, 2008.
  * [2] M. S. Alam, F. S. Al-Ismail, A. Salem, and M. A. Abido, “High-level penetration of renewable energy sources into grid utility: Challenges and solutions,” _IEEE Access_ , vol. 8, pp. 190 277–190 299, 2020.
  * [3] J. McDonald, “Solar power impacts power electronics in the smart grid,” _Power Electronics_ , 2013.
  * [4] K. E. Antoniadou-Plytaria, I. N. Kouveliotis-Lysikatos, P. S. Georgilakis, and N. D. Hatziargyriou, “Distributed and decentralized voltage control of smart distribution networks: Models, methods, and future research,” _IEEE Trans. on Smart Grid_ , vol. 8, no. 6, pp. 2999–3008, 2017.
  * [5] P. Srivastava, R. Haider, V. J. Nair, V. Venkataramanan, A. M. Annaswamy, and A. K. Srivastava, “Voltage regulation in distribution grids: A survey,” _Annual Reviews in Control_ , vol. 55, pp. 165–181, 2023.
  * [6] A. Colot, Y. Chen, B. Cornélusse, J. Cortés, and E. Dall’Anese, “Optimal power flow pursuit via feedback-based safe gradient flow,” _arXiv:2312.12267_ , 2023.
  * [7] L. Ortmann, J. Maeght, P. Panciatici, F. Dörfler, and S. Bolognani, “Online feedback optimization for subtransmission grid control,” _arXiv:2212.07795_ , 2023.
  * [8] H. Zhu and H. J. Liu, “Fast local voltage control under limited reactive power: Optimality and stability analysis,” _IEEE Trans. on Power Syst._ , vol. 31, no. 5, pp. 3794–3803, 2016.
  * [9] X. Zhou, M. Farivar, Z. Liu, L. Chen, and S. H. Low, “Reverse and forward engineering of local voltage control in distribution networks,” _IEEE Trans. on Automatic Control_ , vol. 66, no. 3, pp. 1116–1128, 2021.
  * [10] K. Baker, A. Bernstein, E. Dall’Anese, and C. Zhao, “Network-cognizant voltage droop control for distribution grids,” _IEEE Trans. on Power Syst._ , vol. 33, no. 2, pp. 2098–2108, 2018.
  * [11] P. Jahangiri and D. C. Aliprantis, “Distributed Volt/VAr control by PV inverters,” _IEEE Trans. on Power Syst._ , vol. 28, no. 3, pp. 3429–3439, 2013.
  * [12] M. Farivar, X. Zho, and L. Chen, “Local voltage control in distribution systems: An incremental control algorithm,” in _2015 IEEE International Conference on Smart Grid Communications (SmartGridComm)_ , 2015, pp. 732–737.
  * [13] J. Wei, S. Gupta, D. C. Aliprantis, and V. Kekatos, “A chance-constrained optimal design of Volt/VAR control rules for distributed energy resources,” in _2023 North American Power Symposium (NAPS)_ , 2023, pp. 1–6.
  * [14] F. U. Nazir, B. C. Pal, and R. A. Jabr, “A two-stage chance constrained volt/var control scheme for active distribution networks with nodal power uncertainties,” _IEEE Trans. on Power Syst._ , vol. 34, no. 1, pp. 314–325, 2019.
  * [15] S. Karagiannopoulos, P. Aristidou, and G. Hug, “Data-driven local control design for active distribution grids using off-line optimal power flow and machine learning techniques,” _IEEE Trans. on Smart Grid_ , vol. 10, no. 6, pp. 6461–6471, 2019.
  * [16] A. Eggli, S. Karagiannopoulos, S. Bolognani, and G. Hug, “Stability analysis and design of local control schemes in active distribution grids,” _IEEE Trans. on Power Syst._ , vol. 36, no. 3, pp. 1900–1909, 2020.
  * [17] S. Gupta, A. Mehrizi-Sani, S. Chatzivasileiadis, and V. Kekatos, “Deep learning for scalable optimal design of incremental Volt/VAR control rules,” _IEEE Control Systems Letters_ , vol. 7, pp. 1957–1962, 2023.
  * [18] S. Gupta, V. Kekatos, and S. Chatzivasileiadis, “Optimal design of Volt/VAR control rules of inverters using deep learning,” _arXiv:2211.09557_ , 2024.
  * [19] Z. Yuan, G. Cavraro, M. K. Singh, and J. Cortés, “Learning provably stable local Volt/Var controllers for efficient network operation,” _IEEE Trans. on Power Syst._ , vol. 39, no. 1, pp. 2066–2079, 2024.
  * [20] P. Serna Torre and P. Hidalgo-Gonzalez, “Decentralized optimal power flow for time-varying network topologies using machine learning,” _Electric Power Systems Research_ , vol. 212, p. 108575, 2022.
  * [21] A. Nemirovski and A. Shapiro, “Convex approximations of chance constrained programs,” _SIAM Journal on Optimization_ , vol. 17, no. 4, pp. 969–996, 2007.
  * [22] G. Scutari, F. Facchinei, L. Lampariello, and P. Song, “Parallel and distributed methods for nonconvex optimization,” in _2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2014, pp. 840–844.
  * [23] S. Bolognani and S. Zampieri, “On the existence and linear approximation of the power flow solution in power distribution networks,” _IEEE Trans. on Power Syst._ , vol. 31, no. 1, pp. 163–172, 2015.
  * [24] C. Wang, A. Bernstein, J.-Y. Le Boudec, and M. Paolone, “Explicit conditions on existence and uniqueness of load-flow solutions in distribution networks,” _IEEE Trans. on Smart Grid_ , vol. 9, no. 2, pp. 953–962, 2016.
  * [25] ——, “Existence and uniqueness of load-flow solutions in three-phase distribution networks,” _IEEE Trans. on Power Syst._ , vol. 32, no. 4, pp. 3319–3320, 2017.
  * [26] E. Dall’Anese and A. Simonetto, “Optimal power flow pursuit,” _IEEE Trans. on Smart Grid_ , vol. 9, no. 2, pp. 942–952, 2016.
  * [27] N. Li, G. Qu, and M. Dahleh, “Real-time decentralized voltage control in distribution networks,” in _2014 52nd Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , 2014, pp. 582–588.
  * [28] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Trans. on Power delivery_ , vol. 4, no. 2, pp. 1401–1407, 1989.
  * [29] M. Farivar, L. Chen, and S. Low, “Equilibrium and dynamics of local voltage control in distribution systems,” in _52nd IEEE Conference on Decision and Control_ , 2013, pp. 4329–4334.
  * [30] V. Kekatos, L. Zhang, G. B. Giannakis, and R. Baldick, “Voltage regulation algorithms for multiphase power distribution grids,” _IEEE Trans. on Power Syst._ , vol. 31, no. 5, pp. 3913–3923, 2015.
  * [31] D. Sarajlić and C. Rehtanz, “Low voltage benchmark distribution network models based on publicly available data,” in _2019 IEEE PES Innovative Smart Grid Technologies Europe (ISGT-Europe)_ , 2019, pp. 1–5.
  * [32] “Requirements for generating plants to be connected in parallel with distribution networks – Part 2: connection to a MV distribution network – generating plants up to and including type B,” _CENELEC_ , pp. 1–84, 2019.


[inst1]organization=Ansys Inc, CTO Office,addressline=2400 Ansys Dr.,
city=Canonsburg, postcode=15317, state=PA, country=USA

[inst2]organization=Ansys Inc, CTO Office,addressline=2645 Zanker Rd.,
city=San Jose, postcode=95134, state=CA, country=USA

[inst3]organization=Ansys Inc, ESOBU,addressline=2645 Zanker Rd., city=San
Jose, postcode=95134, state=CA, country=USA

# A Thermal Machine Learning Solver For Chip Simulation

Rishikesh Ranade  Haiyang He  Jay Pathak  Norman Chang  Akhilesh Kumar  Jimin
Wen

###### Abstract

Thermal analysis provides deeper insights into electronic chips’ behavior
under different temperature scenarios and enables faster design exploration.
However, obtaining detailed and accurate thermal profile on chip is very time-
consuming using FEM or CFD. Therefore, there is an urgent need for speeding up
the on-chip thermal solution to address various system scenarios. In this
paper, we propose a thermal machine-learning (ML) solver to speed-up thermal
simulations of chips. The thermal ML-Solver is an extension of the recent
novel approach, CoAEMLSim (Composable Autoencoder Machine Learning Simulator)
with modifications to the solution algorithm to handle constant and
distributed HTC. The proposed method is validated against commercial solvers,
such as Ansys MAPDL, as well as a latest ML baseline, UNet, under different
scenarios to demonstrate its enhanced accuracy, scalability, and
generalizability.

##  1 Introduction

It is well known that thermal issues can severely degrade the performance and
reliability of chips [17, 3, 10, 5, 9]. Overlarge peak temperatures and stiff
thermal gradients can fatally impact transistor performance, stress, aging,
electro-migration (EM), voltage drops and timing [19, 11]. Hence accurate
prediction of the maximum temperature and thermal gradient on the chip becomes
important for the performance and reliability of chip-packaging systems used
in several applications such as 5G, automobiles and computational hardware for
Artificial Intelligence. Conventional Finite Element Analysis (FEA) or
Computational Fluid Dynamics (CFD) based thermal analysis is computationally
expensive due to the enormous system parameter space in the form of stiff
powermaps and wide range of Heat Transfer Coefficients (HTCs), die thicknesses
and chip sizes. As a result, batches of simulations are required to be solved
from scratch every time new system parameters of electronic chips are
considered.

Recently, a multitude of machine learning methods have been proposed to
enhance and accelerate physics based numerical solvers in the context of
electronic chip simulations. For example, a Deep Neural Networks (DNN) based
fast static thermal solver has been proposed in Wen et al. [18] to generate a
high-resolution Delta T map. In a new branch of study, researchers have
employed Physics-Informed Neural Networks (PINNs) [12] to solve chip problems
[2]. Hennigh et al. [8] use PINNs to predict temperature profiles resulting
due to variations in the FGPA heat sink geometry for a uniform power map.
Similarly, Ranade et al. [14] and He and Pathak [7] use discretization-based
techniques in combination with neural networks to predict temperatures on
chips for powermaps sampled from a Gaussian distribution. Stipsitz and
Sanchis-Alepuz [16] presented a proof-of-concept approach to provide
approximate predictions of steady-state temperatures using convolutional
neural networks. Chhabria et al. [4] use a convolutional encoder-decoder
network to learn a mapping between powermap and temperature. Most studies in
relation to thermal simulation of chips make simplifying assumptions about the
power map or other system parameters such as HTC. The high dimensional
parameter space involved in on-chip thermal simulations usually makes it
challenging for conventional deep learning models to predict accurate
temperatures for unseen input parameters.

In this paper, we propose a thermal ML-Solver, which is based on a recent work
[13], for accurately simulating temperatures on electronic chips subjected to
a wide range of system parameters. Our approach extends the previous method to
modeling the complicated use case of electronic chips as well as proposes a
slight modification to the solution algorithm to account for global system
parameters as well as distributed HTCs. The input distributions considered in
this paper have noticeable discontinuities and wider ranges of magnitudes and
spatial distributions, which are extremely challenging to model. This paper
demonstrates how the proposed approach can accurately model the temperature
prediction widely ranging distributions of powermap and HTCs. The ability of
this approach to accurately model such high-dimensional power maps and HTCs
makes it unique and different from other ML approaches.

Significant contributions:

  * 1.

We propose a thermal ML-simulator based on the work by Ranade et al. [13] for
electronic chip simulations over a wide range of parameters with varying
spatial distributions.

  * 2.

We modify the existing solution algorithm to handle global system parameters.
This work can be considered as one of the first implementations of ML to
handle distributed HTCs.

  * 3.

Our approach entails solving a constrained equilibrium during evaluation,
similar to traditional PDE solvers. This results in better stability,
robustness and accuracy as opposed to blackbox ML methods.

###  1.1 Chip Thermal Analysis

Chip thermal analysis may involve multiple components of the chip and various
system parameters with high dimensionality and huge parameter space. This
study focuses on a simplified three-layer chip model consisting of a Si
substrate, an insulation layer and an interconnection layer. The schematic of
the three-layer model is illustrated in Fig. 1.

Figure 1: Schematic of the simplified three-layer chip model.

An arbitrary tile-based power map, which consists of rectangular regions of
different heat source, each representing a different functional block on the
chip, is sandwiched between the insulation layer and interconnection layer.
For instance, a $4000\mu m$ x $4000\mu m$ chip with $200\mu m$x$200\mu m$ tile
size will lead to a 20x20 array power map. An example of such power map is
shown in Fig. 2. Temperature distribution of the chip, especially the heating
plane where the power map is applied, is of significant interest. In practice,
such tile-based power could be random and incur drastic power gradient, which
makes it challenging for conventional data-driven based approaches to learn,
not to mention generalization to any unseen random power maps. Additionally,
the temperature distribution on the heating plane is greatly affected by
system parameters such as HTCs distribution on the boundary or die thickness,
the effect of which is difficult to capture. Hence, in this study we only
focus on accurately modeling chip temperature predictions on the heating plane
across a wide range of parameters. In a system-level simulation, the three-
layer die is cooled either by natural or forced convection. The cooling effect
of air is modeled here using HTC distributions on top and bottom surfaces.

Figure 2: Example of powermap and temperature on chip.

###  1.2 Data Generation

In this paper, we present two use cases with varying levels of complexity. In
the first case, the input parameters consist of a powermap distribution on the
heating plane, HTCs on boundaries and die thickness. Constant HTC values are
applied on the top and bottom boundaries of the die, whereas other boundaries
have adiabatic boundary conditions. The HTC varies in between $1e^{-6}$ and
$1e^{-5}$, while the die thickness varies between $20\mu m$ and $200\mu m$.
$5000$ numerical simulations are carried out utilizing a high-fidelity
numerical solver, Ansys MAPDL, for simulating the heat transfer on a $4000\mu
m$ chip with a grid resolution of $256$x$256$x$32$ for random power maps, heat
transfer coefficients and Die thicknesses. Examples of the power map and the
corresponding temperature on the chip are shown in Fig. 2.

The second use case presents a more challenging scenario close to industrial
applications, where the HTCs specified on the top and bottom boundaries also
have a spatial tile-based distribution similar to that of the power map and
the tile size is smaller. The HTC distribution on each tile are varied
randomly from $-1e^{-6}$ to $1e^{-6}$. Additionally, we reduce the tile size
of powermap and HTC from 200 $\mu m$ down to $50\mu m$. In this case, $2500$
numerical simulations are carried out using Ansys MAPDL for random power maps
and HTCs.

In both cases, other parameters such as interconnection layer thickness,
insulation layer thickness, Si substrate thermal conductivity, insulation
layer thermal conductivity and interconnection layer conductivity are kept
constant. The power on each tile is randomly sampled between $0$ and $300mW$
and the total power on the chip is conserved among all the training samples.
Even though the data is generated on a 3-D die, the solutions are retained
only on the chip where the power map is applied. Hence, the temperature is
predicted only on a surface but its spatial distribution and magnitude are
affected by system parameters such as HTC and die thickness. As a result,
these parameters are accounted for in the thermal ML-Solver along with power
map. The division between the training and testing datasets is $80/20$.

##  2 Thermal ML-Solver Details

A schematic of our model is described in Fig. 3. Here we describe the
components of our thermal ML-solver and the modifications to the existing
algorithm found in the work by Ranade et al. [13].

  1. 1.

Decomposition of computational domain: The computational domain corresponding
to the 2-D chip (or plane of interest in Fig. 1) is decomposed into subdomains
of equal physical sizes, where each subdomain consists of $m^{2}$
computational elements. As prescribed in Ranade et al. [13], based on an
extensive ablation study , the value of $m$ is set to $16$ resulting in around
$256$ subdomains for the $4000\mu m$ chip.

  2. 2.

Encoding on subdomains: An initial temperature of $T=300K$ and a given power
map, $P$ are encoded into lower-dimensional encodings $\eta$ and $s$ on every
subdomain using pre-trained encoders. The global system parameters,
$\gamma=HTC,Die_{z}$, are kept constant across all the subdomains in the first
use case. The conditioning of all the subdomains with global parameters
improves the predictive capability of the solver. In the second use case, HTCs
on top and bottom surfaces are decomposed in to subdomains and encoded similar
to the powermap but $Die_{z}$ is used as a global parameter. Even though the
HTC is applied on the top and bottom surface of the die, they are considered
to be applied on the plane of the chip for the sake of modeling a 2-D problem.

  3. 3.

Constrained Thermal solver: A constrained equilibrium problem is solved using
fixed point iterations where the encodings ($\eta,s,\gamma$) on each subdomain
and its neighbors are evaluated iteratively using a pretrained flux
conservation autoencoder. In each iteration, the output encoding corresponding
to temperature, $\eta^{\prime}$, evaluated by the flux conservation
autoencoder, is used to replace the input temperature encoding, $\eta$ for the
next iteration. Other encodings such as $s$ and $\gamma$ are kept fixed. The
fixed encodings serve as hard constraints and steer the temperature encoding
towards an equilibrium, which resembles the encoding of converged temperature.
The fixed point iteration is stopped when the following condition is met,
$L_{2}(\eta-\eta^{\prime})<1e^{-8}$.

  4. 4.

Decoding on subdomains: The converged temperature encodings on subdomains at
the end of the fixed point iterations are decoded into temperatures and post-
processed on the entire computational domain.

It may be observed that our approach strongly draws ideas and inspiration from
traditional PDE solvers. The novel iterative inferencing approach is designed
to provide stability and robustness to the solution methodology in comparison
to blackbox ML methods. Being the workhorse of our approach, the flux
conservation autoencoder is useful in exchanging information amongst
neighboring subdomains and from the boundary to the interior. This is similar
to traditional solvers where flux conservation is employed between
computational elements using numerical approximations for the same purpose.
Finally, it is important to note that the training portion in this method
corresponds to simply training $3$ autoencoders in use case 1 and $4$ in use
case 2. Autoencoders are an unsupervised ML technique to obtain lower-
dimensional encodings from higher-dimensional fields [6]. It is important to
note that we don’t explicitly train a model to learn a mapping between the
inputs (power map, HTC and thickness) and outputs (Temperature). Instead, we
solve a constrained equilibrium problem in the lower-dimensional latent space
during inference to predict the temperature for a fixed set of inputs, such as
power map, HTC and die thickness, which are specified by the user. The
implementation of this component is unique and differentiates our approach
from common ML methods.

Figure 3: Solution algorithm

###  2.1 Autoencoders in Thermal ML-Solver

In Figure 4, we describe the network architectures of autoencoders used in our
algorithm. The autoencoders in Fig. 4A are CNN-based encoder-decoder networks.
The encoder part of this network downsamples the subdomain field into a
corresponding lower-dimensional encoding. The decoder part of the network
upsamples the encoding back to the original input field. Separate autoencoders
are trained for temperature and powermap for both use cases. Additionally, an
autoencoder for distributed HTC is also trained in use case 2. On the other
hand, the flux conservation autoencoders are DNN-based encoder-decoder
networks. The temperature, powermap and system parameters on a 2-D stencil of
$5$ subdomains are encoded using the pretrained encoders and stacked together
as shown in Fig. 4B to form the input for this network. All the autoencoders
are trained in Tensorflow 2.2.0 using an Adam optimizer with mean absolute
error loss. We choose autoencoders over other compression techniques because
non-linear autoencoders have powerful generalization and denoising
capabilities with reasonable compression ratios [6]. This improves the
stability, robustness and ultimately the convergence of the constrained
equilibrium solve employed during the evaluation of our approach.

####  2.1.1 Specific architecture details:

Solution and Condition autoencoder: These are CNN based encoder-decoder
networks. The encoder part of the network has a series of 3 convolution layers
followed by max pooling. The number of filters in the convolution layers are
16, 32 and 64 respectively. The output of the convolution layer is flattened
and passed through 2 dense layers of size 1024 and 21, where 21 is the size of
the latent vector. In this work, the solution, power map and distributed HTC
fields are encoded to the same latent size. The decoder part of the network is
exactly symmetric to the encoder.

Flux conservation autoencoder: The flux conservation autoencoder has an input
size of $230$ in the case of constant HTC and die thickness and $330$ with
distributed HTC. This is a fully-connected encoder-decoder based network. The
encoder part of the network has 2 hidden layers with sizes 1024 and 512. The
latent vector for this network has a size of 35. The decoder part of the
network is exactly symmetric to the encoder.

Figure 4: A) Temperature and power map autoencoder architecture and B) Flux
conservation autoencoder architecture

##  3 Results and Discussion

In this section, we present results obtained from our approach for unseen test
cases across both the use cases and compare it with Ansys MAPDL. For the first
use case with constant HTC, we provide additional experiments to demonstrate
the generalizability of our approach for varying chip sizes as well as tile
size and out-of-range HTC values. Finally, we compare our solver with a
popular ML network, UNet [15], across all the experiments. It must be noted
that different instances of the solver are trained for the 2 use cases due to
the differences in HTC inputs.

###  3.1 Use case 1: Constant HTC

####  3.1.1 Unseen test cases

In this experiment we test the thermal ML-Solver for unseen, randomly sampled
powermaps, HTC and die thickness. The testing is carried out on $500$ samples
not included in the training set. The relative $L_{2}$ error, defined as
$\epsilon=\frac{L_{2}(T_{true}-T_{pred})}{L_{2}(T_{true})}$, is $0.078$. The
relative $L_{2}$ error serves as the most suitable metric in this case because
the system parameters can result in a wide range of chip temperature
magnitudes. Randomly selected samples are plotted in Fig. 5 and show that the
thermal ML-Solver matches well with results obtained from Ansys MAPDL.

Figure 5: Use case 1: Comparison of ML-Solver with Ansys Mechanical. Figure 6:
Use case 2: Comparison for distributed HTC on top and bottom die surfaces and
50 $\mu m$ tile size with Ansys Mechanical.

####  3.1.2 Generalization to larger chips

In this experiment, we test the generalizability of our approach to chips of
size $8000\mu m$ and $16000\mu m$. It should be noted that all models are
trained on a chip of physical size of $4000\mu m$. The testing is carried out
on $5$ testing samples for each chip size with randomly varying powermaps, HTC
and die thickness. The relative L-2 error for the $8000\mu m$ and $16000\mu m$
chips are $0.058$ and $0.089$ respectively. The prediction accuracy of our
thermal solver to larger chips may be attributed to the local learning
approaches coupled with iterative inferencing schemes. Additionally, in Fig. 7
we compare the contour plots of our approach with Ansys MAPDL for randomly
selected testing samples. It may be observed that for both chip sizes, the
thermal ML-Solver accurately models the temperature. Although not shown here,
we expect our solver to scale to even larger chip sizes because of the local
learning approach employed in our approach.

Figure 7: Use case 1: Comparison for bigger chip sizes 8000$\mu m$ (left) and
16000$\mu m$ (right) with Ansys Mechanical.

####  3.1.3 Generalization to a different tile size

Next, we demonstrate generalization to a different tile size of the powermap
equal to $250\mu m$. It should be noted that all models are trained on tile
sizes of $200\mu m$. The testing is carried out on $5$ testing samples for
each chip size with randomly varying powermaps, HTC and die thickness. The
relative L-2 error obtained is $0.14$. It may be observed from the contour
plots in Fig. 8 that the results match reasonably well with respect to Ansys
MAPDL.

####  3.1.4 Generalization to out-of-range HTCs

Figure 8: Use case 1: Comparison for higher HTC values (left) and bigger power
map tile sizes (right) with Ansys Mechanical.

In this section, we test the performance of the thermal ML-Solver in the
extrapolation regime of HTC. We generate $5$ samples for HTC $1e^{-4}$,
$2e^{-4}$, $3e^{-4}$, $4e^{-4}$ and $5e^{-4}$ for a die thickness of $150\mu
m$ and random power distribution. The average relative $L_{2}$ error is
$0.098$. It may be observed from the contour plots in Fig. 8 that the results
match reasonably well with respect to Ansys MAPDL and the error progressively
increases as we go further away from the training HTCs.

###  3.2 Use case 2: Distributed HTC and 50 $\mu m$ tile

####  3.2.1 Unseen test cases:

In this section, we test the thermal solver on $250$ testing samples with
unseen HTC and power map spatial distributions. Distributed HTCs can have a
greater impact on temperature distribution making it more challenging to
model. Additionally, the smaller tile sizes results in more local heating and
stiffer peaks in temperature which are also difficult to model. Despite the
challenges, our approach performs reasonably well with average relative L-2
errors of $0.045$ across all samples. It may also be observed from the contour
plots in Fig. 6 that the temperature predictions between our thermal ML-Solver
and Ansys MAPDL match well. Although not shown here, but our solver continues
to generalize to bigger chip sizes for this use case as well with the same
accuracy as the previous case.

###  3.3 Comparison with other ML baselines

Next, in table 1, we report the relative $L_{2}$ comparisons between our
approach and another popular approach, UNet [15] for all the test cases
considered in this paper. The significantly better accuracy of results
observed in table 1 demonstrates the superior generalization capability of our
approach. Moreover, since our approach operates like a solver it can
successfully scale to bigger chip sizes where traditional ML approaches
cannot.

Table 1: Comparison between thermal ML-Solver and UNet Experiment | Thermal ML-Solver | UNet
---|---|---
Section 3.1.1 | $0.042$ | $0.049$
Section 3.1.2 ($8000\mu m$) | $0.058$ | $1.05$
Section 3.1.2 ($16000\mu m$) | $0.089$ | $1.18$
Section 3.1.3 ($250\mu m$) | $0.14$ | $0.2$
Section 3.1.4 | $0.098$ | $0.35$

###  3.4 Computational time comparison

The computational time required by MAPDL on a $4000\mu m$ chip is about $30$
min on a single CPU. On the other hand, the thermal ML-Solver converges in
less than $10$ seconds on a single CPU. Moreover, since our approach has
similarities with traditional solvers it can be easily scaled to multiple CPUs
and GPUs. The savings in time increases proportionally as the chip size
increases. Finally, it must be noted that our approach employs iterative
evaluation and can be slower than black-box ML models but it compensates by
providing better accuracy, generalizability and scalability. For example, UNet
requires less than 1 second to compute the solution but has a degraded
accuracy and generalizability.

##  4 Conclusion

In this work, we introduce a thermal ML-Solver, which is based on the recently
proposed CoAEMLSim approach [13], to accurately predict temperature for
electronic chips across the high-dimensional system parameters in the form of
power map, HTC, die thickness etc. In this paper, we show two use cases with
constant and distributed HTCs, and demonstrate the accuracy of temperature
predictions with respect to Ansys MAPDL on unseen test cases. Moreover, we
provide additional experiments to test the generalizability of our approach
across varying chip sizes, tile sizes and out-of-range HTCs and demonstrate
our superior performance in comparison to a state-of the-art ML baseline.
Although the thermal solver demonstrated in this work is trained and tested
for powermap tile sizes of $200\mu m$ x $200\mu m$, the same approach can be
easily extended to smaller tile sizes ranging from of $1-10\mu m$. In future,
we would like to extend the approach to transient chip problems and chip
packages with geometric complexities and material variations.

## References

  * [1]
  * Cai et al. [2021] Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Em Karniadakis. 2021\.  Physics-informed neural networks for heat transfer problems.  _Journal of Heat Transfer_ 143, 6 (2021).
  * Chandra [2006] R Chandra. 2006\.  It’s Time To Consider Temperature Gradients In IC Design.  _Electronic Design_ (2006).
  * Chhabria et al. [2021] Vidya A Chhabria, Vipul Ahuja, Ashwath Prabhu, Nikhil Patil, Palkesh Jain, and Sachin S Sapatnekar. 2021\.  Thermal and IR drop analysis using convolutional encoder-decoder networks. In _Proceedings of the 26th Asia and South Pacific Design Automation Conference_. 690–696.
  * Emilio [2019] M Emilio. 2019.  Hybrid Chips may Solve Thermal Efficiency and Integration Challenges in 5G Mobile Devices.  _Power Electronics News_ (2019).
  * Goodfellow et al. [2016] Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio. 2016.  _Deep learning_. Vol. 1.  MIT press Cambridge.
  * He and Pathak [2020] Haiyang He and Jay Pathak. 2020.  An unsupervised learning approach to solving heat equations on chip based on auto encoder and image gradient.  _arXiv preprint arXiv:2007.09684_ (2020).
  * Hennigh et al. [2021] Oliver Hennigh, Susheela Narasimhan, Mohammad Amin Nabian, Akshay Subramaniam, Kaustubh Tangsali, Zhiwei Fang, Max Rietmann, Wonmin Byeon, and Sanjay Choudhry. 2021.  NVIDIA SimNet™: An AI-accelerated multi-physics simulation framework. In _International Conference on Computational Science_. Springer, 447–461.
  * Kumar et al. [2022] Akhilesh Kumar, Norman Chang, David Geb, Haiyang He, Stephen Pan, Jimin Wen, Saeed Asgari, Mehdi Abarham, and Chris Ortiz. 2022.  ML-based Fast On-Chip Transient Thermal Simulation for Heterogeneous 2.5 D/3D IC Designs. In _2022 International Symposium on VLSI Design, Automation and Test (VLSI-DAT)_. IEEE, 1–8.
  * Mutschler [2018] A Mutschler. 2018\.  New Thermal Issues Emerge.  _Semiconductor Engineering_ (2018).
  * Peach et al. [2019] C Peach, Y Zhang, S Makovejev, S Olsen, V Kilchytska, and J Raskin. 2019\.  Protecting AI Chips from Thermal Challenges during ATE Test.  In _Evaluation Engineering Magazine_.
  * Raissi and Karniadakis [2018] Maziar Raissi and George Em Karniadakis. 2018.  Hidden physics models: Machine learning of nonlinear partial differential equations.  _J. Comput. Phys._ 357 (2018), 125–141.
  * Ranade et al. [2021b] Rishikesh Ranade, Chris Hill, Haiyang He, Amir Maleki, Norman Chang, and Jay Pathak. 2021b.  A composable autoencoder-based iterative algorithm for accelerating numerical simulations.  _arXiv preprint arXiv:2110.03780_ (2021).
  * Ranade et al. [2021a] Rishikesh Ranade, Chris Hill, and Jay Pathak. 2021a.  DiscretizationNet: A machine-learning based solver for Navier–Stokes equations using finite volume discretization.  _Computer Methods in Applied Mechanics and Engineering_ 378 (2021), 113722\.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. 2015\.  U-net: Convolutional networks for biomedical image segmentation. In _International Conference on Medical image computing and computer-assisted intervention_. Springer, 234–241.
  * Stipsitz and Sanchis-Alepuz [2022] Monika Stipsitz and Hèlios Sanchis-Alepuz. 2022\.  Approximating the Steady-State Temperature of 3D Electronic Systems with Convolutional Neural Networks.  _Mathematical and Computational Applications_ 27, 1 (2022), 7\.
  * Sun et al. [2015] Yongsheng Sun, Canhui Zhan, Jianping Guo, Yiwei Fu, Guangming Li, and Jun Xia. 2015\.  Localized thermal effect of sub-16nm FinFET technologies and its impact on circuit reliability designs and methodologies. In _2015 IEEE International Reliability Physics Symposium_. IEEE, 3D–2.
  * Wen et al. [2020] Jimin Wen, Stephen Pan, Norman Chang, Wen-Tze Chuang, Wenbo Xia, Deqi Zhu, Akhilesh Kumar, En-Cih Yang, Karthik Srinivasan, and Ying-Shun Li. 2020\.  DNN-based fast static on-chip thermal solver. In _2020 36th Semiconductor Thermal Measurement, Modeling & Management Symposium (SEMI-THERM)_. IEEE, 65–75.
  * Zhong and Wong [2008] Yu Zhong and Martin DF Wong. 2008.  Thermal-aware IR drop analysis in large power grid. In _9th International Symposium on Quality Electronic Design (isqed 2008)_. IEEE, 194–199.