# Resilient Scheduling of Control Software Updates
in Radial Power Distribution Systems

Kin Cheong Sou and Henrik Sandberg K.C. Sou is with the Department of
Electrical Engineering, National Sun Yat-sen University, Taiwan
<EMAIL_ADDRESS>K.C. Sou was funded by the Ministry of Science and
Technology (MOST) of Taiwan: Optimal dispatch of flexibility resources in
power systems with high penetration of renewables, MOST
109-2221-E-110-016-MY2.H. Sandberg is with the Division of Decision and
Control Systems, EECS, KTH Royal Institute of Technology, Stockholm, Sweden
<EMAIL_ADDRESS>H. Sandberg was funded by the Swedish Energy Agency and ERA-Net
Smart Energy Systems (project RESili8, grant agreement No 883973).

###### Abstract

In response to newly found security vulnerabilities, or as part of a moving
target defense, a fast and safe control software update scheme for networked
control systems is highly desirable. We here develop such a scheme for
intelligent electronic devices (IEDs) in power distribution systems, which is
a solution to the so-called software update rollout problem. This problem
seeks to minimize the makespan of the software rollout, while guaranteeing
safety in voltage and current at all buses and lines despite possible worst-
case update failure where malfunctioning IEDs may inject harmful amounts of
power into the system. Based on the nonlinear DistFlow equations, we derive
linear relations relating software update decisions to the worst-case voltages
and currents, leading to a decision model both tractable and more accurate
than previous models based on the popular linearized DistFlow equations. Under
reasonable protection assumptions, the rollout problem can be formulated as a
vector bin packing problem and instances can be built and solved using
scalable computations. Using realistic benchmarks including one with 10,476
buses, we demonstrate that the proposed method can generate safe and effective
rollout schedules in real-time.

##  I INTRODUCTION

Critical infrastructures such as power and water distribution networks, and
cyber-physical systems (CPSs) in general, have in the past decade become
targets of cyberattacks [1, 2]. Simultaneously, cybersecurity features in
these systems are lacking. Control systems have been designed to meet safety
requirements, often meaning a preference for low complexity solutions that
minimize time delay and prevents the adoption of many standard computer
security features. It has now become urgent to develop security solutions that
take the special safety requirements of CPSs into account.

An essential difference between CPS and regular computer security is that
“software patching and frequent updates, are not well suited for control
systems” [3]. Software updates are problematic in that they may require a
complex reboot, which may come with an operational cost and a safety risk. It
is not uncommon that control software is updated infrequently, and systems may
run for long time even with known vulnerabilities. A practically relevant
research problem is to develop control algorithms that can be safely updated
in real time, to patch newly discovered security weaknesses or as part of
moving target or software rejuvenation defense [4, 5, 6].

This paper considers a control software update rollout problem for power
distribution systems introduced in [7], and extends the results in [8]. The
rollout problem seeks to arrange the software updates of IEDs (e.g., smart
inverters) into a minimum-time schedule, while guaranteeing power system
operational safety despite worst-case update failure. The problem is not
solvable using traditional maintenance planning techniques [9] due to its
intricate cyber-physical relationships involved. A central challenge here is
an accurate yet tractable relation between the software update decisions and
their worst-case consequences to the system states (e.g., voltages and
currents). We derive the desired relation using the _nonlinear_ DistFlow
equations [10, 11] as opposed to the popular linearized DistFlow equations
adopted in [7]. The nonlinear DistFlow equations more accurately describe
voltages and explicitly model line currents absent in the linearized
equations. We numerically demonstrate the ramifications of the DistFlow models
in maintaining safety standard. Compared with [8], this paper presents major
upgrades to the current and voltage safety limit constraints to reduce the
conservatism of the rollout problem. Our numerical case studies indicate that
the rollout schedules obtained in this paper are effective over a wide variety
of test cases. In addition, this paper features streamlined computation
procedures to set up and solve the rollout problem in real-time. These include
fixed point iterations to compute voltage and current bounds (required to
build a rollout problem instance) and a modified heuristic algorithm to solve
the associated bin-packing problem. These innovations eliminate the need for
mixed integer (bi)linear programming in [8], enabling practical real-time
large-scale rollout scheduling. For instance, our procedure can find a non-
trivial schedule in a 10,476-bus case in less than 3 seconds. This was
impossible in [8]. Our novel voltage and current bounds are potentially useful
in robust and contingency-based AC optimal power flow problems, which has
received significant interest in the past few years. See, for example, [12,
13, 14, 15, 16, 17].

Smart distribution grids play an essential role in society’s transition into a
net-zero energy system and thereby achieve highly set climate goals. The
transition requires integration of a wide variety of controllable energy
devices, both for consumption and for supply. The interface between power
supply (grid) and demand (consumer) is sometimes referred to as the _grid
edge_ [18]. The grid edge will turn distribution power grids, traditionally
not very automated, into complex distributed control systems. The many local
control loops in the grid edge simultaneously present an increased attack
surface, and will require resilient and systematic controller update/patching
schemes to handle newly discovered vulnerabilities or to respond to
contingencies [19, 20, 21], while accounting for possible adverse effects of
failed control updates.

This paper is organized as follows. Section II introduces notations,
distribution system model and its operational requirements. In Section III the
software update rollout problem is defined. The exact mathematical model is
unfit for real-time applications. Thus, a tractable approximation is derived.
Section IV summarizes the solution procedure for the rollout problem. Section
V presents case studies to demonstrate that the proposed procedure is time-
efficient and is able to deliver effective and safe rollout schedules for
large systems.

##  II System Modeling

###  II-A Mathematical Notations

We define the following notations: the vector ${\bf 1}$ (resp., ${\bf 0}$) is
the all-one (resp., all-zero) vector, and $e_{i}$ for $i=1,2,\ldots$ is the
$i^{\text{th}}$ unit vector. By default, the symbol $j$ means $\sqrt{-1}$. For
a set of positive integers $\mathcal{N}$ (details in Section II-B) and
$\mathcal{I}\subseteq\mathcal{N}$, the symbol ${\bf
1}_{\mathcal{I}}\in\\{0,1\\}^{\mathcal{N}}$ is the 0-1 binary indicator vector
with support $\mathcal{I}$. For any vector $v$, $\text{diag}(v)$ or $D_{v}$
are the diagonal matrix with the diagonal entries defined by $v$. For any
vector $v$ and scalar $a$ (positive or negative), the symbol $v^{a}$ denotes
the entry-wise exponentiation (i.e., the $i$-th entry of $v^{a}$ is
$v_{i}^{a}$). For any two vectors $x$ and $y$ of the same length, the symbol
$x\odot y$ denotes the entry-wise (Hadamard) product of $x$ and $y$. For any
two matrices $A$ and $B$, the symbol $A\otimes B$ denotes the Kronecker
product of $A$ and $B$.

###  II-B Power Distribution System Modeling

We consider single-phase radial power distribution systems. This can
potentially be interpreted as the positive sequence approximation of a three-
phase system. We assume that the system operates in steady-state and hence all
electrical quantities (e.g., current, voltage) can be represented using per
unit phasors. We use the following notations to describe the system:

  * •

$N$: number of non-reference buses (also number of lines)

  * •

$\mathcal{N}=\\{1,\ldots,N\\}$: the set of all non-reference buses. Bus 0 is
the reference bus (also called the slack bus).

  * •

$\mathcal{L}=\\{1,\ldots,N\\}$: the set of all lines. Each line has a
reference direction for line current and line power flow. Reference direction
points away from bus 0. A line is labeled by the bus that it points to (e.g.,
line $n$ points to bus $n$ and no other line points to $n$ due to radiality).

  * •

$\pi_{n}\in\mathcal{N}\cup\\{0\\}$: the “parent” bus of $n\in\mathcal{N}$ (so
that line $n$ goes from $\pi_{n}$ to $n$).

  * •

$d(n)\subseteq\mathcal{N}$: for any $n\in\mathcal{N}$, $d(n)$ denotes the set
of descendants of $n$ including $n$ itself. That is, $m\in d(n)$ if and only
if the (only) path from bus 0 to $m$ traverses $n$.

  * •

$\tilde{A}$: $N\times(N+1)$ line-bus incidence matrix. For
$(m,n)\in\mathcal{L}\times(\mathcal{N}\cup\\{0\\})$, $\tilde{A}_{mn}=1$ if
line $m$ leaves bus $n$, $\tilde{A}_{mn}=-1$ if line $m$ enters bus $n$ and
$\tilde{A}_{mn}=0$ otherwise.

  * •

$A$: $N\times N$ submatrix of $\tilde{A}$ with the first column $a_{0}$
removed (i.e., $\tilde{A}=\begin{bmatrix}a_{0}&A\end{bmatrix}$).

  * •

$\mathcal{D}$: $\\{0,1\\}^{N\times N}$ descendant matrix so that
$\mathcal{D}_{nm}=1$ if and only if $m\in d(n)$. It can be verified that
$\mathcal{D}^{\top}=-A^{-1}$ when row $k$ of $A$ corresponds to line $k$ for
all $k\in\mathcal{L}$.

  * •

All shunt elements are ignored and $z=r+jx$ is an $N$-vector of line series
impedance; $r>{\bf 0}$ is vector of series resistances and $x>{\bf 0}$ is
vector of series reactances.

  * •

$D_{r}$, $D_{x}$: $N\times N$ positive diagonal matrices with
$D_{r}=\text{diag}(r)$ and $D_{x}=\text{diag}(x)$. Also,
$D_{z}:=D_{r}+jD_{x}$.

  * •

$P$, $Q$, $S$: $N$-vectors of active, reactive and apparent power flows for
all lines. $P$, $Q$ and $S$ are defined at the sending ends of the lines
(e.g., $P_{n}$ is defined at bus $\pi_{n}$).

  * •

$p$, $q$, $s$: $N$-vectors of net injections of active power, reactive power,
apparent power for non-reference buses.

  * •

$p=p^{G}-p^{L}$, $q=q^{G}-q^{L}$ with $p^{G}$, $q^{G}$, $p^{L}$, $q^{L}$
vectors of active power generation, reactive power generation, active power
load, reactive power load for non-reference buses.

  * •

$\ell$: $N$-vector of squared magnitude line currents.

  * •

$v$: $N$-vector of magnitude voltage of non-reference buses.

  * •

$\nu$: $N$-vector of squared magnitude voltage (i.e., $\nu=v^{2}$). At bus 0,
squared magnitude voltage $\nu_{0}$ is constant.

For all vectors above, subscript $n$ means the $n^{\text{th}}$ entry of the
vector. For instance, $p^{G}_{n}$ is the active power generation injection to
bus $n\in\mathcal{N}$, while $\ell_{n}$ is the squared current for line
$n\in\mathcal{L}$. Fig. 1 illustrates a segment of the distribution system
with the relevant electrical quantities.

Figure 1: A segment of distribution system with the relevant quantities

For a radial distribution system with shunt elements ignored, the electrical
quantities $\nu,\ell,P,Q,p,q$ are related by the (nonlinear) DistFlow equation
[10, 11] as follows:

$\displaystyle p$ $\displaystyle=A^{\top}P+D_{r}\ell$ (1a) $\displaystyle q$
$\displaystyle=A^{\top}Q+D_{x}\ell$ (1b) $\displaystyle A\nu+\nu_{0}a_{0}$
$\displaystyle=2D_{r}P+2D_{x}Q-(D_{r}^{2}+D_{x}^{2})\ell$ (1c)
$\displaystyle\ell_{n}$
$\displaystyle=\frac{P_{n}^{2}+Q_{n}^{2}}{\nu_{\pi_{n}}},\quad\forall
n=1,2,\ldots,N$ (1d)

Note that $A^{-1}=-\mathcal{D}^{\top}$ and $A^{-1}a_{0}=-{\bf 1}$ (since
$\tilde{A}{\bf 1}={\bf 0}$). The variables $P$ and $Q$ can be eliminated from
(1) to obtain

$\displaystyle\nu$ $\displaystyle=\nu_{0}{\bf 1}+2Rp+2Xq+M\ell$ (2a)
$\displaystyle\ell$
$\displaystyle=(\nu)^{-1}\odot\big{|}\mathcal{D}(p+jq)-(\mathcal{D}-I)D_{z}\ell\;\big{|}^{2}$
(2b)

where

$\displaystyle\begin{split}R&=\mathcal{D}^{\top}D_{r}\mathcal{D}\\\
X&=\mathcal{D}^{\top}D_{x}\mathcal{D}\\\
M&=\mathcal{D}^{\top}D_{r}(I-2\mathcal{D})D_{r}+\mathcal{D}^{\top}D_{x}(I-2\mathcal{D})D_{x}\end{split}$
(3)

Equation (2) implies that $\nu$ and $\ell$ are functions of net power
injections $(p,q)$, motivating the shorthand $\nu=\nu(p,q)$ and
$\ell=\ell(p,q)$. Also, since $\mathcal{D}\in\\{0,1\\}^{N\times N}$ with unit
diagonal entries and $r$, $x$ are positive, (3) implies that all entries of
$R$ and $X$ are nonnegative, while all entries of $M$ are non-positive.

###  II-C Distribution System Setting and Grid Code

We assume that every non-reference bus in the distribution system is equipped
with a smart inverter to control and monitor its distributed generation. The
inverter-interfaced distributed generation is $p^{G}_{n}\geq 0$ for bus
$n\in\mathcal{N}$. In addition, the inverter can provide reactive power
support with reactive power injection $q^{G}_{n}$, which is not sign
restricted. We assume the rating of the inverter at bus $n$ is $C_{n}$
restricting its apparent power injection. If bus $n$ has no generation, then
$C_{n}=0$. The vector of all inverter ratings is denoted by $C$. Thus, the
active and reactive power injections of the inverters satisfy

$p^{G}_{n}\geq 0,\quad\big{|}p^{G}_{n}+jq^{G}_{n}\big{|}\leq
C_{n},\quad\forall n\in\mathcal{N}$ (4)

It is assumed that the operator knows the inverter generation set points
denoted by $\hat{p}_{n}^{G}$ and $\hat{q}_{n}^{G}$ and the load estimates
denoted by $\hat{p}^{L}_{n}$ and $\hat{q}^{L}_{n}$. In normal operating
conditions, inverter set points $\hat{p}^{G}$ and $\hat{q}^{G}$ are chosen so
that the grid code is satisfied. In this paper, the grid code specifies that
the steady-state $\nu$ and $\ell$ should satisfy the following safety limits:

$\displaystyle\underline{\nu}\leq\nu\leq\overline{\nu},\quad\text{with
given}\;\;\underline{\nu}\in\mathbb{R}^{N}\;\text{and
given}\;\overline{\nu}\in\mathbb{R}^{N},$ (5a) $\displaystyle
0\leq\ell\leq\overline{\ell},\quad\text{with
given}\;\;\overline{\ell}\in\mathbb{R}^{N}$ (5b)

We note that $\ell\geq 0$ is not explicitly enforced as it is implied by (5a)
and (1d). Typical values of the safety limits are: 0.9 pu for
$(\underline{\nu})^{1/2}$, 1.1 pu for $(\overline{\nu})^{1/2}$ and a few
hundred amperes (e.g., 600 A) for $(\overline{\ell})^{1/2}$.

##  III Problem Statement and Formulation

###  III-A Rollout Problem Setting and Statement

As explained in the introduction, we envision a scenario where, due to
security or operational considerations, the operator needs to remotely update
the software or firmware of all inverters through some communication network.
To minimize possible disruption, the operator wishes to schedule the updates
to finish as soon as possible (e.g., all at once if possible). However, some
or all updates may fail. Software update failure for the inverter at a bus is
modeled by the condition that, instead of following its power set point
command, the inverter may inject an uncontrollable amount of power subjected
to its physical limit in (4). We assume that each bus is equipped with a relay
to promptly detect out-of-range voltage or current (usually within a few
milliseconds). Then the breaker or other protection will act to isolate the
faulty inverter from the rest of the system and the update for the inverter is
rolled back. We assume that $\Delta\tau$ is the universal fault clearing time
for the system. In other words, if an inverter experiences software update
failure, it may inject uncontrolled amount of power for up to $\Delta\tau$
seconds, and then its power output is back to normal.

It is not desirable to schedule too many updates at the same time lest the
updates may fail simultaneously and inject a dangerous level of uncontrolled
power causing power outage or even risk to human lives. Thus, the operator is
faced with an update scheduling problem where s/he seeks to minimize the
makespan of the entire process while guaranteeing the safe operation of the
system despite the worst-case update failure due to the schedule. This problem
is referred to as the _software update rollout problem_ or the _rollout
problem_ for short.

Safe operation means that the steady-state squared voltage $\nu$ and steady-
state squared line current $\ell$ satisfy the grid code constraints in (5).
The transient voltage and current from a normal operating condition to the
situation after update failure may be different from their respective steady-
state values. However, we assume that the system possesses sufficient ride-
through capability so that if (5) is satisfied in steady states then the
transient voltage and current will satisfy their respective transient safety
requirements.

Formally, the safety requirement of the rollout problem can be stated as
follows: For any $\mathcal{I}\subseteq\mathcal{N}$, define
$\mathcal{S}(\mathcal{I})$ by

$\begin{split}\mathcal{S}(\mathcal{I}):=&\big{\\{}(p,q)\mid
p=p^{G}-\hat{p}^{L},\;q=q^{G}-\hat{q}^{L}\\\
&\big{|}p^{G}_{n}+jq^{G}_{n}\big{|}\leq C_{n},\;p^{G}_{n}\geq 0,\;\text{if
$n\in\mathcal{I}$}\\\
&p^{G}_{n}=\hat{p}^{G}_{n},\;q^{G}_{n}=\hat{q}^{G}_{n},\;\text{if
$n\in\mathcal{N}\setminus\mathcal{I}$}\big{\\}}\end{split}$ (6)

as the set of all possible $(p,q)$ due to update failures at buses in
$\mathcal{I}\subseteq\mathcal{N}$. Then, it is required that all solutions
$\nu(p,q)$ and $\ell(p,q)$ of (1) due to $(p,q)\in\mathcal{S}(\mathcal{I})$
satisfy the safety limits in (5). Subsequently, we derive linear relationships
between ${\bf 1}_{\mathcal{I}}$ (i.e., the decision indicator vector of the
rollout problem) and the worst-case $\nu(p,q)$ and $\ell(p,q)$ for
$(p,q)\in\mathcal{S}(\mathcal{I})$ (used to describe the safety constraints in
the rollout problem).

###  III-B Universal voltage lower bound and current upper bound

To derive the safety constraints of the rollout problem, we need bounds of
$\nu$ and $\ell$ as simple functions of $(p,q)$. Since $\ell\geq 0$ and $M$
has non-positive entries in (2), an immediate upper bound of $\nu$ arises when
the term $M\ell$ in (2a) is ignored. This will be discussed in Section III-D.
On the other hand, lower bound of $\nu$ and upper bound of $\ell$ are
intertwined. To obtain the bounds, for any $n\in\mathcal{N}$ we denote squared
voltage lower bound $\nu^{L}_{n}$ and squared current upper bound
$\ell^{U}_{n}$ as

$\displaystyle\nu^{L}_{n}$
$\displaystyle:=\min\limits_{p,q,P,Q,\nu,\ell}\;\nu_{n}\quad\text{subject to
\eqref{eqn:DF} and \eqref{eqn:inverter_constraints}}$ (7a)
$\displaystyle\ell^{U}_{n}$
$\displaystyle:=\max\limits_{p,q,P,Q,\nu,\ell}\;\ell_{n}\quad\text{subject to
\eqref{eqn:DF} and \eqref{eqn:inverter_constraints}}$ (7b)

Then, Sections III-C and III-E discuss how $\nu^{L}$ and $\ell^{U}$ can be
used in (2) to obtain the desired voltage and current bounds with respect to
$(p,q)$. Problem (7) can be solved as a non-convex bilinear program (due to
(1d)). However, this may not be suitable for real-time large-scale
applications. Instead, we propose to quickly estimate lower bounds of
$\nu^{L}$ and upper bounds of $\ell^{U}$ using the fixed point iterations in
Algorithm 1.

Algorithm 1 Fixed point iterations to bound $\nu^{L}$ and $\ell^{U}$

1:Matrices $\mathcal{D}$ and $D_{z}$ in Section II-B, vectors $\hat{p}^{L}$,
$\hat{q}^{L}$, $C$ in Section II-C, $R$, $X$, $M$ in (3) and tolerance
$\epsilon>0$

2:$\hat{\nu}^{L}\leq\nu^{L}$ and $\hat{\ell}^{U}\geq\ell^{U}$

3:Denote $a:=\mathcal{D}\hat{p}^{L}$, $b:=\mathcal{D}\hat{q}^{L}$,
$r:=\mathcal{D}C$, $\theta\in\mathbb{R}^{N}$ s.t.
$\displaystyle\theta_{n}:=\begin{cases}\tan^{-1}(b_{n}/a_{n}),&\text{if
$a_{n}\neq 0$ or $b_{n}\neq 0$}\\\ 0,&\text{if $a_{n}=b_{n}=0$}\end{cases}$

4:Define  $\displaystyle\bar{S}:=\max\Big{\\{}$
$\displaystyle|a\\!+\\!j(b\\!+\\!r)|,|a\\!+\\!j(b\\!-\\!r)|,|a\\!+\\!jb\\!-\\!r\\!\odot\\!e^{j\theta}|\Big{\\}}$
where “max” is interpreted row-wise

5:Initialize $v^{0}\in\mathbb{R}^{N}$ and $i^{0}\in\mathbb{R}^{N}$ such that
$v^{0}\leq(\nu^{L})^{1/2}$ and $i^{0}\geq(i^{U})^{1/2}$ and $k\leftarrow 0$

6:while
$\|v^{k}-v^{k-1}\|>\epsilon\;\;\text{or}\;\;\|i^{k}-i^{k-1}\|>\epsilon$ do

7: $v^{k+1}=\big{(}\nu_{0}{\bf
1}\\!-\\!2R\hat{p}^{L}-2X(\hat{q}^{L}\\!+\\!C)\\!+\\!M(i^{k})^{2}\big{)}^{1/2}$

8:
$i^{k+1}=\big{(}v^{k})^{-1}\odot(\bar{S}+|(\mathcal{D}-I)D_{z}(i^{k})^{2}|\big{)}$

9: $k\leftarrow k+1$

10:end while

11:Return $\hat{\nu}^{L}=(v^{k})^{2}$ and $\hat{\ell}^{U}=(i^{k})^{2}$

Three questions pertain to Algorithm 1: (a) Does Algorithm 1 converge? (b) Is
it true that $\hat{\nu}^{L}\approx\nu^{L}$ and
$\hat{\ell}^{U}\approx\ell^{U}$? (c) Is it true that
$\hat{\nu}^{L}\leq\nu^{L}$ and $\hat{\ell}^{U}\geq\ell^{U}$? Question (a)
regarding convergence can be verified online. Further, our extensive
simulations indicate quick convergence as long as the initial guess $i^{0}$ is
not too large (so $v^{k}$ stays real-valued). As explained in Section V-B,
question (b) regarding approximation quality can also be checked online.
However, question (c) regarding bounding property cannot be verified
numerically and it must be proven analytically. Next, we provide an intuition
to justify (b) and (c). Then, a statement is presented to prove (c).

Suppose the while-loop of line 4-8 in Algorithm 1 converge with limits
$\ell^{\star}:=(\lim\limits_{k\to\infty}i^{k})^{2}$ and
$\nu^{\star}:=(\lim\limits_{k\to\infty}v^{k})^{2}$. Then, $i^{\star}$ and
$\nu^{\star}$ are characterized by

$\displaystyle\nu^{\star}$ $\displaystyle=\nu_{0}{\bf
1}-2R\hat{p}^{L}-2X(\hat{q}^{L}+C)+M\ell^{\star}$ (8a)
$\displaystyle\ell^{\star}$
$\displaystyle={(\nu^{\star})}^{-1}\odot\big{(}\bar{S}+|(\mathcal{D}-I)D_{z}\ell^{\star}|\big{)}^{2}$
(8b)

which resembles (2), except that the right-hand-sides of (8) are optimized
with respect to $(p,q)$ subject to (4). Indeed, due to the fact that $R$ and
$X$ in (3) have nonnegative entries and Lemma 1 in Appendix A the following
holds: for any $n\in\mathcal{N}$, entries $\nu^{\star}_{n}$ in (8a) and
$\ell^{\star}_{n}$ in (8b) can be rewritten as

$\displaystyle\begin{split}\nu^{\star}_{n}&=\underset{(p,q)\;\text{in}\;\text{\eqref{eqn:inverter_constraints}}}{\min}e_{n}^{\top}\left(\nu_{0}{\bf
1}+2Rp+2Xq+M\ell^{\star}\right)\\\
\ell^{\star}_{n}&=\underset{(p,q)\;\text{in}\;\text{\eqref{eqn:inverter_constraints}}}{\max}e_{n}^{\top}\big{(}{(\nu^{\star})}^{-1}\odot\big{(}|\mathcal{D}(p+jq)|\\!+\\!|(\mathcal{D}-I)D_{z}\ell^{\star}|\big{)}^{2}\big{)}\\\
&\geq\underset{(p,q)\;\text{in}\;\text{\eqref{eqn:inverter_constraints}}}{\max}e_{n}^{\top}\big{(}{(\nu^{\star})}^{-1}\odot\big{|}\mathcal{D}(p+jq)-(\mathcal{D}-I)D_{z}\ell^{\star}\big{|}^{2}\big{)}\end{split}$
(9)

where “$(p,q)$ in (4)” means $(p,q)$ satisfying (4) in (9). The relations in
(9) have two implications. Firstly, the resemblance of (9) to (2) (equivalent
to (1)) and the optimization suggest that $\nu^{\star}\approx\nu^{L}$ and
$\ell^{\star}\approx\ell^{U}$ (hence $\hat{\nu}^{L}\approx\nu^{L}$ and
$\hat{\ell}^{U}\approx\ell^{U}$). Secondly, the optimization in (9) suggests
the inequalities $\nu^{\star}\leq\nu^{L}$ and $\ell^{\star}\geq\ell^{U}$ and
hence $\hat{\nu}^{L}\leq\nu^{L}$ and $\hat{\ell}^{U}\geq\ell^{U}$. Indeed,
this can be verified by the following statement, whose proof can be found in
Appendix B.

######  Proposition 1.

For Algorithm 1, it holds that $v^{k}\leq(\nu^{L})^{1/2}$ and
$i^{k}\geq(\ell^{U})^{1/2}$ for all $k\geq 0$.

In particular, Proposition 1 specifies that $(\hat{\nu}^{L},\hat{\ell}^{U})$
as returned by Algorithm 1 satisfy the desired bounding inequalities.
Proposition 1 also guarantees that, in case Algorithm 1 does not converge
(though we never encounter this in our extensive simulations), any premature
solution from Algorithm 1 still satisfies the bounding inequalities. In
practice, Algorithm 1 converges quickly even if the initialization conditions
$v^{0}\leq(\nu^{L})^{1/2}$ and $i^{0}\geq(i^{U})^{1/2}$ are not followed. As
is standard in load flow analysis, we adopt a “flat start” of $v^{0}={\bf 1}$
and $i^{0}={\bf 0}$ in here. For convenience, the shorthand
$\hat{v}^{L}:=(\hat{\nu}^{L})^{1/2}$ and $\hat{i}^{U}:=(\hat{\ell}^{U})^{1/2}$
will be used subsequently.

###  III-C Expressions Related to Line Current Upper Limit

This section develops two separate results – (a) an upper bound of line
current affinely dependent on ${\bf 1}_{\mathcal{I}}$ and (b) given any line
current upper limit $\bar{\ell}\in\mathbb{R}^{N}$ a set of linear inequalities
on ${\bf 1}_{\mathcal{I}}$ to ensure that $\ell(p,q)\leq\bar{\ell}$ for all
$(p,q)\in\mathcal{S}(\mathcal{I})$.

Affine current upper bound: Let $\ell=\ell(p,q)$ from (2). Then, $\ell^{1/2}$
is the vector of magnitude line currents and hence by Kirchhoff’s current law

$\displaystyle\sqrt{\ell_{n}}=\left|\sum\limits_{m\in
d(n)}i^{\prime}_{m}\right|\leq\\!\\!\\!\sum\limits_{m\in
d(n)}\\!\\!\\!\left|i^{\prime}_{m}\right|\leq\\!\\!\\!\sum\limits_{m\in
d(n)}\\!\\!\\!\frac{|p_{m}+jq_{m}|}{\hat{v}^{L}_{m}},\;\forall
n\in\mathcal{L}$ (10)

where $i^{\prime}_{m}$ is the net current injection at bus $m$, $d(n)$ is the
set of ($n$-inclusive) descendants of $n$ defined in Section II-B and
$\hat{v}^{L}$ is the voltage lower bound by Algorithm 1 (so that
$\hat{\nu}^{L}\leq\nu^{L}$ in (7a)). In (10), the last inequality holds
because $|p_{m}+jq_{m}|=|v_{m}||i^{\prime}_{m}|$ and hence
$|p_{m}+jq_{m}|\geq\hat{v}^{L}_{m}|i^{\prime}_{m}|$. Given power injection
$(p,q)$, (10) provides an upper bound on $\ell(p,q)^{1/2}$. Let
$\mathcal{I}\subseteq\mathcal{N}$ denote the set of buses with software update
failure implying that $(p,q)\in\mathcal{S}(\mathcal{I})$ in (6). Then, the
maximum line current $\sqrt{\ell_{n}}$ subject to $\mathcal{I}$ is upper
bounded by

$\displaystyle\begin{split}&\max\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}\sum\limits_{m\in
d(n)}\frac{|p_{m}+jq_{m}|}{\hat{v}^{L}_{m}}\\\ =&\sum\limits_{m\in
d(n)\setminus\mathcal{I}}\frac{|(\hat{p}^{G}_{m}-\hat{p}^{L}_{m})+j(\hat{q}^{G}_{m}-\hat{q}^{L}_{m})|}{\hat{v}^{L}_{m}}+\sum\limits_{m\in
d(n)\cap\mathcal{I}}\frac{{|s^{\star}|}_{m}}{\hat{v}^{L}_{m}}\end{split}$ (11)

where

$\displaystyle\begin{split}{|s^{\star}|}_{m}=&\max\limits_{p^{G}_{m},q^{G}_{m}}\quad|(\hat{p}^{L}_{m}+j\hat{q}^{L}_{m})-(p^{G}_{m}+jq^{G}_{m})|\\\
&\text{subject to}\quad p^{G}_{m}\geq 0,\;\;|p^{G}_{m}+jq^{G}_{m}|\leq
C_{m}\end{split}$ (12)

By Lemma 1 in Appendix A, ${|s^{\star}|}_{m}$ in (12) is equal to

$\displaystyle\begin{split}{|s^{\star}|}_{m}=&\max\left\\{|\hat{p}^{L}_{m}+j(\hat{q}^{L}_{m}+C_{m})|,|\hat{p}^{L}_{m}+j(\hat{q}^{L}_{m}-C_{m})|,\right.\\\
&\left.|\hat{p}^{L}_{m}+j\hat{q}^{L}_{m}-C_{m}(\cos(\theta_{sm})+j\sin(\theta_{sm})|\right\\}\end{split}$
(13)

where $\theta_{sm}=\tan^{-1}(\hat{q}^{L}_{m}/\hat{p}^{L}_{m})$ if at least one
of $\hat{p}^{L}_{m}$ and $\hat{q}^{L}_{m}$ is nonzero, and $\theta_{sm}=0$
otherwise. Define

$\displaystyle\begin{split}&\text{$|s^{\star}|\in\mathbb{R}^{N}$ with entries
specified in \eqref{eqn:sm}}\\\
&|s^{\text{nom}}|:=\left|(\hat{p}^{G}-\hat{p}^{L})+j(\hat{q}^{G}-\hat{q}^{L})\right|\end{split}$
(14)

Then, the current upper bound described by (11), (13) and (14) can be written
compactly as

$\displaystyle\begin{split}\ell^{1/2}\\!&\leq\mathcal{D}\Big{(}(\hat{v}^{L})^{-1}\\!\odot\\!|s^{\text{nom}}|\\!+\\!\text{diag}\big{(}(\hat{v}^{L})^{-1}\\!\odot(|s^{\star}|\\!-\\!|s^{\text{nom}}|)\big{)}{\bf
1}_{\mathcal{I}}\Big{)}\\\ &:=I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}}\end{split}$ (15)

In summary, if $\mathcal{I}\subseteq\mathcal{N}$ denotes the set of bus(es)
with software update(s), then despite the worst-case update failure, the line
current vector is upper bounded by the expression in (15).

Linear inequalities guaranteeing current upper limit: For any net power
injections $(p,q)$, (2b) specifies that

$\displaystyle\sqrt{\ell_{n}}$
$\displaystyle=(\nu_{n})^{-1/2}e_{n}^{\top}\left|\mathcal{D}(p+jq)-(\mathcal{D}-I)D_{z}\ell\right|$
$\displaystyle\leq(\hat{v}^{L}_{n})^{-1}e_{n}^{\top}\big{(}|\mathcal{D}(p\\!+\\!jq)|\\!+\\!|(\mathcal{D}\\!-\\!I)D_{z}(\hat{i}^{U}\\!\odot\ell^{1/2})|\big{)}$

where $\mathcal{D}$ and $D_{z}$ are defined in Section II-B, and
$\hat{v}^{L}=(\hat{\nu}^{L})^{1/2}$, $\hat{i}^{U}=(\hat{\ell}^{U})^{1/2}$ from
Algorithm 1. The inequality above is due to the following four facts: (a)
$\hat{v}^{L}\leq(\nu(p,q))^{1/2}$, (b) triangular inequality, (c) ${\bf
0}\leq\ell(p,q)=\ell^{1/2}\odot\ell^{1/2}\leq\hat{i}^{U}\odot\ell^{1/2}$ and
(d) entries of $(\mathcal{D}-I)$, $D_{r}$, $D_{x}$ are nonnegative. Let
$\mathcal{I}\subseteq\mathcal{N}$ denote the set of buses with software update
failure implying that $(p,q)\in\mathcal{S}(\mathcal{I})$ in (6). Then, for any
$n\in\mathcal{N}$

$\displaystyle\begin{split}\max\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}\sqrt{\ell_{n}}\leq&\;(\hat{v}^{L}_{n})^{-1}\max\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}e_{n}^{\top}\left|\mathcal{D}(p+jq)\right|\\\
&+(\hat{v}^{L}_{n})^{-1}e_{n}^{\top}|(\mathcal{D}-I)D_{z}(\hat{i}^{U}\odot\ell^{1/2})|\end{split}$
(16)

The two terms in the right-hand-side of (16) can be analyzed separately. By
Lemma 2 in Appendix A the first term is
$(\hat{v}^{L}_{n})^{-1}\max\\{|a+j(b+r)|,|a+j(b-r)|,|a+jb-re^{j\theta}|\\}$
with $a$, $b$, $r$, $\theta$ defined in (49) in Lemma 2. Due to $\theta$, the
expression is a nontrivial function of ${\bf 1}_{\mathcal{I}}$ and hence we
consider its simplifying upper bound $(\hat{v}^{L}_{n})^{-1}(|a+jb|+r)$, which
is written as

$\displaystyle\begin{split}&(\hat{v}^{L}_{n})^{-1}\Big{|}e_{n}^{\top}\mathcal{D}(\hat{p}^{L}\\!-\\!\hat{p}^{G}\\!+\\!j(\hat{q}^{L}\\!-\\!\hat{q}^{G}))\\!+\\!e_{n}^{\top}\mathcal{D}D_{(\hat{p}^{G}+j\hat{q}^{G})}{\bf
1}_{\mathcal{I}}\Big{|}\\\
&+(\hat{v}^{L}_{n})^{-1}e_{n}^{\top}\mathcal{D}D_{C}{\bf
1}_{\mathcal{I}}\end{split}$ (17)

Because of the absolute value, (17) is not affine with respect to ${\bf
1}_{\mathcal{I}}$. Thus, it is further upper bounded by polygonal
approximation. In particular, for any $x+jy\in\mathbb{C}$ it holds that

$\displaystyle\begin{split}&|x+jy|\\\
\leq&\frac{1}{\cos(\pi/N^{\prime})}\max\limits_{1\leq k\leq
N^{\prime}}\big{\\{}\cos(\frac{(2k-1)\pi}{N^{\prime}})x\\!+\\!\sin(\frac{(2k-1)\pi}{N^{\prime}})y\big{\\}}\end{split}$
(18)

for any integer $N^{\prime}\geq 4$. For example, if $N^{\prime}=4$ then

$\displaystyle|x+jy|\leq\max\\{x+y,x-y,-x+y,-x-y\\}=|x|+|y|$

Therefore,

$\displaystyle\begin{split}&\text{\eqref{eqn:iub_inexact_max}}\\\
\leq&\;(\hat{v}^{L}_{n})^{-1}\frac{1}{\cos(\pi/N^{\prime})}\\!\max\limits_{1\leq
k\leq
N^{\prime}}\\!\Big{\\{}\cos(\frac{(2k\\!-\\!1)\pi}{N^{\prime}})\big{(}e_{n}^{\top}\mathcal{D}(\hat{p}^{L}\\!-\\!\hat{p}^{G})\\\
&+e_{n}^{\top}\mathcal{D}D_{\hat{p}^{G}}{\bf
1}_{\mathcal{I}}\big{)}+\sin(\frac{(2k-1)\pi}{N^{\prime}})\big{(}e_{n}^{\top}\mathcal{D}(\hat{q}^{L}-\hat{q}^{G})\\\
&+e_{n}^{\top}\mathcal{D}D_{\hat{q}^{G}}{\bf
1}_{\mathcal{I}}\big{)}\Big{\\}}+(\hat{v}^{L}_{n})^{-1}e_{n}^{\top}\mathcal{D}D_{C}{\bf
1}_{\mathcal{I}}\end{split}$ (19)

Now we consider the second term in the right-hand-side of (16). Using (15) and
the fact that all entries of $(\mathcal{D}-I)$, $D_{r}$, $D_{x}$ are
nonnegative, the term can be upper bounded by

$\displaystyle(\hat{v}^{L}_{n})^{-1}\Big{|}e_{n}^{\top}(\mathcal{D}-I)\text{diag}(z\odot\hat{i}^{U})(I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}})\Big{|}$ (20)

Using polygonal approximation in (18) with integer $N^{\prime\prime}\geq 4$,
it holds that

$\displaystyle\begin{split}&\text{\eqref{eqn:iub_inexact_loss}}\\\
\leq&\;(\hat{v}^{L}_{n})^{-1}\frac{1}{\cos\left(\pi/N^{\prime\prime}\right)}\max\limits_{1\leq
k\leq
N^{\prime\prime}}\Big{\\{}\cos(\frac{(2k\\!-\\!1)\pi}{N^{\prime\prime}})e_{n}^{\top}(\mathcal{D}-I)\\\
&\times\text{diag}(r\odot\hat{i}^{U})(I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}})+\sin(\frac{(2k-1)\pi}{N^{\prime\prime}})\\\ &\times
e_{n}^{\top}(\mathcal{D}-I)\text{diag}(x\odot\hat{i}^{U})(I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}})\Big{\\}}\end{split}$ (21)

Hence, the sum of the RHS of (19) and (21) is an upper bound of
$\sqrt{\ell_{n}(p,q)}$ for all $(p,q)\in\mathcal{S}(\mathcal{I})$. Thus,
$\forall(p,q)\in\mathcal{S}(\mathcal{I})$

$\displaystyle\text{RHS \eqref{eqn:iub_inexact_max_pa}}+\text{RHS
\eqref{eqn:iub_inexact_loss_pa}}\leq(\bar{\ell}_{n})^{1/2}\\!\implies\\!\sqrt{\ell_{n}}\leq(\bar{\ell}_{n})^{1/2}$
(22)

Utilizing the fact that for any $x\in\mathbb{R}^{N^{\prime}}$,
$y\in\mathbb{R}^{N^{\prime\prime}}$, $z\in\mathbb{R}$

$\displaystyle\max\limits_{i}x_{i}+\max\limits_{j}y_{j}\leq z\iff{\bf
1}_{N^{\prime\prime}}\otimes x+y\otimes{\bf 1}_{N^{\prime}}\leq z{\bf
1}_{N^{\prime}N^{\prime\prime}}$

with ${\bf 1}_{N^{\prime}}$, ${\bf 1}_{N^{\prime\prime}}$, ${\bf
1}_{N^{\prime}N^{\prime\prime}}$ denoting all-one vectors of appropriate
dimensions, the condition in (22) can be “vectorized” as a set of linear
constraints with respect to ${\bf 1}_{\mathcal{I}}$ to enforce
$\sqrt{\ell_{n}}\leq(\bar{\ell}_{n})^{1/2}$ for a specific $n\in\mathcal{N}$.
By “vectorizing” the conditions in (22) for all $n$, the desired linear
constraints can be summarized as

$I^{\text{iub2}}+W^{\text{iub2}}{\bf 1}_{\mathcal{I}}\leq{\bf
1}_{N^{\prime\prime}}\otimes({\bf 1}_{N^{\prime}}\otimes(\bar{\ell})^{1/2})$
(23)

where

$\displaystyle\begin{split}I^{\text{iub2}}&={\bf
1}_{N^{\prime\prime}}\otimes\Big{(}c^{\prime}\otimes\big{(}(\hat{v}^{L})^{-1}\odot\mathcal{D}(\hat{p}^{L}-\hat{p}^{G})\big{)}\\\
&+s^{\prime}\otimes\big{(}(\hat{v}^{L})^{-1}\odot\mathcal{D}(\hat{q}^{L}-\hat{q}^{G})\big{)}\Big{)}\\\
&+c^{\prime\prime}\\!\otimes\\!{\bf
1}_{N^{\prime}}\otimes\big{(}{(\hat{v}^{L})}^{-1}\odot(\mathcal{D}\\!-\\!I)\text{diag}(r\odot\hat{i}^{U})I^{\text{iub1}}\big{)}\\\
&+s^{\prime\prime}\\!\otimes\\!{\bf
1}_{N^{\prime}}\otimes\big{(}{(\hat{v}^{L})}^{-1}\odot(\mathcal{D}\\!-\\!I)\text{diag}(x\odot\hat{i}^{U})I^{\text{iub1}}\big{)}\\\
W^{\text{iub2}}&={\bf 1}_{N^{\prime\prime}}\otimes\Big{(}c^{\prime}\otimes
D_{\hat{v}^{L}}^{-1}\mathcal{D}D_{\hat{p}^{G}}+s^{\prime}\otimes
D_{\hat{v}^{L}}^{-1}\mathcal{D}D_{\hat{q}^{G}}\\\ &+{\bf
1}_{N^{\prime}}\otimes D_{\hat{v}^{L}}^{-1}\mathcal{D}D_{C}\Big{)}\\\
&+c^{\prime\prime}\otimes{\bf
1}_{N^{\prime}}\otimes\big{(}D_{\hat{v}^{L}}^{-1}(\mathcal{D}-I)\text{diag}(r\odot\hat{i}^{U})W^{\text{iub1}}\big{)}\\\
&+s^{\prime\prime}\otimes{\bf
1}_{N^{\prime}}\otimes\big{(}D_{\hat{v}^{L}}^{-1}(\mathcal{D}-I)\text{diag}(x\odot\hat{i}^{U})W^{\text{iub1}}\big{)}\end{split}$
(24)

and $c^{\prime},s^{\prime}\in\mathbb{R}^{N^{\prime}}$,
$c^{\prime\prime},s^{\prime\prime}\in\mathbb{R}^{N^{\prime\prime}}$ with
entries defined by

$\displaystyle\begin{split}c^{\prime}_{k}&=\frac{\cos((2k-1)\pi/N^{\prime})}{\cos(\pi/N^{\prime})},\quad\forall
k\in\\{1,\ldots,N^{\prime}\\}\\\
s^{\prime}_{k}&=\frac{\sin((2k-1)\pi/N^{\prime})}{\cos(\pi/N^{\prime})},\quad\forall
k\in\\{1,\ldots,N^{\prime}\\}\\\
c^{\prime\prime}_{k}&=\frac{\cos((2k-1)\pi/N^{\prime\prime})}{\cos(\pi/N^{\prime\prime})},\quad\forall
k\in\\{1,\ldots,N^{\prime\prime}\\}\\\
s^{\prime\prime}_{k}&=\frac{\sin((2k-1)\pi/N^{\prime\prime})}{\cos(\pi/N^{\prime\prime})},\quad\forall
k\in\\{1,\ldots,N^{\prime\prime}\\}\end{split}$ (25)

In summary, if $\mathcal{I}\subseteq\mathcal{N}$ denotes the set of bus(es)
with software update(s), then despite all possible worst-case update failure,
if (23) holds then the (squared) line currents are no greater than the upper
limits specified by $\bar{\ell}$.

We note that (15) can also be used to construct line current upper limit
constraints by imposing $I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}}\leq(\bar{\ell})^{1/2}$. This is in fact adopted by our
previous work in [8], requiring $N$ instead of $NN^{\prime}N^{\prime\prime}$
linear constraints in (23) used in this paper. However, (15) is derived by
replacing the current phasors with their magnitudes in (10). This tends to be
conservative especially when there is significant reverse power flows in the
system. Our numerical experiment (not shown in this paper) indeed indicates
that (23) are generally less conservative than (15). In addition, the total
time to construct the constraints in (23) and to solve the corresponding
rollout problem remains modest (see Section V for more detail). This motivates
our proposed current upper limit constraints using (23). Nevertheless, (15)
retains its distinct advantage that it is a upper bound of $(\ell)^{1/2}$
linearly dependent on ${\bf 1}_{\mathcal{I}}$. Usually for
$\mathcal{I}\subset\mathcal{N}$, $I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}}\leq\hat{i}^{U}$ because $\hat{i}^{U}$ assumes update failure
at all buses in $\mathcal{N}$. This observation is utilized to obtain a less
conservative second term in (16) (as in (20)).

###  III-D Voltage Upper Limit Constraint

Note that $M\ell\leq 0$ in (2a) because all entries of $M$ in (3) are non-
positive. Therefore, for any $(p,q)$, $\nu(p,q)\leq\nu_{0}{\bf 1}+2Rp+2Xq$.
Let $\mathcal{I}\subseteq\mathcal{N}$ denote the set of buses with software
updates implying that $(p,q)\in\mathcal{S}(\mathcal{I})$ in (6). Then, for any
$n\in\mathcal{N}$,

$\displaystyle\begin{split}&\max\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}\nu_{n}\\\
\leq&\max\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}\big{\\{}\nu_{0}+2e_{n}^{\top}Rp+2e_{n}^{\top}Xq\big{\\}}\\\
=&\;e_{n}^{\top}\Big{(}\nu_{0}{\bf
1}+2R(\hat{p}^{G}-\hat{p}^{L})+2X(\hat{q}^{G}-\hat{q}^{L})\Big{)}\\\
&+\sum\limits_{i\in\mathcal{I}}\max\limits_{p^{G}_{i}\geq
0,\;q^{G}_{i}}\Big{\\{}2R_{ni}p^{G}_{i}+2X_{ni}q^{G}_{i}\mid|p^{G}_{i}+jq^{G}_{i}|\leq
C_{i}\Big{\\}}\\\
&-\sum\limits_{i\in\mathcal{I}}(2R_{ni}\hat{p}^{G}_{i}+2X_{ni}\hat{q}^{G}_{i})\end{split}$
(26)

where $R_{ni}$ and $X_{ni}$ denote the $(n,i)$ entry of $R$ and $X$,
respectively. Since $R_{ni}\geq 0$, $X_{ni}\geq 0$, it holds in (26) that

$\displaystyle\begin{split}&\max\limits_{p^{G}_{i}\geq
0,\;q^{G}_{i}}\Big{\\{}2R_{ni}p^{G}_{i}+2X_{ni}q^{G}_{i}\mid|p^{G}_{i}+jq^{G}_{i}|\leq
C_{i}\Big{\\}}\\\ =&\;2C_{i}\sqrt{R_{ni}^{2}+X_{ni}^{2}}\\\
:=&\;2C_{i}Z_{ni}\end{split}$ (27)

with

$\displaystyle Z\in\mathbb{R}^{N\times
N},\;\;Z_{ni}:=\sqrt{R_{ni}^{2}+X_{ni}^{2}},\;\;\forall(n,i)$ (28)

Note that in (26) it holds (for the last term) that

$\sum\limits_{i\in\mathcal{I}}(2R_{ni}\hat{p}^{G}_{i}+2X_{ni}\hat{q}^{G}_{i})=e_{n}^{\top}(2RD_{\hat{p}^{G}}+2XD_{\hat{q}^{G}}){\bf
1}_{\mathcal{I}}$ (29)

Thus, (27) and (29) imply that the inequalities in (26) can be assembled for
all $n$ as

$\displaystyle\nu(p,q)\leq\hat{\nu}^{\text{nom}}+W^{\text{vub}}{\bf
1}_{\mathcal{I}},\quad\forall(p,q)\in\mathcal{S}(\mathcal{I})$ (30)

where

$\displaystyle\begin{split}\hat{\nu}^{\text{nom}}&:=\nu_{0}{\bf
1}+2R(\hat{p}^{G}-\hat{p}^{L})+2X(\hat{q}^{G}-\hat{q}^{L})\\\
W^{\text{vub}}&:=2ZD_{C}-(2RD_{\hat{p}^{G}}+2XD_{\hat{q}^{G}})\end{split}$
(31)

In summary, if $\mathcal{I}$ denotes the set of buses with software updates
the maximum squared voltages are upper bounded by
$\hat{\nu}^{\text{nom}}+W^{\text{vub}}{\bf 1}_{\mathcal{I}}$ as specified in
(30). Therefore, if the constraints

$\hat{\nu}^{\text{nom}}+W^{\text{vub}}{\bf 1}_{\mathcal{I}}\leq\bar{\nu}$ (32)

are imposed on ${\bf 1}_{\mathcal{I}}$, then the worst-case possible squared
voltages due to software update failure associated with inverters at buses in
$\mathcal{I}$ cannot exceed $\bar{\nu}$.

We remark that $\hat{\nu}^{\text{nom}}$ in (31) is the vector of nominal
squared voltage (i.e., no software update failure) of the _linearized_
DistFlow equations (i.e., (2a) with the loss term $M\ell$ ignored).

###  III-E Voltage Lower Limit Constraint

Let $\mathcal{I}\subseteq\mathcal{N}$ denotes the set of bus(es) with software
update(s) implying that $(p,q)\in\mathcal{S}(\mathcal{I})$. Then, for any
$n\in\mathcal{N}$, (2a) specifies that

$\displaystyle\begin{split}\nu_{n}&=\nu_{0}+2e_{n}^{\top}Rp+2e_{n}^{\top}Xq+e_{n}^{\top}M\ell\\\
&\geq\nu_{0}+e_{n}^{\top}\Big{(}2Rp+2Xq+M\big{(}\hat{i}^{U}\odot(I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}})\big{)}\Big{)}\end{split}$

where the inequality holds because of the following four facts: (a) all
entries of $M$ are non-positive, (b) $\ell=\ell^{1/2}\odot\ell^{1/2}$, (c) for
all $(p,q)$, it holds that $\ell^{1/2}\leq\hat{i}^{U}$ by (7b), (d) for all
$(p,q)\in\mathcal{S}(\mathcal{I})$, it holds that $\ell^{1/2}\leq
I^{\text{iub1}}+W^{\text{iub1}}{\bf 1}_{\mathcal{I}}$. Therefore, the minimum
$\nu_{n}$ due to $(p,q)\in\mathcal{S}(\mathcal{I})$ satisfies

$\displaystyle\begin{split}&\min\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}\nu_{n}\\\
\geq&\;\nu_{0}+\min\limits_{(p,q)\in\mathcal{S}(\mathcal{I})}\big{\\{}2e_{n}^{\top}Rp+2e_{n}^{\top}Xq\big{\\}}\\\
&+e_{n}^{\top}M\big{(}\hat{i}^{U}\odot(I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}})\big{)}\\\ =&\;e_{n}^{\top}\Big{(}\nu_{0}{\bf
1}_{\mathcal{I}}+2R(\hat{p}^{G}-\hat{p}^{L})+2X(\hat{q}^{G}-\hat{q}^{L})\Big{)}\\\
&+\sum\limits_{i\in\mathcal{I}}\min\limits_{p^{G}_{i}\geq
0,q^{G}_{i}}\Big{\\{}2R_{ni}p^{G}_{i}+2X_{ni}q^{G}_{i}\mid|p^{G}_{i}+jq^{G}_{i}|\leq
C_{i}\Big{\\}}\\\
&-\sum\limits_{i\in\mathcal{I}}(2R_{ni}\hat{p}^{G}_{i}\\!+\\!2X_{ni}\hat{q}^{G}_{i})\\!+\\!e_{n}^{\top}M\big{(}\hat{i}^{U}\\!\odot\\!(I^{\text{iub1}}+W^{\text{iub1}}{\bf
1}_{\mathcal{I}})\big{)}\end{split}$ (33)

Since $R_{ni}\geq 0$, $X_{ni}\geq 0$, in (33) it holds that

$\displaystyle\begin{split}&\min\limits_{p^{G}_{i}\geq
0,q^{G}_{i}}\Big{\\{}2R_{ni}p^{G}_{i}+2X_{ni}q^{G}_{i}\mid|p^{G}_{i}+jq^{G}_{i}|\leq
C_{i}\Big{\\}}\\\ =&-\sum\limits_{i\in\mathcal{I}}2X_{ni}C_{i}\\\
=&-2e_{n}^{\top}XD_{C}{\bf 1}_{\mathcal{I}}\end{split}$ (34)

Thus, (29), (31) and (34) imply that the inequalities in (33) can be assembled
for all $n$ as

$\nu(p,q)\geq\hat{\nu}^{\text{nom}}+M(\hat{i}^{U}\odot
I^{\text{iub1}})-W^{\text{vlb}}{\bf
1}_{\mathcal{I}},\;\;\forall(p,q)\in\mathcal{S}(\mathcal{I})$ (35)

where

$W^{\text{vlb}}:=2XD_{C}+2RD_{\hat{p}^{G}}+2XD_{\hat{q}^{G}}-MD_{\hat{i}^{U}}W^{\text{iub1}}$
(36)

In summary, if $\mathcal{I}$ denotes the set of buses with software updates
then the minimum squared voltages are lower bounded by
$\hat{\nu}^{\text{nom}}+M(\hat{i}^{U}\odot I^{\text{iub1}})-W^{\text{vlb}}{\bf
1}_{\mathcal{I}}$ as in (35). Therefore, if the constraints

$\hat{\nu}^{\text{nom}}+M(\hat{i}^{U}\odot I^{\text{iub1}})-W^{\text{vlb}}{\bf
1}_{\mathcal{I}}\geq\underline{\nu}$ (37)

are imposed on ${\bf 1}_{\mathcal{I}}$, then the worst-case squared voltages
due to software update failure associated with inverters at buses in
$\mathcal{I}$ cannot be lower than $\underline{\nu}$.

###  III-F Software Update Rollout Problem Formulation

The decision variables of the rollout problem are the software update instants
$t_{i}\geq 0$ for bus $i\in\mathcal{N}$. Since the updates should be scheduled
to cause the least interruption to normal operation, the makespan (i.e., the
time when the first update starts to the time when the last update finishes)
should be minimized. Furthermore, the safety constraints should be satisfied.
For given $t_{1},t_{2},\ldots,t_{N}$ and $t$ we define

$\mathcal{I}(t):=\\{i\in\mathcal{N}\mid
t\in\begin{bmatrix}t_{i},t_{i}+\Delta\tau\end{bmatrix}\\}$ (38)

as the set of buses such that time interval
$\begin{bmatrix}t_{i},t_{i}+\Delta\tau\end{bmatrix}$ contains $t$. Recall that
$\Delta\tau$ is the fault clearing time, so $\mathcal{I}(t)$ is the set of
buses whose adverse effect of software update failure can be seen at time $t$
(for given $t_{1},t_{2},\ldots$). Then, to maintain the grid code in (5) in
the worst-case scenario (all scheduled updates fail and the corresponding
nodes suffer the worst possible power injections), constraints (23), (32) and
(37) should hold for $\mathcal{I}(t)$ at all $t$. Hence, the rollout problem
is summarized as

$\displaystyle\begin{split}\underset{t_{i}}{\text{minimize}}&\quad\underset{i\in\mathcal{N}}{\max}\;t_{i}\\\
\text{subject to}&\quad{\bf H}{\bf 1}_{\mathcal{I}(t)}\leq{\bf b},\quad\forall
t\\\ &\quad\text{$\mathcal{I}(t)$ satisfies \eqref{eqn:It}}\end{split}$ (39)

where

$\begin{split}\\!\\!{\bf H}:=\begin{bmatrix}W^{\text{iub2}}\\\
W^{\text{vub}}\\\
W^{\text{vlb}}\end{bmatrix},\;\mathbf{b}:=\begin{bmatrix}{\bf
1}_{N^{\prime\prime}}\otimes({\bf
1}_{N^{\prime}}\otimes(\bar{\ell})^{1/2})-I^{\text{iub2}}\\\
\overline{\nu}-\hat{\nu}^{\text{nom}}\\\
\hat{\nu}^{\text{nom}}+M(\hat{i}^{U}\odot
I^{\text{iub1}})-\underline{\nu}\end{bmatrix}\end{split}$ (40)

Due to the indicator function $\mathcal{I}(t)$, problem (39) is nonlinear and
non-convex. However, as shown in [7, Proposition 1], $t_{i}$ in (39) can in
fact be restricted to nonnegative integer multiples of $\Delta\tau$. Then,
(39) can be reinterpreted as follows. The time axis can be divided into
intervals called time slots of the form $[k\Delta\tau,(k+1)\Delta\tau]$ with
$0\leq k\leq N-1$. The decision to make is the assignment of updates to the
time slots, so that each update is assigned to exactly one time slot. For the
shortest update schedule, the number of time slots used is minimized. Further,
instead of imposing the constraint in (39) for all $t$, it suffices to
consider only integer multiples of $\Delta\tau$. By defining the “discrete-
time” version of $\mathcal{I}$ in (38) with time slot index starting from
$j=1$

$\mathcal{I}_{j}:=\mathcal{I}\big{(}(j-1)\Delta\tau\big{)},\quad
j=1,2,\ldots,N$

problem (39) can be reformulated as follows:

$\displaystyle\begin{split}\text{minimize}&\quad T\\\ \text{subject
to}&\quad\text{$\mathcal{I}_{1},\mathcal{I}_{2},\ldots,\mathcal{I}_{T}$
partition $\mathcal{N}$}\\\ &\quad{\bf H}{\bf 1}_{\mathcal{I}_{j}}\leq{\bf
b},\quad j\in\\{1,\ldots,T\\}\end{split}$ (41)

with ${\bf H}$ and ${\bf b}$ defined in (40). Problem (41) is known as the
_vector bin packing problem_. It can be modeled as an integer program (e.g.,
[7]) or approximately solved using greedy algorithms to be detailed in the
next section.

##  IV Solving Software Update Rollout Problem

To quickly solve large instances of (41), we adopt the best-fit-decreasing
greedy algorithm. First we check if it is possible to let each time slot hold
exactly one update, since the instance is feasible if and only if this is
possible. If feasible, the updates are sorted by their “sizes” specified by
the row vector ${\bf 1}^{\top}{\bf H}$ (large to small). Then, the sorted
updates are considered one-by-one. For each update, we first try to assign it
to the time slot with the minimum “sum-up residual capacity” (to be defined in
Algorithm 2), among the used slots with some update already assigned. Only if
no used slot can hold the update, will a new slot be appended to the schedule.

Algorithm 2 summaries our customized best-fit-decreasing algorithm.

Algorithm 2 Best-fit-decreasing algorithm to solve (41)

1:${\bf H}\in\mathbb{R}^{K\times N}$, ${\bf b}\in\mathbb{R}^{K}$

2:${\bf x}\in\\{0,1\\}^{N\times N}$, ${\bf y}\in\\{0,1\\}^{N}$ if instance is
feasible; ${\bf x}_{ij}=1$ iff update $i$ is assigned to time slot $j$, and
${\bf y}_{j}=1$ iff time slot $j$ is needed in the rollout schedule

3:if ${\bf b}\ngeq{\bf H}(:,i)$ for any $i\in\mathcal{N}$ then
$\triangleright$ feasibility check

4: report infeasibility, QUIT

5:end if

6:Define ${\bf s}={\bf 1}^{\top}{\bf H}$ $\triangleright$ “size” of each
update

7:Define update indices ${\bf u}_{1},\ldots,{\bf u}_{N}:{\bf s}({\bf
u}_{1})\geq{\bf s}({\bf u}_{2})\geq\ldots$

8:Initialize ${\bf x}\leftarrow{\bf 0}_{N\times N}$, ${\bf y}\leftarrow{\bf
0}_{N}$

9:Initialize $N_{y}\leftarrow 0$ $\triangleright$ number of time slots used

10:Initialize ${\bf C}\leftarrow{\bf 0}_{K\times N}$ (residual capacity
vectors for all time slots), ${\bf c}\leftarrow{\bf 1}^{\top}{\bf b}$ (sum-up
residual capacity)

11:Initialize time slot indices ${\bf s}_{1}\leftarrow 1$, …, ${\bf
s}_{N}\leftarrow N$

12:for $i=1,\ldots,N$ do $\triangleright$ loop over sorted updates

13: for $j=1,\ldots,N$ do $\triangleright$ loop over sorted time slots

14: if ${\bf C}(:,{\bf s}_{j})\geq{\bf H}(:,{\bf u}_{i})$ then
$\triangleright$ check used slot

15: $s^{*}\leftarrow{\bf s}_{j}$, quit for-loop over $j$

16: end if

17: if $j>N_{y}$ then $\triangleright$ open new slot

18: $s^{*}\leftarrow j$, $N_{y}+\\!+$, ${\bf C}(:,N_{y})\leftarrow{\bf b}$,
${\bf y}(N_{y})\leftarrow 1$

19: quit for-loop over $j$

20: end if

21: end for

22: Update ${\bf x}({\bf u}_{i},s^{*})\leftarrow 1$

23: Update ${\bf C}(:,s^{*})\mathrel{-}=H(:,{\bf u}_{i})$, ${\bf
c}(s^{*})\mathrel{-}={\bf s}({\bf u}_{i})$

24: while $j\geq 2$ do $\triangleright$ update used time slots ordering

25: if ${\bf c}({\bf s}_{j})<{\bf c}({\bf s}_{j-1})$ then

26: Swap ${\bf s}_{j}$ with ${\bf s}_{j-1}$

27: end if

28: $j-\\!-$

29: end while

30:end for

Assume that the instance is feasible. The assignment ${\bf
C}(:,N_{y})\leftarrow{\bf b}$ in line 16 specifies that whenever a new time
slot is used, its capacity vector is initialized to be ${\bf b}$. Whenever an
update is assigned to the time slot, the modification ${\bf
C}(:,s^{*})\mathrel{-}=H(:,{\bf u}_{i})$ in line 21 ensures that the residual
capacity vector of the slot is properly deducted. Then the condition ${\bf
C}(:,{\bf s}_{j})\geq{\bf H}(:,{\bf u}_{i})$ in line 12 specifies that the
time slot can accept an update only if its residual capacity vector remains
nonnegative after the assignment. Together, the assignment returned by the
algorithm satisfies the second constraints in (41). In addition, the for-loop
over $i$ implies that each update must be assigned to exactly one time slot (a
new slot if necessary). Hence, the first constraint of (41) is satisfied by
the algorithm output. Therefore, the algorithm returns a feasible rollout
whenever the instance is feasible. Though the returned rollout is not
necessarily optimal, it typically requires very few time slots (see Table II
in Section V for detail).

The main loop over $i$ is executed $N$ times. Each iteration of the inner loop
over $j$ requires $O(K)$ operations. The remaining steps in an iteration of
the main loop require $O(N+K)$ operations. Hence, the time-complexity of
Algorithm 2 is $O(N^{2}K)$. In Section V, we demonstrate that Algorithm 2 as
stated can be executed very efficiently even in MATLAB. For example, running
Algorithm 2 for the 10,476-bus instance from [22] requires less than 2
seconds. The overall software rollout scheduling procedure is summarized in
Algorithm 3.

Algorithm 3 Software update rollout scheduling procedure

1:(a) DistFlow equation (1), (b) load $(\hat{p}^{L},\hat{q}^{L})$ and nominal
generation $(\hat{p}^{G},\hat{q}^{G})$, (c) inverter capacity $C$ in (4) and
(d) safety limits $\underline{\nu}$, $\overline{\nu}$ and $\overline{\ell}$ in
(5)

2:Compute universal voltage lower bound $\hat{\nu}^{L}=(\hat{v}^{L})^{2}$ and
current upper bound $\hat{\ell}^{U}=(\hat{i}^{U})^{2}$ according to Algorithm
1

3:Compute $W^{\text{iub2}}$ and $I^{\text{iub2}}$ according to (24), (25) and
(15)

4:Compute $W^{\text{vub}}$ and $\hat{\nu}^{\text{nom}}$ according to (31), (3)
and (28)

5:Compute $W^{\text{vlb}}$, $M$ and $I^{\text{iub1}}$ using (36), (3) and (15)

6:Assemble matrix ${\bf H}$ and vector ${\bf b}$ according to (40)

7:Apply Algorithm 2 to decide if instance $({\bf H},{\bf b})$ is feasible. If
feasible, rollout schedule is defined by output ${\bf x}$

##  V Numerical Demonstration

We demonstrate the practical performance of the proposed software update
rollout scheduling procedure with distribution system benchmarks. All
computation is performed on a Mac Studio (2022) with M1 Ultra CPU and 128 GB
of RAM, running MATLAB through the Rosetta 2 translation environment.

###  V-A Distribution System Benchmarks

The benchmarks considered in this paper include: CIGRE 44-bus low-voltage (LV)
radial network from [23], radial networks from Matpower 7.1 [24] and the REDS
repository [22]. The CIGRE network contains one substation busbar, three
feeder heads and 40 non-reference buses. In our study, the substation and the
feeder heads are merged into one slack bus with squared voltage fixed at
$\nu_{0}=1$ pu. The descriptions (topology, line parameters, load, etc.) of
all other benchmarks can be found in [24] and [22]. We assume that for each
benchmark all non-reference buses are equipped with smart inverters with the
same (benchmark dependent) power rating.

###  V-B Quality of Universal Voltage and Current Bounds

We numerically demonstrate the quality of the voltage lower bound
$\hat{\nu}^{L}$ and current upper bound $\hat{\ell}^{U}$ obtained by Algorithm
1. Due to (7) and Proposition 1, for any $(p,q)$ it holds that
$\nu(p,q)\geq\nu^{L}\geq\hat{\nu}^{L}$ and
$\ell(p,q)\leq\ell^{U}\leq\hat{\ell}^{U}$. Hence, if the difference
$\nu(p,q)-\hat{\nu}^{L}$ is small then the error $\nu^{L}-\hat{\nu}^{L}$ must
be small. A similar relation holds for $\ell(p,q)$, $\ell^{U}$ and
$\hat{\ell}^{U}$. A suitable choice of $(p,q)$ is
$(\tilde{p},\tilde{q})=(-\hat{p}^{L},-\hat{q}^{L}-C)$ because
$(\tilde{p},\tilde{q})$ minimizes $\nu$ in (2a) if the loss term $M\ell$ is
ignored (also making $\ell$ large in in (2b) consequently). For the CIGRE LV
benchmark we compare $\hat{\nu}^{L}$ and $\nu(\tilde{p},\tilde{q})$ (resp.
$\hat{\ell}^{U}$ and $\ell(\tilde{p},\tilde{q})$) in two cases. In both cases
the total inverter power rating is the total active power load. However, for
the first case the default load is considered resulting in typical forward
power flows. On the other hand, for the second case “background” distributed
generation is added so that the total net active power load is zero. This is a
situation with significant distributed generation and reverse power flows.
Fig. 2 indicates that $\hat{\nu}^{L}$ and $\hat{\ell}^{U}$ obtained by
Algorithm 1 are of acceptable quality for the CIGRE benchmark.

Figure 2: Quality demonstration of universal voltage and current bounds with
CIGRE LV network. Solid blue lines are $\hat{\nu}^{L}$ obtained by Algorithm
1. Dotted blue lines are $\nu(\tilde{p},\tilde{q})$, upper bounds of
$\hat{\nu}^{L}$. Solid orange lines are $\hat{\ell}^{U}$ obtained by Algorithm
1. Dotted orange lines are $\ell(\tilde{p},\tilde{q})$, lower bounds of
$\hat{\ell}^{U}$. Left: typical load with forward power flows. The average
relative error in voltage is $0.17\%$ and the relative error in current is
$1.1\%$. Right: zero net active power load with significant reverse power
flows. The average relative error in voltage is $1.4\%$ and the relative error
in current is $15\%$.

The aforementioned experiment is repeated for selected benchmarks from [24]
and [22]. Table I suggests similar conclusions. Except for the approximately
10% average relative error for current upper bound in scenarios with
significant reverse power flows, generally the error is negligible. We note
that, in the experiment above, on average Algorithm 1 converges in 6.5 fixed
point iterations and the maximum number of iterations is 9 over all cases. The
relative tolerance in Algorithm 1 is $10^{-6}$.

TABLE I: Average (over buses or lines) relative error of universal voltage lower bounds and current upper bounds case name | forward power flows | reverse power flows
---|---|---
rel err $\nu$ | rel err $\ell$ | rel err $\nu$ | rel err $\ell$
Matpower 33bw | 0.006% | 0.21% | 0.059% | 5.6%
Matpower 85 | 0.050% | 0.38% | 0.064% | 1.3%
Matpower 118zh | 0.0050% | 0.22% | 0.017% | 11%
Matpower 141 | 0.0018% | 0.047% | 0.023% | 3.9%
REDS 135+8 | $\approx 0$ | 0.022% | 0.0070% | 7.3%
REDS 201+3 | $\approx 0$ | 0.012% | $\approx 0$ | 4.5%
REDS 873+7 | $\approx 0$ | 0.011% | $\approx 0$ | 4.0%
REDS 10476+84 | $\approx 0$ | 0.010% | $\approx 0$ | 4.2%

###  V-C CIGRE 44-bus LV Distribution Grid – Heavy Load

In this case study, we assume that the power rating of each smart inverter is
$0.7/N=0.015$ pu (the system’s total inverter rating is 0.7 pu). The (linear)
voltage limits in (5) are $0.85$ pu and $1.1$ pu for all buses. For each line,
the (linear) current upper limit is the maximum of 600 A and 1.5 times the
nominal current of the line in the default load setting.

To simulate the scenario with heavy load, the nominal distributed generation
is zero (i.e., $\hat{p}^{G}=0$ and $\hat{q}^{G}=0$) and the load is 130% of
the default load (both for $\hat{p}^{L}$ and $\hat{q}^{L}$) so that the total
active power load is 0.89 pu. With all relevant data specified, the rollout
problem in (41) is set up and solved using Algorithm 2. The Gantt diagram of
the obtained update schedule is shown in Fig. 3 (left), utilizing nine time
slots.

Figure 3: Results of CIGRE heavy load scenario case study concerning (41).
Left: Gantt diagram. Right: voltage and current margins.

Fig. 3 (right) also shows the voltage and current margins for each of the time
slots defined as follows: for any time slot with assigned update
$\mathcal{I}\subseteq\mathcal{N}$, the following problems are solved

$\displaystyle\nu^{L}_{n}(\mathcal{I})$
$\displaystyle:=\min\limits_{p,q,P,Q,\nu,\ell}\;\nu_{n}\;\;\text{s.t.~{}\eqref{eqn:DF},
$(p,q)\in\mathcal{S}(\mathcal{I})$ in \eqref{eqn:SI}}$ (42a)
$\displaystyle\nu^{U}_{n}(\mathcal{I})$
$\displaystyle:=\max\limits_{p,q,P,Q,\nu,\ell}\;\nu_{n}\;\;\text{s.t.~{}\eqref{eqn:DF},
$(p,q)\in\mathcal{S}(\mathcal{I})$ in \eqref{eqn:SI}}$ (42b)
$\displaystyle\ell^{U}_{n}(\mathcal{I})$
$\displaystyle:=\max\limits_{p,q,P,Q,\nu,\ell}\;\ell_{n}\;\;\text{s.t.~{}\eqref{eqn:DF},
$(p,q)\in\mathcal{S}(\mathcal{I})$ in \eqref{eqn:SI}}$ (42c)

to obtain the worst-case voltage for each non-reference bus $n\in\mathcal{N}$
and the worst-case current for each line $n\in\mathcal{L}$ for $\mathcal{I}$
associated with the time slot. Then, the worst-case quantities from (42) are
compared with the safety limits in (5) to obtain the relevant margins
(nonnegative values mean limit satisfaction). It can be seen that the rollout
schedule obtained by solving the proposed problem in (41) with Algorithm 2
indeed satisfies all safety limits.

As a comparison, the linearized DistFlow version of (41) is considered. The
following problem (i.e., an extension of the linearized DistFlow rollout
problem in [7]) is solved to obtain the corresponding rollout schedule:

$\displaystyle\begin{split}\text{minimize}&\quad T\\\ \text{subject
to}&\quad\text{$\mathcal{I}_{1},\mathcal{I}_{2},\ldots,\mathcal{I}_{T}$
partition $\mathcal{N}$}\\\ &\quad{\bf\tilde{H}}{\bf
1}_{\mathcal{I}_{j}}\leq{\bf\tilde{b}},\quad j\in\\{1,\ldots,T\\}\end{split}$
(43)

where

$\displaystyle{\bf\tilde{H}}=\begin{bmatrix}W^{\text{vub}}\\\
\tilde{W}^{\text{vlb}}\end{bmatrix},\quad{\bf\tilde{b}}=\begin{bmatrix}\overline{\nu}-\hat{\nu}^{\text{nom}}\\\
\hat{\nu}^{\text{nom}}-\underline{\nu}\end{bmatrix}$

with $W^{\text{vub}}$ and $\hat{\nu}^{\text{nom}}$ defined in (31) and
$\tilde{W}^{\text{vlb}}:=2XD_{C}+2RD_{\hat{p}^{G}}+2XD_{\hat{q}^{G}}$ (i.e.,
the matrix in (36) with the last term removed). Fig. 4

Figure 4: Results of CIGRE heavy load scenario case study concerning (43).
Left: Gantt diagram. Right: voltage and current margins.

shows the Gantt diagram and the margins of the rollout schedule obtained by
solving (43) using Algorithm 2. The linearized DistFlow schedule requires four
versus nine time slots required by the proposed problem in (41). However, Fig.
4 shows negative voltage lower limit margins for the first three time slots,
as well as negative current upper limit margin for the second time slot. This
indicates voltage and current limit violations, as illustrated in Fig. 5.

Figure 5: Voltage and current profiles showing safety limit violations of the
rollout schedule due to (43). Left: voltage profile of time slot 1. Right:
current profile of time slot 2. In both subfigures, ACPF denotes the voltage
and current profiles obtained from AC power flow simulation using Matpower
[24] with some power injections $(p,q)\in\mathcal{S}(\mathcal{I})$ for
$\mathcal{I}$ conforming to the time slot. In the left, LinDistFlow denotes
the minimum voltage due to linearized DistFlow equation (i.e.,
$(\hat{\nu}^{\text{nom}}+W^{\text{vub}}{\bf 1}_{\mathcal{I}})^{1/2}$), which
suggests false limit satisfaction.

This justifies the apparently more complicated safety constraints proposed in
(40) for the rollout scheduling problem, as opposed to the simpler linearized
DistFlow version introduced earlier in [7].

###  V-D CIGRE 44-bus LV Distribution Grid – Significant DG

Opposite to Section V-C, this section considers the scenario with significant
distributed generation and reverse power flows. In this scenario, the loads
$\hat{p}^{L}$ and $\hat{q}^{L}$ retain their default values. The total
inverter power rating is increased to 0.9 pu (versus total active power load
of 0.69 pu). In addition, all inverters nominally generate active power
according to their ratings (i.e., $\hat{p}^{G}=0.9/N{\bf 1}$ pu and
$\hat{q}^{G}={\bf 0}$ pu). The voltage lower limits and current limits and the
current upper limits are all uniform over their elements. They are
respectively 0.9 pu, 1.2 pu and 800 A. Fig. 6

Figure 6: Results of CIGRE significant DG scenario case study concerning (41).
Left: Gantt diagram. Right: voltage and current margins.

shows the Gantt diagram and the voltage and current margins for the rollout
schedule obtained by the proposed method. In comparison, Fig. 7

Figure 7: Results of CIGRE significant DG scenario case study concerning (43).
Left: Gantt diagram. Right: voltage and current margins.

shows the analogous results for the rollout schedule obtained by solving the
linearized DistFlow rollout scheduling problem in (43). The linearized version
requires one fewer time slot. However, voltage lower limits are again
violated.

###  V-E REDS 10,476-bus/83-feeder Distribution Grid

This benchmark is the largest example from the repository in [22], containing
10,476 non-reference buses and 83 feeders. We assume that inverters are
installed in all non-reference buses with a uniform rating. The total rating
is 68.4 pu, while the total active power load is 38.2 pu. The safety limits
are, respectively, 0.9 pu, 1.1 pu and 600 A for voltage and current (uniform
over all relevant elements). Solving the proposed rollout scheduling problem
in (41) with this benchmark yields a schedule with Gantt diagram and voltage
and current margins shown in Fig. 8.

Figure 8: Results of 10,476-bus benchmark case study concerning (41). Left:
Gantt diagram. Right: voltage and current margins.

The schedule requires four time slots to arrange all 10,476 updates, and the
total time to set up and solve the rollout problem is less than 3 seconds
(with Algorithm 2 requiring about 50% of the time). It can be verified that it
is impossible to assign all updates to one single time slot, excluding the
trivial solution. To obtain a rollout schedule with four time slots using
enumeration would require checking the second constraints in (41) more than
$5\times 10^{14}$ times. Suppose each such examination requires $2\times
10^{-5}$ seconds (based on the time recording in our experiment). The total
enumeration would require more than 300 years, which is not suitable for real-
time applications.

###  V-F Summary for all cases

Instances of the rollout problem in (41) derived from benchmarks from the REDS
repository [22] and Matpower 7.1 [24] are solved using Algorithm 2. In all
these cases, the voltage limits are 0.9 pu and 1.1 pu for all buses and the
current upper limits are 600 A or above (depending on the nominal current due
to default load setting). Typical load and inverter rating scenarios are
assumed, which could be different from the scenario in Section V-E. Table II
summarizes the results of the rollout schedules.

TABLE II: Computation time and number of time slots of rollout schedules obtained by the proposed method case name | time | slots | case name | time | slots
---|---|---|---|---|---
REDS 13+3 | 4.98 ms | 1 | 33bw | 4.35 ms | 2
REDS 29+1 | 4.77 ms | 2 | 51he | 6.08 ms | 1
REDS 32+1 | 4.84 ms | 2 | 69 | 8.85 ms | 4
REDS 83+11 | 6.35 ms | 4 | 85 | 11.6 ms | 6
REDS 135+8 | 10.9 ms | 3 | 94pi | 20.5 ms | 11
REDS 201+3 | 30.7 ms | 4 | 118zh | 11.3 ms | 8
REDS 873+7 | 151 ms | 2 | 136ma | 10.2 ms | 3
REDS 10476+84 | 1.97 s | 2 | 141 | 20.1 ms | 5

We note that the computation time for all cases, including the 10476-bus
example, remains modest. In addition, the number of time slots used by the
schedules are not excessive in comparison to the “sequential schedules” where
each time slot holds exactly one update. This suggests, despite all
approximations, that our proposed method is able to obtain safe and effective
rollout schedules in real-time, even for distribution networks of substantial
size.

##  VI CONCLUSIONS

Rapid system-level deployment and reconfiguration of control software is
gaining prominence as power systems evolve to integrate large amounts of
inverter-based distributed resources with little or zero inertia, thus
requiring more dynamic control. As the power system is part of the critical
infrastructure, its control software update should be carefully executed to
prevent major reliability incidents. This requires proper accounting of the
cyber-physical relationship between software behavior (e.g., update failure)
and power system operational state (e.g., voltage and current). For radial
distribution systems, the DistFlow equation accurately describes the
underlying physics. However, the intrinsic nonlinearity renders the equation
unfit for any real-time decision-making. The linearized DistFlow equation is
generally a good balance between model complexity and fidelity. Nevertheless,
as demonstrated in this paper, it may not be accurate enough to maintain
system safety standard when conflicting operational goals are present (e.g.,
makespan minimization for software update rollout). This necessitates
innovations in distribution system modeling exemplified by the proposed
linearized relationships between software update failure and the worst-case
voltage/current, which are theoretically justifiable with our detailed
analysis and practically viable as demonstrated by benchmarks with complexity
far exceeding the norm in the literature (e.g., 10476-bus vs 33-bus). Besides
the software update rollout problem studied here, these linearized
relationships can potentially enable online distribution system contingency
monitoring and re-dispatch. In addition, to enable realistic large-scale
applications computational innovations are needed. The reformulation of the
rollout problem into a bin-packing problem amenable to efficient best-fit
decreasing algorithm is one such example. In addition, the fixed point
iteration procedure to quickly estimate the worst-case voltages and currents
crucially enables the real-time large-scale application of the proposed
rollout scheduling procedure. The investigation to apply the fixed point
procedure to other situations (e.g., quick power flow analysis, contingency
evaluation, etc.) appears to be worthwhile.

## Appendix A – Auxiliary Statements

######  Lemma 1.

Let $a+jb\in\mathbb{C}$ and $r>0$ be given. Also, let

$\displaystyle\theta:=\begin{cases}\tan^{-1}(b/a),&\text{if $a\neq 0$ or
$b\neq 0$}\\\ 0,&\text{if $a=b=0$}\end{cases}$

Then, the maximum objective value of the following problem

$\displaystyle\begin{split}\underset{p,q}{\text{maximize}}&\quad|(a+jb)-(p+jq)|\\\
\text{subject to}&\quad p\geq 0,\;\;|p+jq|\leq r\end{split}$ (44)

is

$\displaystyle\max\big{\\{}$ $\displaystyle|a+j(b+r)|,|a+j(b-r)|,|a+jb-
re^{j\theta}|\big{\\}}$ (45)

###### Proof.

First we consider the special case where $a=b=0$. In this case, the maximum of
(44) is $r$ and all three terms in (45) are equal to $r$. Thus, the statement
is trivially true when $a=b=0$. Next, we consider the case where at least one
of $a$ or $b$ is not zero. Let $(p^{*},q^{*})$ be the argument of maximum of
(44). We claim that $|p^{*}+jq^{*}|=r$. Suppose not, then there exists $\Delta
q\neq 0$ such that $|p^{*}+j(q^{*}+\Delta q)|=r$ and $|b-q^{*}-\Delta
q|>|b-q^{*}|$. This implies that $(p^{*},q^{*}+\Delta q)$ is a better feasible
solution than $(p^{*},q^{*})$, leading to a contradiction. Hence, (44) can be
re-parameterized as

$\displaystyle\underset{\alpha\in[-\pi/2,\;\pi/2]}{\text{maximize}}\quad|r\cos(\alpha)-a+j(r\sin(\alpha)-b)|$
(46)

It can be verified that the only stationary point of the objective for
$-\pi/2\leq\alpha\leq\pi/2$ is at $\tan^{-1}(b/a)$ (well-defined because at
least one of $a$ or $b$ is nonzero). Therefore, the maximum of (46) can be
attained only at one of the three possibilities: $-\pi/2$, $\pi/2$ and
$\tan^{-1}(b/a)$, which is $\theta$ according to the statement definition.
Thus, the maximum of (46) must be one of the following: $|a+j(b+r)|$,
$|a+j(b-r)|$ or $|a+jb-r(\cos(\theta)+j\sin(\theta))|$, same as (45). ∎

######  Lemma 2.

Let $\mathcal{D}$ be the descendant matrix defined in Section II-B. Let
$\hat{p}^{L},\hat{q}^{L},\hat{p}^{G},\hat{q}^{G}\in\mathbb{R}^{N}$ be given.
Then, for any $n\in\mathcal{N}$ and $\mathcal{I}\subseteq\mathcal{N}$, the
optimal objective value of

$\displaystyle\begin{split}\underset{p^{G},q^{G}}{\text{maximize}}\quad&e_{n}^{\top}\big{|}\mathcal{D}(p+jq)\big{|}\\\
\text{subject to}\quad&(p,q)\in\mathcal{S}(\mathcal{I})\;\text{in
\eqref{eqn:SI}}\end{split}$ (47)

is

$\displaystyle\max\Big{\\{}$
$\displaystyle\\!|a\\!+\\!j(b\\!+\\!r)|,|a\\!+\\!j(b\\!-\\!r)|,|a\\!+\\!jb\\!-\\!re^{j\theta}|\\!\Big{\\}}$
(48)

where

$\displaystyle\begin{split}&a:=e_{n}^{\top}\mathcal{D}(\hat{p}^{L}-({\bf
1}-{\bf 1}_{\mathcal{I}})\odot\hat{p}^{G})\\\
&b:=e_{n}^{\top}\mathcal{D}(\hat{q}^{L}-({\bf 1}-{\bf
1}_{\mathcal{I}})\odot\hat{q}^{G})\\\ &r:=e_{n}^{\top}\mathcal{D}({\bf
1}_{\mathcal{I}}\odot C)\\\ &\theta:=\begin{cases}\tan^{-1}(b/a),&\text{if
$a\neq 0$ or $b\neq 0$}\\\ 0,&\text{if $a=b=0$}\end{cases}\end{split}$ (49)

###### Proof.

Note that for any $x\in\mathbb{R}^{N}$,
$e_{n}^{\top}\mathcal{D}x=\sum\limits_{m\in d(n)}x_{m}$. Hence, with
$(p,q)\in\mathcal{S}(\mathcal{I})$ the objective function can be written as

$\displaystyle\left|\sum\limits_{m\in
d(n)}(p^{G}_{m}-\hat{p}^{L}_{m})+j\sum\limits_{m\in
d(n)}(q^{G}_{m}-\hat{q}^{L}_{m})\right|$ $\displaystyle=$
$\displaystyle\left|\big{(}\sum\limits_{m\in
d(n)}\hat{p}^{L}_{m}-\\!\\!\\!\sum\limits_{m\in
d(n)\setminus\mathcal{I}}\hat{p}^{G}_{m}\big{)}\\!+\\!j\big{(}\sum\limits_{m\in
d(n)}\hat{q}^{L}_{m}-\\!\\!\\!\sum\limits_{m\in
d(n)\setminus\mathcal{I}}\hat{q}^{G}_{m}\big{)}\right.$
$\displaystyle\left.-\;\sum\limits_{m\in
d(n)\cap\mathcal{I}}(p^{G}_{m}+jq^{G}_{m})\right|$

According to the definition of $\mathcal{S}(\mathcal{I})$, for each
$m\in\mathcal{I}$, the set of all possible $p^{G}_{m}+jq^{G}_{m}$ is the half-
disk $\Omega(C_{m}):=\\{(x,y)\mid x\geq 0,|x+jy|\leq C_{m}\\}$. Thus, the set
of all possible sum $\sum\limits_{m\in
d(n)\cap\mathcal{I}}(p^{G}_{m}+jq^{G}_{m})$ is the Minkowski sum
$\sum\limits_{m\in d(n)\cap\mathcal{I}}\Omega(C_{m})$, which can be verified
to be the half-disk $\Omega(\sum\limits_{m\in d(n)\cap\mathcal{I}}C_{m})$. By
noting that

$\displaystyle\sum\limits_{m\in
d(n)}\hat{p}^{L}_{m}-\\!\\!\\!\sum\limits_{m\in
d(n)\setminus\mathcal{I}}\hat{p}^{G}_{m}=e_{n}^{\top}\mathcal{D}(\hat{p}^{L}-({\bf
1}-{\bf 1}_{\mathcal{I}})\odot\hat{p}^{G})=a$ $\displaystyle\sum\limits_{m\in
d(n)}\hat{q}^{L}_{m}-\\!\\!\\!\sum\limits_{m\in
d(n)\setminus\mathcal{I}}\hat{q}^{G}_{m}=e_{n}^{\top}\mathcal{D}(\hat{q}^{L}-({\bf
1}-{\bf 1}_{\mathcal{I}})\odot\hat{q}^{G})=b$ $\displaystyle\sum\limits_{m\in
d(n)\cap\mathcal{I}}C_{m}=e_{n}^{\top}\mathcal{D}({\bf 1}_{\mathcal{I}}\odot
C)=r$

problem (47) can be rewritten as maximizing $|a+jb-(x+jy)|$ with respect to
$x$ and $y$, such that $x\geq 0$ and $|x+jy|\leq r$. According to Lemma 1, the
maximum objective value is given in (48). ∎

## Appendix B – Proof of Proposition 1

It suffices to show the following induction: suppose for some $k$ it holds
that $v^{k}\leq(\nu^{L})^{1/2}$ and $i^{k}\geq(\ell^{U})^{1/2}$, then
$v^{k+1}\leq(\nu^{L})^{1/2}$ and $i^{k+1}\geq(\ell^{U})^{1/2}$.

For any $n\in\mathcal{N}$, according to (7) and (2), there exist some
optimizing pair $(p^{G(n)},q^{G(n)})\in\mathbb{R}^{N}\times\mathbb{R}^{N}$
satisfying (4) with $(p^{G},q^{G})$ interpreted as $(p^{G(n)},q^{G(n)})$ such
that $\nu^{L}_{n}=v_{n}^{2}$ with $v$ (together with $i$) satisfying

$\displaystyle\begin{split}v&=\big{(}\nu_{0}{\bf
1}+2R(p^{G(n)}-\hat{p}^{L})+2X(q^{G(n)}-\hat{q}^{L})+Mi^{2}\big{)}^{1/2}\\\
i&=v^{-1}\odot\big{|}\mathcal{D}((p^{G(n)}-\hat{p}^{L})\\!+\\!j(q^{G(n)}-\hat{q}^{L}))\\!-\\!(\mathcal{D}\\!-\\!I)D_{z}i^{2}\big{|}\end{split}$

We note that all entries of $M$ are non-positive. This, together with
$i^{k}\geq(\ell^{U})^{1/2}$ by the induction assumption and
$(\ell^{U})^{1/2}\geq i$ by (7), implies that $Mi^{2}\geq M(i^{k})^{2}$. Thus,
for any $n\in\mathcal{N}$

$\displaystyle(\nu^{L}_{n})^{1/2}$ $\displaystyle=$
$\displaystyle\big{(}\nu_{0}+e_{n}^{\top}\big{(}2R(p^{G(n)}-\hat{p}^{L})+2X(q^{G(n)}-\hat{q}^{L})+Mi^{2}\big{)}\big{)}^{1/2}$
$\displaystyle\geq$
$\displaystyle\big{(}\nu_{0}\\!+\\!e_{n}^{\top}\big{(}2R(p^{G(n)}-\hat{p}^{L})\\!+\\!2X(q^{G(n)}-\hat{q}^{L})\\!+\\!M(i^{k})^{2}\big{)}\big{)}^{1/2}$
$\displaystyle\geq$ $\displaystyle\min_{\text{$(p^{G},q^{G})$ satisfying
\eqref{eqn:inverter_constraints}}}\big{\\{}\big{(}\nu_{0}+e_{n}^{\top}\big{(}2R(p^{G}-\hat{p}^{L})$
$\displaystyle+2X(q^{G}-\hat{q}^{L})+M(i^{k})^{2}\big{)}\big{)}^{1/2}\big{\\}}$
$\displaystyle=$ $\displaystyle e_{n}^{\top}\big{(}\nu_{0}{\bf
1}\\!-\\!2R\hat{p}^{L}-2X(\hat{q}^{L}\\!+\\!C)\\!+\\!M(i^{k})^{2}\big{)}^{1/2}$
$\displaystyle=$ $\displaystyle v_{n}^{k+1}$

Similarly, there exists some
$(p^{g(n)},q^{g(n)})\in\mathbb{R}^{N}\times\mathbb{R}^{N}$ satisfying (4) with
$(p^{G},q^{G})$ interpreted as $(p^{g(n)},q^{g(n)})$ such that
$\ell^{U}_{n}=i_{n}^{2}$ with $i$ (together with $v$) satisfying

$\displaystyle\begin{split}v&=\big{(}\nu_{0}{\bf
1}+2R(p^{g(n)}-\hat{p}^{L})+2X(q^{g(n)}-\hat{q}^{L})+Mi^{2}\big{)}^{1/2}\\\
i&=v^{-1}\odot\big{|}\mathcal{D}((p^{g(n)}-\hat{p}^{L})\\!+\\!j(q^{g(n)}-\hat{q}^{L}))\\!-\\!(\mathcal{D}\\!-\\!I)D_{z}i^{2}\big{|}\end{split}$

Since all entries of $(\mathcal{D}-I)D_{r}$ and $(\mathcal{D}-I)D_{x}$ are
nonnegative, $i^{k}\geq(\ell^{U})^{1/2}$ and $v^{k}\leq(\nu^{L})^{1/2}$ by the
induction assumption, $(\ell^{U})^{1/2}\geq i$ and $(\nu^{L})^{1/2}\leq v$ by
(7), it holds that (a) $v\geq v^{k}$, (b)
$(\mathcal{D}-I)D_{r}i^{2}\leq(\mathcal{D}-I)D_{r}(i^{k})^{2}$ and (c)
$(\mathcal{D}-I)D_{x}i^{2}\leq(\mathcal{D}-I)D_{x}(i^{k})^{2}$. Hence, for any
$n\in\mathcal{N}$

$\displaystyle(\ell_{n}^{U})^{1/2}$ $\displaystyle=$ $\displaystyle
v_{n}^{-1}e_{n}^{\top}\big{|}\mathcal{D}((p^{g(n)}-\hat{p}^{L})\\!+\\!j(q^{g(n)}-\hat{q}^{L}))\\!-\\!(\mathcal{D}\\!-\\!I)D_{z}i^{2}\big{|}$
$\displaystyle\leq$
$\displaystyle{(v_{n}^{k})}^{-1}\\!\\!e_{n}^{\top}\big{(}\big{|}\mathcal{D}((p^{g(n)}\\!\\!-\\!\hat{p}^{L})\\!\\!+\\!\\!j(q^{g(n)}\\!\\!-\\!\hat{q}^{L}))\big{|}\\!\\!+\\!\\!\big{|}(\mathcal{D}\\!\\!-\\!\\!I)D_{z}(i^{k})^{2}\big{|}\big{)}$
$\displaystyle\leq$ $\displaystyle\max_{\text{$(p^{g},q^{g})$ satisfying
\eqref{eqn:inverter_constraints}}}{(v_{n}^{k})}^{-1}\big{|}e_{n}^{\top}\mathcal{D}((p^{g}-\hat{p}^{L})+j(q^{g}-\hat{q}^{L}))\big{|}$
$\displaystyle+{(v_{n}^{k})}^{-1}e_{n}^{\top}\big{|}(\mathcal{D}-I)D_{z}(i^{k})^{2}\big{|}$

According to Lemma 2 with $\mathcal{I}=\mathcal{N}$, the first term above can
be rewritten as
$\bar{S}_{n}:=\max\\{|a_{n}+j(b_{n}+r_{n})|,|a_{n}+j(b_{n}-r_{n})|,|a_{n}+jb_{n}-r_{n}(\cos(\theta_{n})+j\sin(\theta_{n})|\\}$,
where $a_{n}=e_{n}^{\top}\mathcal{D}\hat{p}^{L}$,
$b_{n}=e_{n}^{\top}\mathcal{D}\hat{q}^{L}$, $r_{n}=e_{n}^{\top}\mathcal{D}C$
and $\theta_{n}=\tan^{-1}(b_{n}/a_{n})$ if at least one of $a_{n}$ and $b_{n}$
is nonzero (otherwise $\theta_{n}=0$). Therefore,

$\displaystyle(\ell_{n}^{U})^{1/2}\leq{(v_{n}^{k})}^{-1}\left(\bar{S}_{n}+e_{n}^{\top}\big{|}(\mathcal{D}-I)D_{z}(i^{k})^{2}\big{|}\right)=i_{n}^{k+1}\;\;\;\qed$

## References

  * [1] K. E. Hemsley and R. E. Fisher, “History of industrial control system cyber incidents,” Idaho National Lab. (INL), Idaho Falls, ID (United States), Tech. Rep. INL/CON-18-44411-Rev002, Dec 2018.
  * [2] A. A. Cárdenas, “Cyber-physical systems security,” CyBOK, Knowledge Area Report, 2019.
  * [3] A. A. Cardenas, S. Amin, and S. Sastry, “Research challenges for the security of control systems,” in _3rd USENIX Workshop on Hot Topics in Security (HotSec 08)_. San Jose, CA: USENIX Association, Jul. 2008. [Online]. Available: https://www.usenix.org/conference/hotsec-08/research-challenges-security-control-systems
  * [4] P. Griffioen, S. Weerakkody, and B. Sinopoli, “A moving target defense for securing cyber-physical systems,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 5, pp. 2016–2031, 2021.
  * [5] A. Kanellopoulos and K. G. Vamvoudakis, “A moving target defense control framework for cyber-physical systems,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 3, pp. 1029–1043, 2020.
  * [6] R. Romagnoli, B. H. Krogh, D. de Niz, A. D. Hristozov, and B. Sinopoli, “Software rejuvenation for safe operation of cyber–physical systems in the presence of run-time cyberattacks,” _IEEE Transactions on Control Systems Technology_ , vol. 31, no. 4, pp. 1565–1580, 2023.
  * [7] M. G. de Medeiros, K. C. Sou, and H. Sandberg, “Minimum-time secure rollout of software updates for controllable power loads,” _Electric Power Systems Research_ , vol. 189, p. 106797, 2020.
  * [8] K. C. Sou and H. Sandberg, “Resilient scheduling of control software updates in power distribution systems,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 6146–6153.
  * [9] A. Froger, M. Gendreau, J. E. Mendoza, Éric Pinson, and L.-M. Rousseau, “Maintenance scheduling in the electricity industry: A literature review,” _European Journal of Operational Research_ , vol. 251, no. 3, 2016.
  * [10] M. E. Baran and F. F. Wu, “Optimal capacitor placement on radial distribution systems,” _IEEE Transactions on Power Delivery_ , vol. 4, no. 1, pp. 725–734, Jan 1989.
  * [11] ——, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Transactions on Power Delivery_ , vol. 4, no. 2, pp. 1401–1407, April 1989.
  * [12] C. Coffrin, H. L. Hijazi, and P. Van Hentenryck, “Strengthening the sdp relaxation of ac power flows with convex envelopes, bound tightening, and valid inequalities,” _IEEE Transactions on Power Systems_ , vol. 32, no. 5, pp. 3549–3558, 2017.
  * [13] D. K. Molzahn and L. A. Roald, “Towards an ac optimal power flow algorithm with robust feasibility guarantees,” in _2018 Power Systems Computation Conference (PSCC)_ , 2018, pp. 1–7.
  * [14] R. Louca and E. Bitar, “Robust ac optimal power flow,” _IEEE Transactions on Power Systems_ , vol. 34, no. 3, pp. 1669–1681, 2019.
  * [15] X. Wu, A. J. Conejo, and N. Amjady, “Robust security constrained ACOPF via conic programming: Identifying the worst contingencies,” _IEEE Transactions on Power Systems_ , vol. 33, no. 6, 2018.
  * [16] C. Coffrin, R. Bent, B. Tasseff, K. Sundar, and S. Backhaus, “Relaxations of ac maximal load delivery for severe contingency analysis,” _IEEE Transactions on Power Systems_ , vol. 34, no. 2, 2018.
  * [17] J. Liu, B. Cui, D. K. Molzahn, C. Chen, X. Lu, and F. Qiu, “Optimal power flow in dc networks with robust feasibility and stability guarantees,” _IEEE Transactions on Control of Network Systems_ , vol. 9, 2021.
  * [18] “The grid edge revolution: Innovative drivers towards net-zero energy,” 2019, Siemens white paper. [Online]. Available: https://new.siemens.com/global/en/company/topic-areas/smart-infrastructure/grid-edge/white-paper-grid-edge-net-zero-energy-drivers.html
  * [19] F. Kintzler, T. Gawron-Deutsch, S. Cejka, J. Schulte, M. Uslar, E. M. Veith, E. Piatkowska, P. Smith, F. Kupzog, H. Sandberg, M. S. Chong, D. Umsonst, and M. Mittelsdorf, “Large scale rollout of smart grid services,” in _2018 Global Internet of Things Summit (GIoTS)_ , June 2018.
  * [20] D. Shelar and S. Amin, “Security assessment of electricity distribution networks under der node compromises,” _IEEE Transactions on Control of Network Systems_ , vol. 4, no. 1, pp. 23–36, 2016.
  * [21] D. Shelar, S. Amin, and I. A. Hiskens, “Evaluating resilience of electricity distribution networks via a modification of generalized benders decomposition method,” _IEEE Transactions on Control of Network Systems_ , vol. 8, no. 3, pp. 1225–1238, 2021.
  * [22] R. Kavasseri and C. Ababei. Reds: Repository of distribution systems. [Online]. Available: http://www.dejazzer.com/reds.html
  * [23] L. Thurner, A. Scheidler, F. Schäfer, J. Menke, J. Dollichon, F. Meier, S. Meinecke, and M. Braun, “pandapower — an open-source python tool for convenient modeling, analysis, and optimization of electric power systems,” _IEEE Transactions on Power Systems_ , vol. 33, Nov 2018.
  * [24] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power systems research and education,” _IEEE Transactions on power systems_ , vol. 26, no. 1, pp. 12–19, 2010.

Kin Cheong Sou received a Ph.D. degree in Electrical Engineering and Computer
Science at Massachusetts Institute of Technology in 2008. From 2008 to 2010 he
was a postdoctoral researcher at Lund University, Lund, Sweden. From 2010 to
2012 he was a postdoctoral researcher at KTH Royal Institute of Technology,
Stockholm, Sweden. Between 2013 and 2016 he was an assistant professor with
the department of Mathematical Sciences, Chalmers University of Technology and
the University of Gothenburg, Sweden. Dr. Sou is now an associate professor
with the department of Electrical Engineering at the National Sun Yat-sen
University in Taiwan. His research interests include optimization and system
analysis with power systems applications.
---
Henrik Sandberg is Professor at the Division of Decision and Control Systems,
KTH Royal Institute of Technology, Stockholm, Sweden. He received the M.Sc.
degree in engineering physics and the Ph.D. degree in automatic control from
Lund University, Lund, Sweden, in 1999 and 2004, respectively. From 2005 to
2007, he was a Post-Doctoral Scholar at the California Institute of
Technology, Pasadena, USA. In 2013, he was a Visiting Scholar at the
Laboratory for Information and Decision Systems (LIDS) at MIT, Cambridge, USA.
He has also held visiting appointments at the Australian National University
and the University of Melbourne, Australia. His current research interests
include security of cyber-physical systems, power systems, model reduction,
and fundamental limitations in control. Dr. Sandberg was a recipient of the
Best Student Paper Award from the IEEE Conference on Decision and Control in
2004, an Ingvar Carlsson Award from the Swedish Foundation for Strategic
Research in 2007, and a Consolidator Grant from the Swedish Research Council
in 2016. He has served on the editorial boards of IEEE Transactions on
Automatic Control and the IFAC Journal Automatica. He is Fellow of the IEEE.
---