TSO

    transmission system operator
AC

    alternating current
HVDC

    high voltage direct current
API

    application programming interface
SCOPF

    security constrained optimal power flow
OPF

    optimal power flow
ED

    economic dispatch
IMML

    inverse matrix multiplication lemma
LP

    linear program
OPS

    optimal probabilistic security
SCOPF

    security-constrained optimal power flow
P-SCOPF

    preventive SCOPF
C-SCOPF

    corrective SCOPF
RMAC

    reliability management approach and criterion
IMML

    inverse matrix modification lemma
SM

    Sherman-Morrison
SMW

    Sherman-Morrison-Woodbury
DRC

    deterministic reliability criterion
MINLP

    mixed-integer nonlinear program
TSO

    transmission system operator
PTDF

    power transfer distribution factors
LODF

    line outage distribution factors
CENS

    cost of energy not served
ENS

    energy not served
LOL

    loss of load
FOR

    forced outage rate
EOC

    expected operational cost
FCR

    frequency containment reserves
aFRR

    automatic frequency restoration reserves
mFRR

    manual frequency restoration reserves
FFR

    fast frequency reserves
VOLL

    value of lost load
VRE

    variable renewable energy
RTS

    reliability test system

# A Fast and Scalable Iterative Solution of a Socio-Economic Security-
Constrained Optimal Power Flow with Two-Stage Post-Contingency Control
††thanks: The research leading to results in this paper has received funding
from the Research Council of Norway through the project “Resilient and
Probabilistic reliability management of the transmission grid” (RaPid) (Grant
No. 294754), The Norwegian Water Resources and Energy Directorate, and
Statnett.

Matias Vistnes1, Vijay Venu Vadlamudi1, and Oddbjørn Gjerde2 Trondheim, Norway
1 Department of Electric Energy, Norwegian University of Science and
Technology  2 Department of Energy Systems, SINTEF Energy Research

###### Abstract

A fast and scalable iterative methodology for solving the security-constrained
optimal power flow (SCOPF) problem is proposed using problem decomposition and
the inverse matrix modification lemma. The SCOPF formulation tackles system
security operational planning by using short- and long-term post-contingency
limits, probability of branch outages, and preventive and corrective actions,
a probabilistic corrective-SCOPF problem formulation. Using two post-
contingency states and contingency probabilities, the SCOPF could provide good
system security at a lower cost than only preventive actions as the typical
‘N-1’-formulation does. Additional security is ensured using a post-
contingency load-shedding limit constraint based on system operator policy.
The proposed methodology is applied to a range of test systems containing up
to 10,000 buses with a computational time of up to 3375 s for all 12,706
branch contingencies. Calculating contingency power flows takes 1.3% of the
total solution time using the proposed methodology exploiting the inverse
matrix modification lemma.

###### Index Terms:

Contingency Analysis, Optimal Power Flow, Probabilistic Methods, Inverse
Matrix Modification Lemma, System Security Planning.

##  I Introduction

Reliability studies have computational bottlenecks in the fast-growing power
systems where the system’s state changes faster due to higher penetration of
renewable energy, new loads, and more frequent extreme weather. The system is
operated closer to the equipment’s limits due to increased power needs from
the energy transition and electrification, while power system investments are
not keeping up with demand growth. Consequently, it is vital to have accurate
solutions to ensure the reliability of electricity supply for customers at an
affordable price; society is more dependent on it than ever before.

In a traditional reliability management approach, reliability can be
accomplished through grid development solutions that consistently meet a
deterministic reliability criterion. Operators often use preventive
rescheduling of the economic dispatch to meet the criterion. However, such
solutions may be prohibitively expensive and not socio-economically efficient
when not including corrective actions (generator redispatch, load shedding,
and other transmission system operator (TSO) actions after a contingency).
More probabilistic and risk-based methods, where reliability and costs are
balanced, have been proposed in literature; an example is from the EU-funded
project GARPUR: “taking into account probabilities of contingencies and using
a compound socio-economic objective function blending the costs of TSO
preventive and corrective actions with a monetization of the risk of service
interruptions” [1]. It is shown with a socio-economic cost assessment that
reducing the power flow on critical lines increase power market costs, while
interruption costs decrease [2].

###  I-A Related work

SCOPF is an extensively used tool in reliability studies, whose output could
be used to obtain quantitative reliability metrics [1, 3]. The seminal works
on security-constrained optimal power flow (SCOPF) were first published by
Alsac and Stott in [4] using preventive actions (P-SCOPF) and by Monticelli,
Pereira and Granville in [5] using both preventive and corrective actions
(C-SCOPF), and since then variants of SCOPF has found several applications in
the planning and operation of power systems. However, accounting for
contingencies $\mathcal{C}$ increases the problem size dramatically. Each
contingency leads to a different post-contingency state, which ideally is
modeled with the same detail as the pre-contingency system. That results in a
problem size $|\mathcal{C}|$ times bigger than without accounting for
contingencies. Furthermore, if corrective post-contingency actions are
available, the number of decision variables also multiplies by
$|\mathcal{C}|$. This is a computational challenge, especially with a large
system size. In addition, risk and chance constraints could be added to
account for uncertainties in demand, generation, or outages, which also
increase the complexity of the problem [6, 7].

Corrective actions in probabilistic methods need new types of statistical data
compared to deterministic methods. These data points are contingency
probabilities and the failure rates of the different types of corrective
actions. Increasingly good data is available for the former, though it is hard
to come by any data for the latter. Corrective action failures could be
important to model [8], though the present paper will not discuss this.
However, the present paper will study the impact of contingency probabilities.
Moreover, how the post-contingency system is modeled impacts the outcome. A
post-contingency state could be set immediately, minutes, or hours after the
contingency, all with different operating limits and possible control actions.
At least two post-contingency states are proposed to maintain good security
[9].

Recent research shows the need for better methods for reliability studies to
deal with the computational challenges, as a direct approach of SCOPF of large
power systems is computationally intractable [1]. There is a lot of variation
in formulations of SCOPF within the literature, with many using an iterative
solution procedure. Ensuring post-contingency feasibility while modeling the
primary response of generators in a DC preventive SCOPF (P-SCOPF) is shown for
a large power system [10]; the primary response is here modeled with binary
variables. Corrective actions and outages could be modeled as nodal injections
for use in corrective SCOPF (C-SCOPF) [11]. A probabilistic C-SCOPF is
proposed in [12], referred to as an optimal probabilistic security (OPS)
approach, which uses an iterative procedure with a particle swarm solution
method on a small system. In that approach, weather conditions were
incorporated. Analyzing a large power system is possible using a risk-based DC
C-SCOPF with Lagrangian relaxation and Benders decomposition [6]. Several
works identify the necessary and sufficient conditions to describe the
feasible set of solutions for a DC P-SCOPF using pre-processing of the
problem, significantly reducing the size of the problem at the expense of the
pre-processing time [13, 14]. To cope with the computational challenges,
machine learning has been employed in reliability management research to
predict the system state and generation schedule [15]. After a longer initial
training of the machine learning model, such methods use a small fraction of
the time to find a solution compared to the original SCOPF optimization. In
the last few years, they have seen great improvement in accuracy and variation
in their use within SCOPF formulations. Machine learning will not be further
discussed in the present paper. There is a lack of scalable methods of a
probabilistic C-SCOPF type which would aid operators in system security
operational planning optimizing for socio-economic cost benefit.

###  I-B Contributions and paper organization

Combining excellent previous work on mathematical methods could help balance
computational tractability and accuracy of SCOPF and form the motivation for
the propositions made in the present paper. 1) The inverse matrix modification
lemma (IMML)111Also called the Sherman-Morrison formula, but has been
published separately by many authors, so we use the name IMML. efficiently
computes the inverse of a rank-one update of a matrix using the inverse of the
original matrix [16]. The present paper uses the IMML to solve the contingency
power flow problem [17]. 2) The sparse matrix solver KLU proposed in [18] is
used as an alternative calculation method for the contingency power flow. KLU
is a linear system solver based on LU-factorization and specifically made to
exploit the structure of power system matrices. 3) The SCOPF problem is
decomposed into a main problem to which cuts are added, and sub-problems are
analyzed using the IMML or KLU solver. The seminal work on SCOPF used this
approach to divide the problem based on the contingency cases [5].

The present paper contributes towards reliability studies using a socio-
economic C-SCOPF, which has stringent requirements for solution speed. Based
on the output from the proposed methodology, it is possible to obtain
quantitative reliability metrics; the subsequent reliability assessment is not
in the scope of the present paper. An efficient and scalable methodology to
identify the optimal socio-economic operation of a power system with
preventive and corrective actions is proposed. The cutting-plane method is
used to iteratively add control actions to the main formulation. The linear
cuts are calculated using the IMML and the sparse matrix solver KLU.
Preventive actions are applied to the pre-contingency (base case) generation
schedule, and corrective actions are planned and accounted for, even when the
system separates into islands. The proposal extends a previously developed
methodology for DC C-SCOPF [19]. The DC approximation is chosen as it is
linear and the scope of the study is limited to managing the thermal limits of
power system components. The C-SCOPF model could be a part of the daily
operational security planning of a TSO. Before running the C-SCOPF, a unit
commitment and the latest state estimation of the power system need to be
available. In practical applications, further post-processing using AC power
flow is recommended to ensure voltage stability and a feasible operation
state.

The unique contributions of the present paper are:

  1. 1.

an improved algorithm based on the previous paper [19] by the authors, making
a) the method’s optimal value of the objective function equal to that of the
direct approach, and b) the solution time faster,

  2. 2.

a comparison of the methods for contingency power flow calculations used in
the proposed methodology, and

  3. 3.

a detailed case study discussing the two post-contingency states proposed for
more accurate steady-state simulations, where combinations of post-contingency
states are compared to their expected operational cost (EOC), cost of energy
not served (CENS), and run time.

Next, we discuss the post-contingency modeling decisions in Section II. Then,
Section III presents the C-SCOPF formulation, contingency power flow solution
methods, and the proposed methodology. Further, Section IV presents and
discusses the case study. Finally, Section V presents a concluding summary of
the work.

##  II Sequential post-contingency states and actions

CENS serves an active role in probabilistic methods. Its costs and power
system impact are weighted against that of generator rescheduling and other
possible corrective actions. The flexible load in power systems is increasing,
but other loads have grown more reliant on continuous power. There are also
soft constraints on what the public tolerates of load shedding before the
operator’s reputation is impacted. The Norwegian TSO (Statnett) currently has
an operating policy that restricts load shedding after a single credible
contingency to a maximum of 500 MW (restrictions apply) [20]. This is about 2%
of the maximum load in the area of responsibility. A robust constraint is set
up based on this to restrict high load-shedding percentages after a
contingency. In the case study, the limit is set to 2%.

Another aspect of post-contingency actions is the temporal evolution. There is
an inherent danger in working with static models in that dynamic phenomena can
occur between steady-state solutions, leading to unstable trajectories,
tripping of relays, or other control mechanisms that cascade the failures.
This issue has been tackled in different ways: a line loading risk factor [7],
a line failure rate constraint [21], a robust convex restriction on the load
uncertainty set [22], and two post-contingency states [9]. Focusing on the
last, which provides safer control strategies than with one post-contingency
state by better modeling the equipment operating limits using the C-SCOPF [9].
The model used in the present paper adopts this. Operators can take different
actions in the two post-contingency states: the first post-contingency state
is set after circuit breakers have tripped, and minutes later, after the
reserves are active, the second post-contingency state is set.

The post-contingency model is a simplification of the balancing reserve market
as it is organized in the Nordic countries with a short- and a long-term post-
contingency state. The Nordic reserve market is divided into production
rescheduling (mitigating problems such as overloaded branches and over- or
under-voltage) and frequency reserves (after contingencies such as an outage
of production plants or system split). Rescheduling is mainly done through
special regulation, where the operator approves a production bid outside the
merit order. This can be used both as a preventive action and a corrective
action. The fast frequency reserves (FFR) market is not modeled in the steady-
state model used in the present paper. Frequency deviations in the short-term
post-contingency state are managed by the frequency containment reserves (FCR)
and the automatic frequency restoration reserves (aFRR) markets. These markets
balance the system within its emergency limits, usually valid for 15 minutes
of operation. The long-term post-contingency state model results from the
manual frequency restoration reserves (mFRR) market. Compared to markets where
the reserves are automatically activated, it has manual activation and keeps
the system frequency at its nominal value while allowing long-term emergency
ratings (usually valid for 4 hours). These markets, and thus post-contingency
states, have different goals; the fast products require fast ramping but
minimal energy while the manual frequency restoration reserves have less
stringent ramping but more energy needs to fulfill. Moreover, the markets are
cleared separately, and a power plant cannot necessarily provide products for
all markets. Hence, the pre-contingency state is a better basis for the market
in each post-contingency state than using the previous post-contingency state.

##  III Method

The methodology described in the present paper is a continuation of the work
in [19]. Details of the changes made are provided in this section, while a
short summary follows. In the C-SCOPF model, a maximum load-shedding
constraint is added. Calculation of the post-contingency power transfer
distribution factors (PTDF) matrix is done one vector at a time. The iteration
procedure is changed and a pre-interation procedure is added.

The power system consists of sets of buses $\mathcal{N}$, branches
$\mathcal{B}$, generators $\mathcal{G}$, and demands $\mathcal{D}$, where
$N=|\mathcal{N}|$ and $B=|\mathcal{B}|$. Contingencies are denoted $c$ and
taken from the contingency set $\mathcal{K}$. Matrix and vector variables are
in bold.

###  III-A Corrective Security-Constrained Optimal Power Flow

Extending the formulation in [9], we use the short-term limits, load shedding
and decrease in active power production as control variables, do not only
restrict the long-term operational limits to the pre-contingency limits, and
add a maximum load shedding constraint.

$\displaystyle\min_{\bm{x}^{0},\bm{u}^{0},\bm{x}^{k},\bm{u}^{k}}\quad$
$\displaystyle\sum_{g\in\mathcal{G}}c_{g}p_{g}+\sum_{k\in\mathcal{K}}\sum_{i\in\bm{u}^{k}}\pi^{k}u^{k}_{i}$
(1a) $\displaystyle s.t.\quad$
$\displaystyle\bm{g^{0}}(\bm{x}^{0},\bm{u}^{0})\leq\bm{\overline{h}}^{0}$ (1b)
$\displaystyle\bm{g^{k}}(\bm{x}^{k},\bm{u}^{k1})\leq\bm{\overline{h}}^{k1}$
$\displaystyle\quad k\in\mathcal{K}$ (1c)
$\displaystyle|\bm{u}^{k1}-\bm{u}^{0}|\leq\Delta\bm{u}^{k1}$
$\displaystyle\quad k\in\mathcal{K}$ (1d)
$\displaystyle\bm{g^{k}}(\bm{x}^{k},\bm{u}^{k2})\leq\bm{\overline{h}}^{k2}$
$\displaystyle\quad k\in\mathcal{K}$ (1e)
$\displaystyle|\bm{u}^{k2}-\bm{u}^{0}|\leq\Delta\bm{u}^{k2}$
$\displaystyle\quad k\in\mathcal{K}$ (1f)

where $\bm{x}$ is the state variables, $\bm{u}$ is the control variables,
$\bm{x}^{k}$, and $\bm{u}^{k}$, $c_{g}$ is the marginal cost of generator $g$,
$\pi^{k}$ is the probability of contingency $k$, $\bm{g}$ is abstract notation
for the operating limits (including the power flow equations and TSO
constraints), and $\Delta\bm{u}^{k}$ is the maximum allowed change of control
variables (active power generation and active power decrease or increase
between the pre- and post-contingency). Superscripts $0$ and $k$ represent the
base and post-contingency states, respectively. $\bm{\overline{h}}^{0}$,
$\bm{\overline{h}}^{k1}$, and $\bm{\overline{h}}^{k2}$ are operating limits
for normal operation, short-term limits after a contingency, and long-term
limits after a contingency, respectively.

The modeling of the operating limits $\bm{g}$ in (1) is given in the
following. The power flow equations are modeled using PTDFs $\bm{\varphi}$
with the appropriate variables for the specific state. In addition, power
balance is needed in all states.

$\displaystyle-\bm{\overline{h}}_{l}\leq\bm{\varphi}\cdot\left[\bm{u}_{g}-\left(\bm{\overline{h}}_{d}-\bm{u_{d}}\right)\right]\leq\bm{\overline{h}}_{l}$
(2) $\displaystyle\sum_{g\in\mathcal{G}}u_{g}-\sum_{d\in\mathcal{D}}u_{d}=0$
(3)

The post-contingency generation control and load shedding is limited by

$\displaystyle 0\leq u_{g}^{0}+\Delta u_{g}^{k}\leq\overline{h}_{g}^{0}$
$\displaystyle\qquad k\in\mathcal{K},\quad g\in\mathcal{G}$ (4) $\displaystyle
0\leq u_{d}^{0}+u_{d}^{k}\leq\overline{h}_{d}^{0}$ $\displaystyle\qquad
k\in\mathcal{K},\quad d\in\mathcal{D}$ (5)

ensuring that active power production is within the generator limits. In the
pre-contingency state, $\Delta u_{g}^{k}$ and $u_{d}^{k}$ are removed from (4)
and (5).

Further, the maximum allowed post-contingency load shedding $\Gamma$ is
modeled as a TSO constraint (from the Norwegian TSO policy [20]). We use this
method of load shed constraint as a means to ensure that no outages have
higher consequences than are accepted by the operator. It is a form of risk-
averse control, a risk policy, for the operator as the linear value of lost
load (VOLL) is not able to capture all costs in extreme circumstances.

$\sum_{d\in\mathcal{D}}u_{d}\leq\Gamma$ (6)

In summary, the model minimizes the EOC, including preventive and corrective
actions, such as pre-contingency active power production, active power
production re-dispatch, and load-shedding using the probability for each
contingency in the objective function [1]. The generator cost is a piecewise
linear representation of the generator cost curve, where each piece is modeled
as a separate generating unit [23]. The optimality problem is constrained with
all N-1 branch contingencies (transmission lines and transformers). The
operating limits are maximum active power in branches (short- and long-term),
power balance, generator maximum active power output, generator ramping limit
(up and down), and maximum load shedding described by (1b) – (1f) with details
in (2) – (6). Suitable operating limits are selected in the pre- and post-
contingency states.

A DC power flow neglects reactive power and voltage magnitudes to make the
formulation linear. The recent development of convex relaxations of AC power
flow could be used to include reactive power and voltage magnitudes.
Alternatively, active power losses could be included with a compensation
scheme. The present paper does not include power losses, as other parts of the
model and solution methodology are highlighted. However, a compensation scheme
will be needed in an operational setting. This could easily be included.

###  III-B Contingency screening

Modeling the pre- and post-contingency power flow uses the DC power flow
approximation. The general modeling of the power flow is first given before a
more scalable solution method is shown using the IMML.

####  III-B1 DC power flow

The DC power flow equation for node injected active power $\bm{P}$ is

$\displaystyle\bm{P}=\bm{\Phi}^{\mathsf{T}}\cdot\bm{\Psi}\cdot\bm{\Phi}\cdot\bm{\theta}=\bm{H}\cdot\bm{\theta}$
(7) $\displaystyle\bm{\theta}=\bm{H}\setminus\bm{P}=\bm{X}\cdot\bm{P}$ (8)

where $\bm{\Phi}$ is the connectivity matrix, $\bm{\Psi}$ is the diagonal
susceptance matrix, $\bm{\theta}$ is the node voltage angles, $\bm{H}$ is the
susceptance matrix, and $\bm{X}$ its inverse. The power flow $\bm{F}$ can be
found by

$\displaystyle\bm{F}=\bm{\Psi}\cdot\bm{\Phi}\cdot\bm{\theta}$ (9)
$\displaystyle\bm{F}=\bm{\Psi}\cdot\bm{\Phi}\cdot\bm{X}\cdot\bm{P}=\bm{\varphi}\cdot\bm{P}$
(10)

where $\bm{\varphi}$ is the PTDF matrix.

For each contingency, the power flow $\bm{F}^{k}$ is calculated using the
latest iteration of the generation schedule. For a contingency denoted by
superscript $k$:

$\displaystyle\bm{F}^{k}=\bm{\Psi}^{k}\cdot\bm{\Phi}^{k}\cdot\bm{\theta}^{k}$
(11)
$\displaystyle\bm{F}^{k}=\bm{\varphi}^{k}\cdot(\bm{P}+\Delta\bm{P}^{k})=\bm{\varphi}^{k}\cdot\bm{P}^{k}$
(12)

where $\bm{\varphi}^{k}$ is the contingency PTDF matrix (effectively the line
outage distribution factors (LODF)), $\bm{\Psi}^{k}$ and $\bm{\Phi}^{k}$ are
simple two- and one-element-per-branch contingency changes to $\bm{\Psi}$ and
$\bm{\Phi}$, and $\bm{P}^{k}$ is received from the generation schedule.
$\Delta\bm{P}^{k}$ is the manifestation of the corrective actions. It can
often be zero, but is essential to be included when not zero.

Both the IMML and KLU solver can be used to calculate $\bm{\theta}^{k}$ and
$\bm{\varphi}^{k}$. The main differences are highlighted in Table I. Using a
diakoptic solution method, the IMML could also be used with system splitting
when more than one branch is outaged. However, that solution is not
necessarily faster than using a sparse matrix solver [17].

TABLE I: DC contingency power flow methods. An ´X’ in $\Delta\bm{P}$ (or System Separation) means the method can always be used when the power changes (or the system separates into two or more subsystems) after a contingency.  Method | Output | Based on | $\Delta\bm{P}$ | System Separation
---|---|---|---|---
I | $\bm{\theta}^{k}$ | IMML | * | **
II | $\bm{\theta}^{k}$ | KLU | X | X
III | $\bm{\varphi}^{k}$ | IMML | X | **
IV | $\bm{\varphi}^{k}$ | KLU | X | X
*With pre-processing.
**Only implemented for single outage contingencies.

Method II employs the KLU solver to calculate the contingency power flow using
$\bm{H}^{k}$, which is easily found from $\bm{H}$ by reversing the impact of
the branch contingency on two diagonal and two off-diagonal elements.

$\displaystyle\bm{\theta}^{k}=\bm{H}^{k}\setminus\bm{P}^{k}$ (13)

The power flow is then obtained using (11).

Similarly, Method IV employs KLU using $\bm{H}^{k}$ to calculate
$\bm{\varphi}^{k}$ using $\bm{\Psi}^{k}$ and $\bm{\Phi}^{k}$.

$\displaystyle\bm{\varphi}^{k}[m,:]=\bm{H}^{k}\setminus(\bm{\Phi}^{k}[m,:]\cdot\bm{\Psi}^{k})$
(14)

$\bm{\varphi}[i,:]$ and $\bm{\varphi}[:,j]$ are notations for the $i$th row
and $j$th column of matrix $\bm{\varphi}$, respectively. The resulting vector
$\bm{\varphi}^{k}[i,:]$ models the power flow on branch $m$ after the branch
outage described by $k$. To calculate $\bm{\varphi}^{k}$ (14) is repeated for
all branches or $\bm{\phi}^{k}[m,:]$ is replaced by $\bm{\phi}^{k}$ in (14).
Note that $\bm{\varphi}^{k}$ could later be used to calculate the power flow
$\bm{F}^{k}$ using (12), but $\bm{\theta}^{k}$ cannot be used to find
$\bm{\varphi}^{k}$.

####  III-B2 Inverse matrix modification lemma

The contingency screening process does calculate the post-contingency power
flow repeatedly. Thus, efficient performance and scalability are needed for
the calculations. Methods I and III adopt the IMML to calculate the
contingency power flow using the pre-contingency system matrices to achieve
the same values as Methods II and IV, respectively. Method I and III
(initially) assume a connected system after the contingency.

Derivation of the equations are found in [24, 19]; Method I calculate
$\bm{F}^{k}$ of branch $l$ from node $i$ to node $j$ using IMML, the following
are used

$\displaystyle\bm{\theta}^{k}=\bm{\theta}^{0}-\frac{\bm{\delta}\cdot\left(\theta^{0}[i]-\theta^{0}[j]\right)}{1/H^{0}[i,j]+\delta[i]-\delta[j]}$
(15)
$\displaystyle\bm{\delta}=\frac{x_{l}}{H^{0}[i,j]}\cdot\left(\bm{X}^{0}[:,i]-\bm{X}^{0}[:,j]\right)$
(16)

where $x_{l}$ is the reactance of branch $l$. If $\bm{P}^{k}$ is different
than $\bm{P}^{0}$, an intermediate $\bm{\theta}^{\prime}$ must be calculated
using

$\displaystyle\bm{\theta}^{\prime}=\bm{X}^{0}\cdot\bm{P}^{k}$ (17)

which replaces $\bm{\theta}^{0}$ in (15).

Method III adopts the IMML to calculate $\bm{\varphi}^{k}$ using pre-
contingency system matrices and $\bm{\Psi}^{k}$ and $\bm{\Phi}^{k}$. The
intermediate product $\bm{X}^{k}$, the inverse contingency susceptance matrix,
could be found directly by inverting the contingency susceptance matrix,
$\bm{X}^{k}=(\bm{H}+\bm{\Delta H})^{-1}$. However, using the IMML to find
$\bm{X}^{k}$ is more efficient using:

$\displaystyle\bm{X}^{k}=\bm{X}^{0}-\frac{\bm{X}^{0}\cdot\bm{\Phi}^{0}[i,:]\cdot
H^{0}[i,j]\cdot\bm{\delta}}{1+H^{0}[i,j]\cdot\bm{\delta}\cdot\bm{\Phi}^{0}[i,:]}$
(18)
$\displaystyle\bm{\varphi}^{k}=\bm{\Psi}^{k}\cdot\bm{\Phi}^{k}\cdot\bm{X}^{k}$
(19)

The computation can be further streamlined using the knowledge of
$\bm{\Phi}^{k}$ and $\bm{\Psi}^{k}$. The connectivity matrix $\bm{\Phi}$ has
two values for each branch, 1 on node $i$ and -1 on node $j$. $\bm{\Psi}$ only
has values on the diagonal. Using that knowledge and simple arithmetics (18) –
(19) can be written as (superscript 0 is dropped for simplicity):

$\displaystyle\bm{\varphi}^{k}[m,:]=x_{m}\cdot(\bm{X}[:,n]-\bm{X}[:,o]-\alpha\cdot(\bm{X}[:,i]-\bm{X}[:,j]))$
(20)
$\displaystyle\alpha=\frac{(X[i,o]-X[i,n])-(X[j,o]-X[j,n])}{1/x_{l}+(X[i,j]-X[i,i])-(X[j,j]-X[j,i])}$
(21)

for an overload on branch $m$ from node $n$ to node $o$. Eq. (20) is more
efficient because $\alpha$ is a single number and $\bm{\varphi}^{k}[m,:]$ is
calculated using only a few vectors.

Methods I and III are slightly changed when a single contingency cuts off
nodes radially connected to the rest of the system post-contingency. The cut-
off nodes are disconnected from the ground in $\bm{H}^{0}$ by zeroing out the
corresponding rows and columns before applying the described methods. For
higher-order contingencies, a diakoptic solution method could be adopted [17].

###  III-C Algorithm

The solution algorithm, an improvement over the previous version [19], is
shown in Fig. 1. The optimization problem (1) is decomposed into a main
problem, consisting of the pre-contingency state, and sub-problems, each
consisting of a post-contingency state. The total number of system states is
$2C+1$. A constraint (cut) is added to the main problem for each post-
contingency overloaded branch. Please refer to [19] for the deduction of
equations for the cuts.

Figure 1: Flow chart of the proposed algorithm for solving a C-SCOPF.

A particularly important algorithm step is in handling the power flow
constraints. This is changed from [19] and exploits the fast contingency power
flow computation. At the program’s start, no power flow constraints are added,
and the model is similar to an economic dispatch problem, only considering
generation dispatch, and disregarding the system connectivity. After the
initial generation dispatch is found, branch constraints are added to the pre-
contingency state for overloaded branches. Consequently, the optimization
problem is much smaller and thus easier and faster to solve than using all
branch constraints.

A pre-iteration procedure is added to the algorithm that serves a twofold
purpose. It adds corrective actions for every contingency that results in
system separation. The island without the reference node shed all load.
Moreover, the pre-iteration procedure does a contingency screening by
calculating the contingency power flow using Method I (or Method II for
comparison). It ranks the contingencies from maximum branch overload to the
lowest (or no overload). During this pre-iteration procedure the SCOPF is not
re-solved; the system is only re-solved after all variables for corrective
actions have been added to the model, as seen in the middle of Fig. 1. As a
double benefit, the ranking of contingencies results in cuts being added for
the most ‘difficult’ contingencies first; the large changes to the generation
schedule were forced earlier in the main iteration procedure, resolving other
contingencies with a similar but lower impact on branch loading. Both sorting
based on short- and long-term branch limits were tried. Sorting based on
short-term limits was deemed better as these limits more often change the pre-
contingency schedule to the end optimal schedule in our testing. The model is
re-solved right before the main iteration procedure to find the new pre-
contingency generation schedule. In addition, the re-solving balances the
remaining sub-system using post-contingency corrective rescheduling for each
contingency where the system separates.

In the second half of the algorithm (the bottom of Fig. 1) the necessary cuts
are added to the model to enforce the branch flow limits $g$ in (1). First,
(1b) is re-evaluated. Then, (1c) and (1e) are ensured through the iteration
procedure. Intermediate voltage angles are calculated if there are corrective
actions for the post-contingency state, before the contingency power flow is
calculated. Next, suppose there are overloaded branches in the system after a
contingency with the current generation schedule. In that case,
$\bm{\varphi}^{k}$ is calculated by Method III and is used in the addition of
a cut for the short- or long-term branch limits.

As in the pre-iteration procedure, the contingency power flow is always
calculated by Method I, a change from the algorithm presented in [19]
resulting from the implementation of system separation of single outages in
Method I. Moreover, the full dense matrices ($X$ and $\varphi$) are not needed
in the algorithm; only the parts of these matrices corresponding to the
overloaded branches are needed for each contingency. Lastly, the standard
linear solver in Julia was switched to the KLU solver to improve speed.

##  IV Case Study

The efficiency of the proposed methodology is demonstrated by reporting
results from a range of systems, from 24 to 10,000 buses [25, 26], with
accompanying discussions. The systems are set to peak load according to the
system description. VOLL for load shedding are found in [27]. The DC
approximation models the systems using all buses and branches (lines and
transformers). The probability of a contingency $\pi^{k}$ is set as the
average yearly outage probability for an hour for each branch. All N-1 branch
contingencies are considered in each SCOPF for each system. For the 10,000 bus
system, considering all 12,706 branch contingencies using the formulated
C-SCOPF, there are 85 million decision variables and 373 million constraints.
In other terms, this formulation is intractable for direct solution. However,
using the proposed algorithm, the problem starts out with 6655 decision
variables and 16,655 constraints, iteratively growing as needed but never
approaching the size of the naive problem formulation.

Simulations are run for three variants of SCOPF:

  * •

SCOPF: Only a pre-contingency state is optimized.

  * •

P-SCOPF: A pre-contingency and one post-contingency state are optimized. Post-
contingency overloads are based on the long-term limits by using preventive
control actions.

  * •

C-SCOPF: A pre-contingency and two post-contingency states (both short and
long-term limits) with preventive and corrective control actions are
optimized.

All results are gathered on an Intel® CoreTM i9-9920X, 128 GB RAM server. The
methodology is implemented in the Julia programming language v1.9 [28] using
JuMP v1.22.2 [29] and PowerSystems v3.3.0 [30], linear solver KLU v0.6.0 [18],
and problem optimization using the Gurobi v10.0 solver [31]. The reported
solution times include the model preparation time, intermediate calculations,
and LP solver time.

The case study is divided into two parts. First, the computational time of the
proposed solution algorithm is shown for a breadth of different system sizes.
Further, feasibility and the improved computational time for large power
systems enable analysis of the systems using C-SCOPF. A sensitivity analysis
is presented in the second part of the case study.

###  IV-A Computational time

####  IV-A1 DC power flow

Fig. 2 shows the computational time for DC contingency power flow. Method IV
is the slowest method for all system sizes. The difference between methods
increases with the number of nodes in the system. For the 10,000 bus system,
Method III is three times faster than Method IV. Further, Method II is two
orders of magnitude faster, and Method I is over a magnitude faster again. In
total, Method I using IMML is four orders of magnitude faster than Method IV
based on KLU and $\bm{\varphi}^{k}$.

Figure 2: Computational time of DC contingency power flow with change of post-
contingency power $\Delta P$ and system separation.

$\bm{X}$ and $\bm{\varphi}$ are typically dense matrices for all power
systems. They are derived from the $\bm{B}$, $\bm{\Phi}$, and $\bm{\Psi}$
matrices that typically for large power systems are very sparse. Moreover, the
computational power and memory needed to calculate and store the dense
matrices can be considerable. Therefore, if only part of a matrix is needed,
that is the only part that should be calculated. Only one row at a time from
$\bm{\varphi}^{k}$ is needed when iteratively adding cuts. Using Method III
and Method IV for calculation of the $\bm{\varphi}^{k}$ sub-matrix (a vector)
leads to similar computational times as for Method I and Method II calculating
the power flow, respectively.

####  IV-A2 SCOPF

The proposed algorithm is used on system sizes up to 10,000 buses to show the
algorithm’s scalability. Table II shows the computational time using different
system sizes and SCOPF variants. The computational time (time to solution) for
the 10,000 bus system with no contingencies is $0.42\text{\,}\mathrm{s}$,
while considering preventive and corrective security actions for contingencies
gives a computational time of $3375\text{\,}\mathrm{s}$.

TABLE II: Computational time for solving the specified SCOPF. |  |  | Computational time [s]
---|---|---|---
System | Nodes | Branches | SCOPF | P-SCOPF | C-SCOPF
IEEE RTS | 24 | 38 |  $<$0.01 |  $<$0.01 |  $<$0.01
ACTIVSg500 | 500 | 597 | 0.01 | 0.11 | 1.33
ACTIVSg2000 | 2000 | 3206 | 0.05 | 1.73 | 18.85
ACTIVSg10k | 10,000 | 12,706 | 0.42 | 168 | 3375

A breakdown of the computational time into different categories is shown in
Fig 3. As also seen in Table II, setting up and solving the problem using the
SCOPF formulation (initialize pre-contingency in Fig. 3) takes negligible time
for the larger systems. This shows that a focus on post-contingency modeling
is warranted. Moreover, using better approximations for the pre-contingency
could be considered while keeping the DC post-contingency approximation.
Calculation of contingency power flow, a focus area of the present paper, only
uses 0.9–3.2% of the total computational time. Hence, the contingency
screening is not a limiting factor in the algorithm. 19.1–29.1% of the time is
used inside the optimization solver after new cuts are added to the model.

The pre-contingency procedure takes up to 60.6% of the computational time for
the larger systems. This time is mostly spent adding variables and constraints
for the two post-contingency states in the model for the system splitting
contingencies. Note that removing the pre-iteration procedure would not
eliminate the computation time; it just spreads out during the iteration
procedure. Moreover, implementation, choice of optimization solver, solver
interface, and programming language influence this block of time greatly.
Without these variables one would need another scheme to ensure enough ramping
capacity in the correct locations to balance the power. An alternative is to
use a predefined distributed slack that ensures reserve capacity on the
participating generators for all foreseeable contingencies. This would also
greatly reduce the model’s size such that the time spent on re-solving the
model would also be reduced.

Figure 3: Breakdown of the time used for each part of the solution algorithm.

###  IV-B Sensitivity analysis

The sensitivity of EOC from the method C-SCOPF is analyzed with respect to the
input parameters using the synthetic 500-bus SouthCarolina500 system [26]. The
system has 597 branch contingencies, of which 254 branch contingencies split
the system by isolating one bus.

The sensitivity of the input parameters to the output values is found by
random sampling. 1000 random samples were drawn for each parameter (see Table
III) as input to a C-SCOPF of the 500-bus system. The result is summarized in
Fig. 4 which shows the correlation between the input parameters, VOLL, $\pi$,
$\Gamma$, $u^{k}$, $b^{k1}$, and $b^{k2}$, and the output costs, base costs
(pre-contingency costs), short-term corrective load shedding, long-term
corrective load shedding, and EOC. Note that EOC is a combination of the other
output costs, and it is energy not served (ENS) and not CENS in the heatmap.
The input parameters were drawn from uniform distributions and did not show
any correlation in Fig 4. The base costs and short-term ENS show a strong
negative correlation. While short-term ENS also negatively correlates with the
other parameters, it is not as strong. In addition, the base costs and the EOC
positively correlate with $\pi$ and VOLL, with EOC showing the strongest
correlation. Further investigating this relationship, Fig. 5 shows the two
parameters $\pi$ and VOLL multiplied. It is evident from looking at Fig. 5
that both $\pi$ and VOLL need to have low values for the C-SCOPF to schedule
the cheapest pre-contingency generation.

TABLE III: Parameter multipliers with its uniform distribution. Parameter | Range
---|---
VOLL | 0.1–10
$\pi$ | 0.1–10
$\Gamma$ | 0.1–50
$u^{k}$ | 1–10
$\Delta u^{k2}$ | 10
$b^{k1}$ | 1.25–1.5
$b^{k2}$ | 1.0–1.25
Figure 4: Heatmap for the 500-bus system. Figure 5: Scatter plot of base costs
vs VOLL multiplied with $\pi$ for the 500-bus system.

Maximum allowed load shed per contingency $\Gamma$ has a slight negative
correlation with base costs, as shown in Fig 4. The scatter plot in Fig. 6
shows the $\Gamma$ and base cost correlation is limited; for low values of
$\Gamma$ there is a strong negative correlation. However, for higher values
the base cost is not affected by $\Gamma$.

Figure 6: Scatter plot of base costs vs $\Gamma$ for the 500-bus system.

The sensitivity analysis suggests that the main drivers of the input
parameters for low base cost are low $\pi$ and VOLL, while high base cost is
driven by the restriction of post-contingency load shedding ($\Gamma$).
Increasing the generator post-contingency rescheduling cost has less impact on
the base costs in this case study of the 500-bus system.

In the C-SCOPF objective the generation cost is balanced against CENS, both
pre- and post-contingency, where CENS is VOLL times ENS. Thus, VOLL is as
essential as the marginal cost of generators in a SCOPF. However, the VOLL is
much harder to quantify than the marginal generation cost as the ENS could be
from many different customers whose perceptions of VOLL differ. Removing all
system splitting contingencies, Fig. 7 shows how changing the VOLL changes the
pre-contingency generation cost for the SCOPF variants. For up to 30% VOLL the
C-SCOPF makes the same schedule as the SCOPF (without any contingencies).
Further, by increasing the VOLL, the base cost increases as more preventive
measures are scheduled to lessen the CENS post-contingency. Above 3000% VOLL,
the base cost of the C-SCOPF is the same as that from the P-SCOPF.

Figure 7: Base case cost for 500-bus system changing the VOLL. Acronyms are
for the pre-contingency state (P), and for the post-contingency states they
are the short- (C1) and long-term (C2).

With respect to the discussion in this case study, it is essential to note
that these results are gathered on specific test-power systems, and the
conclusions are not necessarily general. However, the trends seem similar for
the systems in this case study.

##  V Conclusion

A fast and scalable iterative solution methodology of a socio-economic C-SCOPF
with two post-contingency states is proposed. Moreover, essential to good
scalability in the methodology are the contingency power flow and optimization
cuts, which are calculated using the IMML. Results from the case study using
the proposed methodology did show the feasibility of a 10,000-bus system
C-SCOPF on modest computer hardware, an Intel® CoreTM i9-9920X with 128 GB
RAM. Still, the method depends on the number of contingencies constraining the
system in the optimum, including contingencies that separate the system.
Optimizing using only preventive actions for the 10,000-bus system took
$168\text{\,}\mathrm{s}$ and also including two-post-contingency states with
corrective actions resulted in a computational time of
$3375\text{\,}\mathrm{s}$. There is still room for further improvements,
especially regarding the pre-iteration procedure where the post-contingency
power balance is ensured. Alternatively, a modeling option is to use
distributed slack with pre-contingency reserved capacity for ramping. The
sensitivity analysis shows cost-correlations with $\pi$, VOLL, and $\Gamma$.
$\pi$ are system state estimations and could be erroneous. VOLL and $\Gamma$
are set parameter quantities based on customer values and system operator
policy, respectively, but could be hard to quantify. These parameters should
be further investigated for other systems using C-SCOPF.

A natural extension of the proposed method includes multi-contingency outages
and islanded operation of the system. Using the results from the proposed
C-SCOPF and a C-SCOPF with only one post-contingency state, a dynamic time-
domain simulation could be conducted to show how the dynamic security changes
with one post-contingency state, compared to two post-contingency states.
Integrating methods for contingency screening and concurrent contingency
calculation will further increase the applications for which the proposed
methodology can be used.

## References

  * [1] GARPUR Consortium  “Generally accepted reliability principle with uncertainty modelling and through probabilistic risk assessment”, 2017  URL: https://cordis.europa.eu/docs/results/608/608540/final1-d11-1d-garpur-final-report.pdf
  * [2] Iver Bakken Sperstad et al.  “Probabilistic socio-economic cost assessment integrating power market and reliability analysis”  In _2018 IEEE International Conference on Probabilistic Methods Applied to Power Systems (PMAPS)_ , 2018, pp. 1–6  DOI: 10.1109/PMAPS.2018.8440342
  * [3] Florin Capitanescu  “Critical review of recent advances and further developments needed in AC optimal power flow”  In _Electric Power Systems Research_ 136, 2016, pp. 57–68  DOI: https://doi.org/10.1016/j.epsr.2016.02.008
  * [4] O. Alsac and B. Stott  “Optimal Load Flow with Steady-State Security”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-93.3, 1974, pp. 745–751  DOI: 10.1109/TPAS.1974.293972
  * [5] A. Monticelli, M… Pereira and S. Granville  “Security-Constrained Optimal Power Flow with Post-Contingency Corrective Rescheduling”  In _IEEE Transactions on Power Systems_ 2.1, 1987, pp. 175–180  DOI: 10.1109/TPWRS.1987.4335095
  * [6] Qin Wang, James D. McCalley, Tongxin Zheng and Eugene Litvinov  “Solving corrective risk-based security-constrained optimal power flow with Lagrangian relaxation and Benders decomposition”  In _International Journal of Electrical Power & Energy Systems_ 75 Elsevier, 2016, pp. 255–264  DOI: 10.1016/J.IJEPES.2015.09.001
  * [7] Line Roald, Maria Vrakopoulou, Frauke Oldewurtel and Göran Andersson  “Risk-based optimal power flow with probabilistic guarantees” The Special Issue for 18th Power Systems Computation Conference.  In _International Journal of Electrical Power & Energy Systems_ 72, 2015, pp. 66–74  DOI: https://doi.org/10.1016/j.ijepes.2015.02.012
  * [8] Efthymios Karangelos and Louis Wehenkel  “Post-contingency corrective control failure: A risk to neglect or a risk to control?”  In _2018 International Conference on Probabilistic Methods Applied to Power Systems, PMAPS 2018 - Proceedings_ Institute of ElectricalElectronics Engineers Inc., 2018  DOI: 10.1109/PMAPS.2018.8440348
  * [9] Florin Capitanescu and Louis Wehenkel  “Improving the Statement of the Corrective Security-Constrained Optimal Power-Flow Problem”  In _IEEE Transactions on Power Systems_ 22.2, 2007, pp. 887–889  DOI: 10.1109/TPWRS.2007.894850
  * [10] Alexandre Velloso, Pascal Van Hentenryck and Emma S. Johnson  “An exact and scalable problem decomposition for security-constrained optimal power flow”  In _Electric Power Systems Research_ 195, 2021, pp. 106677  DOI: https://doi.org/10.1016/j.epsr.2020.106677
  * [11] Pedro Jose Martínez-Lacañina, José Luis Martínez-Ramos, Antonio Villa-Jaén and Alejandro Marano-Marcolini  “DC corrective optimal power flow based on generator and branch outages modelled as fictitious nodal injections”  In _IET Generation, Transmission & Distribution_ 8.3, 2014, pp. 401–409  DOI: https://doi.org/10.1049/iet-gtd.2013.0322
  * [12] J. He, L. Cheng, D.S. Kirschen and Y. Sun  “Optimising the balance between security and economy on a probabilistic basis”  In _IET Generation, Transmission & Distribution_ 4, 2010, pp. 1275–1287(12)
  * [13] Ali Jahanbani Ardakani and François Bouffard  “Identification of Umbrella Constraints in DC-Based Security-Constrained Optimal Power Flow”  In _IEEE Transactions on Power Systems_ 28.4, 2013, pp. 3924–3934  DOI: 10.1109/TPWRS.2013.2271980
  * [14] Richard Weinhold and Robert Mieth  “Fast Security-Constrained Optimal Power Flow Through Low-Impact and Redundancy Screening”  In _IEEE Transactions on Power Systems_ 35.6, 2020, pp. 4574–4584  DOI: 10.1109/TPWRS.2020.2994764
  * [15] Yize Chen, Ling Zhang and Baosen Zhang  “Learning to solve DCOPF: A duality approach”  In _Electric Power Systems Research_ 213, 2022, pp. 108595  DOI: https://doi.org/10.1016/j.epsr.2022.108595
  * [16] Jack Sherman and Winifred J. Morrison  “Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix”  In _The Annals of Mathematical Statistics_ 21.1 Institute of Mathematical Statistics, 1950, pp. 124–127  DOI: 10.1214/aoms/1177729893
  * [17] O. Alsac, B. Stott and W.F. Tinney  “Sparsity-Oriented Compensation Methods for Modified Network Solutions”  In _IEEE Transactions on Power Apparatus and Systems_ 102, 1983, pp. 1050–1060
  * [18] Timothy A. Davis and Ekanathan Palamadai Natarajan  “Algorithm 907: KLU, A Direct Sparse Solver for Circuit Simulation Problems”  In _ACM Transactions on Mathematical Software_ 37.3 New York, NY, USA: Association for Computing Machinery, 2010  DOI: 10.1145/1824801.1824814
  * [19] Matias Vistnes, Vijay Venu Vadlamudi, Sigurd Hofsmo Jakobsen and Oddbjørn Gjerde  “Solving Security-Constrained Optimal Power Flow Using Analytical Benders’ Cuts”  In _2023 International Conference on Smart Energy Systems and Technologies (SEST)_ , 2023, pp. 1–6  DOI: 10.1109/SEST57387.2023.10257488
  * [20] Statnett  “Statnetts driftspolicy – gjennomgang og presisering”, 2019  URL: https://www.statnett.no/globalassets/om-statnett/samarbeidsforum/2019-12-13----driftspolicy---orientering-til-samarbeidsforum.pdf
  * [21] Anirudh Subramanyam, Jacob Roth, Albert Lam and Mihai Anitescu  “Failure Probability Constrained AC Optimal Power Flow”  In _IEEE Transactions on Power Systems_ 37.6, 2022, pp. 4683–4695  DOI: 10.1109/TPWRS.2022.3146377
  * [22] Dongchan Lee, Konstantin Turitsyn, Daniel K. Molzahn and Line A. Roald  “Robust AC Optimal Power Flow With Robust Convex Restriction”  In _IEEE Transactions on Power Systems_ 36.6, 2021, pp. 4953–4966  DOI: 10.1109/TPWRS.2021.3075925
  * [23] Brian Stott and Eric Hobson  “Power System Security Control Calculations Using Linear Programming, Part II”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-97.5, 1978, pp. 1721–1731  DOI: 10.1109/TPAS.1978.354665
  * [24] Max A. Woodbury  “Inverting modified matrices”  Princeton, N. J.: Princeton University, 1950, pp. 4
  * [25] Probability Methods Subcommittee  “IEEE Reliability Test System”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-98.6, 1979, pp. 2047–2054  DOI: 10.1109/TPAS.1979.319398
  * [26] Adam B. Birchfield et al.  “Grid Structural Characteristics as Validation Criteria for Synthetic Networks”  In _IEEE Transactions on Power Systems_ 32.4, 2017, pp. 3258–3265  DOI: 10.1109/TPWRS.2016.2616385
  * [27] M. Fotuhi-Firuzabad and R. Billinton  “Impact of load management on composite system reliability evaluation short-term operating benefits”  In _IEEE Transactions on Power Systems_ 15.2, 2000, pp. 858–864  DOI: 10.1109/59.867185
  * [28] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B. Shah  “Julia: A Fresh Approach to Numerical Computing”  In _SIAM Review_ 59.1, 2017, pp. 65–98  DOI: 10.1137/141000671
  * [29] Miles Lubin et al.  “JuMP 1.0: Recent improvements to a modeling language for mathematical optimization” In press.  In _Mathematical Programming Computation_ , 2023
  * [30] José Daniel Lara et al.  “PowerSystems.jl — A power system data management package for large scale modeling”  In _SoftwareX_ 15, 2021, pp. 100747  DOI: https://doi.org/10.1016/j.softx.2021.100747
  * [31] Gurobi Optimization, LLC  “Gurobi Optimizer Reference Manual”, 2023  URL: https://www.gurobi.com

  *[TSO]: transmission system operator
  *[SCOPF]: security-constrained optimal power flow
  *[P-SCOPF]: preventive SCOPF
  *[C-SCOPF]: corrective SCOPF
  *[OPS]: optimal probabilistic security
  *[IMML]: inverse matrix modification lemma
  *[EOC]: expected operational cost
  *[CENS]: cost of energy not served
  *[FFR]: fast frequency reserves
  *[FCR]: frequency containment reserves
  *[aFRR]: automatic frequency restoration reserves
  *[mFRR]: manual frequency restoration reserves
  *[PTDF]: power transfer distribution factors
  *[PTDFs]: power transfer distribution factors
  *[VOLL]: value of lost load
  *[LODF]: line outage distribution factors
  *[ENS]: energy not served


# Third-order topological insulator induced by disorder

Hugo Lóio1,2, Miguel Gonçalves1, Pedro Ribeiro1,3, Eduardo V. Castro1,3,4
1CeFEMA, LaPMET, Instituto Superior Técnico, Universidade de Lisboa, Av.
Rovisco Pais, 1049-001 Lisboa, Portugal  2Laboratoire de Physique Théorique et
Modélisation, CNRS UMR 8089, CY Cergy Paris Université, 95302 Cergy-Pontoise
Cedex, France  3Beijing Computational Science Research Center, Beijing 100084,
China  4Centro de Física das Universidades do Minho e Porto, LaPMET,
Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do
Porto, 4169-007 Porto, Portugal

###### Abstract

We have found the first instance of a third-order topological Anderson
insulator (TOTAI). This disorder-induced topological phase is gapped and
characterized by a quantized octupole moment and topologically protected
corner states, as revealed by a detailed numerically exact analysis. We also
find that the disorder-induced transition into the TOTAI phase can be
analytically captured with remarkable accuracy using the self-consistent Born
approximation. For a larger disorder strength, the TOTAI undergoes a
transition to a trivial diffusive metal, that in turn becomes an Anderson
insulator at even larger disorder. Our findings show that disorder can induce
third-order topological phases in 3D, therefore extending the class of known
higher-order topological Anderson insulators.

##  I Introduction

In symmetry-protected topological (SPT) phases of matter, such as topological
insulators (TIs), non-trivial bulk topology leads to protected gapless
excitations on the system’s boundary [1, 2, 3, 4]. These edge-states have
exotic, disorder-robust properties with promising applications for quantum
computation [5, 6, 7]. SPT phases of matter are classified in the ten fold way
[8], based on the discrete symmetries (time-reversal, charge-conjugation and
chiral) that constrain the system’s Hamiltonian. Spatial symmetries of
crystalline nature may also be encountered, producing topological crystalline
insulators (TCIs) [9, 10]. Recently, TIs have been generalized to higher-order
topological insulators (HOTIs), where the bulk-boundary correspondence applies
to the $(d-n)$ dimensional boundary, for a $d$-dimensional, $n$th-order
topological insulator [11, 12, 13, 14, 15, 16]. HOTIs were first demonstrated
in the Benalcazar-Bernevig-Hughes (BBH) models [17, 18], where the topological
invariant corresponds to quantized bulk quadrupole or octupole electric
moments respectively in a 2D second-order topological insulator (SOTI) and 3D
third-order topological insulator (TOTI), with protected corner states. In the
BBH models, the topological properties are protected by spatial symmetries,
rendering them an extension of the TCIs.

Many experimental implementations of HOTIs have since been found, first in
classical metamaterial analogues like mechanical metamaterials [19], electric
circuits [20, 21], coupled microwave resonators [22], photonic waveguides
[23]; and later even in solid-state materials [24, 25, 26]. In any practical
realization of a system, disorder is present, e.g. due to defects in
manufacturing and can even be tuned in metamaterials. Disorder has a profound
impact on quantum transport due to Anderson localization of electronic wave
functions [27, 28]. This gives rise to Anderson insulators, that can have
gapless excitations in contrast with conventional (gapped) band insulators
[29]. It is generally known that TIs are robust against symmetry-preserving
disorder. Still, with enough disorder, it is possible to suppress topological
phases. Remarkably, increasing disorder can also induce topological
transitions from trivial to topological phases, giving rise to Topological
Anderson insulators (TAIs) [30, 31], which have been experimentlly realized
recently in different platforms [32, 33, 34].

The concept of TAIs was recently extended to higher-order topological Anderson
insulators (HOTAIs) in Ref. [35], where a 2D SOTI was obtained by adding
chiral-symmetric disorder to the 2D-BBH model. This result establishes chiral
symmetry as a sufficient symmetry to protect the HOTAI phases, even when the
crystalline symmetries are broken by disorder. A full phase diagram was
obtained in Ref. [36] for a system that can be mapped to the 2D-BBH model. It
was found that the disorder-induced SOTI comes in two varieties with
increasing disorder: the gapped and gapless HOTAI phases, followed by a
Griffiths phase. Noteworthy, the classical analogue of a 2D HOTAI was recently
experimentally observed using electric circuits [37], where disorder can be
tuned. A disorder-driven 3D SOTI was also found in amorphous systems with
structural disorder [38, 39, 40].

In this work, we find the first instance of a disorder-induced third-order
topological Anderson insulator (TOTAI). Our conclusions are drawn from from
the numerical analysis of the interplay between topology and chiral symmetry
preserving disorder in the 3D-BBH model. The TOTAI phase is gapped and
undergoes a transition into a trivial (gapless) diffusive metal (DM) with
increasing disorder. At significantly larger disorder, it turns into an
Anderson insulator (AI). The gapless HOTAI phase and the Griffiths phase are
absent, in contrast to the disordered 2D-BBH model [36]. The detailed
topological, spectral, and localization properties of the different phases
found are summarized in Fig. 2 and Tab. 1, and will be justified in detail
below.

This paper is structured as follows. In Sec. II, we present the model and the
topological invariants that characterize non-trivial phases, and which we
compute numerically. Detailed numerical results are presented in Sec. III,
which allowed for the full description of the phase diagram of the model. We
also analytically capture the disorder-induced topological phase transition
using the self-consistent Born approximation. In Sec. IV we discuss our
results and their implications.

##  II Model and methods

(a)(b)(c)

Figure 1:  (a,b) Schematics of the 3D-BBH model with disorder. In (a) only the
inter-cell hoppings are shown, whilst in (b) the intra-cell hoppings are
presented. Dotted lines correspond to negative signs in the clean hopping
amplitudes. (c) Schematics of full system with bulk octupole moment $o_{xyz}$,
boundary quadrupole moments $q_{ij}$ and size $L_{i}$ in directions
$i,j\in\\{x,y,z\\}$.

Model.— The model under consideration is the 3D-BBH model [17], generalized
with disorder in the intra-cell hopping amplitudes, as illustrated in Fig.
1(a,b). The corresponding tight-binding Hamiltonian is given by

$\hat{H}=\sum_{\textbf{r}}\left[\hat{c}_{\textbf{r}}^{\dagger}\Gamma_{\textbf{r}}\hat{c}_{\textbf{r}}+\sum_{i\in\\{x,y,z\\}}\left(\hat{c}^{\dagger}_{\textbf{r}}\Lambda_{i}\hat{c}_{\textbf{r}+\textbf{e}_{i}}+H.c.\right)\right]\,,$
(1)

where $\hat{c}^{\dagger}_{\textbf{r}}=(\hat{c}_{\textbf{r}1}^{\dagger}\
\hat{c}_{\textbf{r}2}^{\dagger}\ \dots\ \hat{c}_{\textbf{r}8}^{\dagger})$,
$\hat{c}_{\textbf{r}\alpha}^{\dagger}$ creates a particle at the $\alpha$-th
site of cell r, and the hopping matrices are given by

$\begin{split}[\Gamma_{\textbf{r}}]_{ij}&=\gamma_{\textbf{r}}^{ij}[\sigma_{z}\otimes(\sigma_{x}\otimes\mathds{1}-\sigma_{y}^{\otimes
2})+\sigma_{x}\otimes\mathds{1}^{\otimes 2}]_{ij}\,,\\\
\Lambda_{x}&=\frac{\lambda}{2}\mathds{1}\otimes(\sigma_{x}\otimes\mathds{1}+i\sigma_{y}\otimes\sigma_{z})\,,\\\
\Lambda_{y}&=\frac{\lambda}{2}\mathds{1}\otimes
i\sigma_{y}\otimes(\sigma_{x}+i\sigma_{y})\,,\\\
\Lambda_{z}&=\frac{\lambda}{2}(\sigma_{x}+i\sigma_{y})\otimes\mathds{1}^{\otimes
2}\,,\end{split}$ (2)

where $\\{\mathds{1},\sigma_{x},\sigma_{y},\sigma_{z}\\}$ is the set of
$2\times 2$ identity and Pauli matrices. We set $\lambda=1$ so that the energy
is measured in units of $\lambda$. The intra-cell hopping amplitudes are (up
to a sign as indicated in Fig. 1(b) and Eq. (2))
$\gamma_{\textbf{r}}^{ij}=\gamma+W\Delta^{ij}_{\textbf{r}}$, where $W$ is the
disorder strength and $\Delta_{\textbf{r}}^{ij}=\Delta_{\textbf{r}}^{ji}$ are
uniformly distributed random variables in the interval
$[-\frac{1}{2},\frac{1}{2}]$ without correlation. In our finite-size
calculations, we consider cubic systems of size $L_{x}=L_{y}=L_{z}=L$.

In the clean limit, $W=0$, we have $H\rightarrow H_{0}$,
$\Gamma_{\textbf{r}}\rightarrow\Gamma_{0}$ and translational invariance allows
us to express the Hamiltonian in reciprocal space as,

H0(k) =σz⊗[σx⊗1 (cos(kx)+ γ) - σy⊗σzsin(kx)]-σz⊗σy⊗[σy(cos(ky)+ γ) +
σxsin(ky)]+[σx(cos(kz)+ γ) - σysin(kz)]⊗1⊗2 .

(3)

The topological properties depend on the value of the parameter $\gamma$, as
discussed next.

Topological properties.— When $|\gamma|<1$ ($|\gamma|>1$), the clean
Hamiltonian in Eq. (3) is in a topological (trivial) phase with quantized
octupole moment $o_{xyz}$. In reciprocal space, $o_{xyz}$ may be computed by
the nested Wilson loop method, where the spatial reflection and inversion
symmetries of the clean system, along with time-reversal, charge-conjugation,
and chiral symmetries were shown to protect the topology [17]. In real space,
$o_{xyz}$ is computed through many-body electric multipole operators [41, 42,
43]. Since this involves finding the ground state of the system, it is
computationally demanding to do it in 3D. However, in the topological phase,
we also expect to find quantized quadrupole moments $q_{xy},q_{xz},q_{yz}$ in
the 2D-boundaries of the insulator, as illustrated in Fig. 1(c), allowing for
the definition of the topological invariant

$Q=8\left|q_{xy}q_{xz}q_{yz}\right|\ ,$ (4)

where each quadrupole moment is expressed as

$q_{ab}=\left[\frac{1}{2\pi}\textnormal{Im}\log\bra{\Psi_{c}}\mathcal{U}_{ab}\ket{\Psi_{c}}-q^{(0)}_{ab}\right]\textnormal{mod
}1\ ,$ (5)

with

$\mathcal{U}_{ab}=\exp\left(\frac{2\pi
i\sum_{j=1}^{N_{\textnormal{occ}}}\hat{r}^{j}_{a}\hat{r}^{j}_{b}}{L_{a}L_{b}}\right)\
,$ (6)

for $c\neq a\neq b$, where $\hat{r}_{a}^{j}$ is the position operator in
direction $a=x,y,z$ for electron $j$ and $N_{\textnormal{occ}}=2L_{a}L_{b}$
the number of occupied states in the boundary $ab$, with $L_{a}$ the number of
unit cells in direction $a$.
$q^{(0)}_{ab}=\frac{1}{2}\sum^{N_{\textnormal{a}}}_{j=1}r_{a}^{j}r_{b}^{j}/(L_{a}L_{b})$
is the contribution from the positive background charge, taking into account
that the sample is electrically neutral with
$N_{\textnormal{a}}=2N_{\textnormal{occ}}$ atomic orbitals in the boundary.
$\ket{\Psi_{c}}$ is the boundary many-body ground state obtained from the
effective Hamiltonian $H_{c}=-{G_{N_{c}}^{c}}(E=0)^{-1}$, with $N_{c}=2L_{c}$.
$G_{N_{c}}^{c}$ is the boundary Green’s function [44] that can be computed by
dividing the Hamiltonian matrix into 2D layers in the direction $c$ and
solving the following Dyson equation,

$G_{n}^{c}=(E-h_{n}^{c}-V_{n-1}^{c}G_{n-1}^{c}{V_{n-1}^{c\,\dagger}})\ ,$ (7)

where $h_{n}^{c}$ is the Hamiltonian of the $n$th-layer, and the $(n-1)$th-
layer couples to the $n$th-layer through matrix $V_{n-1}^{c}$. The reduced
Hilbert space dimensionality of each layer allows for reaching far larger
system sizes when computing $Q$ than by computing the bulk octupole moment
through many-body electric multipole operators.

In the disordered system, spatial crystalline symmetries are broken. However,
the system is still chiral symmetric, since it is decomposable in sub-latices
that do not possess any hopping terms between themselves. We will see that
this symmetry suffices to protect the topology. Chiral symmetry is also
preserved in the effective boundary Hamiltonian. The quadrupole moments are
known to be quantized by chiral symmetry [36], which means that, in each
realization of disorder, $Q$ is quantized to $0$ or $1$.

Spectral properties.— To study the spectral properties of the different
phases, we computed the energy gap using exact diagonalization and the density
of states (DOS), $\rho(E)=\frac{1}{D}\sum_{k=0}^{D-1}\delta(E-E_{k})$, where
$D$ is the Hilbert space dimension and ${E_{k}}$ are the single-particle
eigenenergies. For an efficient calculation of the DOS, we employed the kernel
polynomial method (KPM) [45]. In all our KPM calculations, we evaluated the
trace stochastically over a single random state and used the Jackson kernel. A
related quantity that can also be computed with the KPM is the local density
of states (LDOS),
$\rho(E,\textbf{r})=\sum_{k=0}^{D}\sum_{\alpha}\left|\psi_{k}(\textbf{r},\alpha)\right|^{2}\delta(E-E_{k})$,
where $\psi_{k}(\textbf{r},\alpha)$ is the $k$th eigenfunction evaluated at
unit cell r and orbital $\alpha$. We used this quantity to inspect the
existence of localized corner states, to complement the analysis on the
topological properties.

Localization properties.— Finally, we also study the localization properties
of the eigenstates by evaluating their localization length, the average level-
spacing ratio (LSR), the inverse participation ratio (IPR) and the fractal
dimension.

The normalized localization length $\Lambda=\lambda/L$, where $\lambda$ is the
localization length along the $z$ direction and $L=L_{x}=L_{y}$, was computed
using the transfer matrix method (TMM) [46, 47]. For extended states,
$\Lambda$ increases with $L$, while for localized states, $\Lambda\rightarrow
0$, since $\lambda$ is finite. At critical points, $\lambda\sim L$ and
therefore $\Lambda\sim L^{0}$.

The LSR is given by

$\textnormal{LSR}=\frac{1}{n-2}\sum_{i=1}^{n-2}\frac{\min(\delta_{i},\delta_{i+1})}{\max(\delta_{i},\delta_{i+1})}\
,$ (8)

where $\delta_{i}=E_{i+1}-E_{i}$ are the spacings between $n$ eigenenergies
$E_{i}$ sorted in ascending order. The energy gap spacing is not included. We
expect the energy level spacings of localized eigenstates to follow Poisson
statistics, in which case $\textnormal{LSR}\approx 0.386$. For diffusive
extended states, the level spacings follow the Gaussian Orthonormal Ensemble
(GOE) probability distribution, corresponding to $\textnormal{LSR}\approx
0.530$ [48].

The IPR [49] is expressed as

$\textnormal{IPR}=\frac{1}{n}\sum_{i=1}^{n}\sum_{\textbf{r}}\left(\sum_{\alpha}\left|\psi_{i}(\textbf{r},\alpha)\right|^{2}\right)^{2}\
,$ (9)

where $\psi_{i}(\textbf{r},\alpha)$ is the amplitude of the $i$th
eigenfunction at unit cell r and orbital $\alpha$. The IPR scales with system
size as $\textnormal{IPR}\propto L^{-D_{2}}$, where $D_{2}$ is the (real-
space) fractal dimension given by $D_{2}=3$ for extended states, $D_{2}=0$ for
localized states and $0<D_{2}<3$ for fractal or multifractal states [29].

##  III Results

Figure 2:  Schematic phase diagram as a function of $W$ for $\gamma=1.1$.

(a)(b)(c)

Figure 3:  (a) Topological phase diagram obtained from the topological
invariant $Q$ defined in Eq. (4) with respect to the disorder strength $W$.
For the lines with fixed size, $Q$ was averaged with 40 disorder realizations.
To compute the extrapolated points at some selected values of $W$, $Q$ was
averaged over 400 disorder realizations. In (b), an example of the
extrapolation is shown for $W=3$. (c) Bulk energy gap computed from exact
diagonalization for a system size $L=20$ and averaging over $200$ disorder
realizations, as a function of $W$.

(a)(b)

Figure 4:  Local density of states at zero-energy as function of unit-cell
number r, in a corner of a system with size $L=30$. The kernel polynomial
method was used with a single random state trace approximation and $N=2^{10}$
moments, averaged over 200 disordered samples with disorder weight (a) $W=3$
(phase II) and (b) $W=4$ (phase III).

(a)(b)(c)

Figure 5: Density of states $\rho(E)$, computed with the kernel polynomial
method for a system size $L=80$. (a) $\rho(E=0)$ as function of disorder
strength $W$. The two curves shown are for different choices of the number of
Chebyshev moments $N$. (b) $\rho(E)$ computed with $N=2^{13}$ moments, for
selected values of $W$, with a zoomed-in view around the zero-energy region in
(c).

Phase | I: GI | II: TOTAI | III: DM | IV: AI
---|---|---|---|---
Topology | Trivial | Non-trivial | Trivial | Trivial
Spectrum | Gapped | Gapped | Gapless | Gapless
Zero-energy states | Localized | Localized | Extended | Localized
$W_{c}$ | $-$ | $2.55(20)$ | $3.54(3)$ | $24(2)$

Table 1: Summary of all the phases observed in the model for $\gamma=1.1$:
trivial gapped insulator (GI), third-order topological Anderson insulator
(TOTAI), diffusive metal (DM) and Anderson insulator (AI); with the respective
topological, spectral and localization properties.

Starting from a trivial insulator in the clean limit, $\gamma=1.1$, we found
four different phases as a function of disorder strength $W$, that are
summarized in Fig. 2 and Tab. 1. In the next sections, we detail the
properties of each phase.

###  III.1 Topological phase diagram

In Fig . 3(a), the phase diagram for the topological invariant $Q$ is shown.
Due to the large finite-size effects, extrapolations to $L\rightarrow\infty$
were performed. For the extrapolations, three linear fits were performed for
$Q(L^{-1})$, for the 5 largest values of $L$ (green), for the 10 largest
values of $L$ (orange) and also for all values of $L$ (red), as shown in Fig.
3(b). The extrapolated value of $Q\left(L^{-1}\rightarrow 0\right)$ is the
average result of the three fits. Starting from a topologically trivial phase
I, at the critical disorder $W_{c}^{\textnormal{II}}=2.55(20)$ an abrupt
increase in $Q$ occurs, indicating the start of phase II.
$W_{c}^{\textnormal{II}}$ is precisely determined by the lowest $W$ for which
$Q$ increases and a relatively large error is considered to take into account
finite-size effects. For this phase II, extrapolated values of $Q$ are
compatible with $Q=1$ within error bar, signaling a topologically non-trivial
phase. As shown in Fig. 3(c), the gap closes and reopens at
$W_{c}^{\textnormal{II}}$, further pointing to a transition into a
tropological phase. Since this phase was induced by disorder, we dubbed it a
TOTAI. However, it is important to note that the system is gapped in this
phase, as is evident from Fig. 3(c).

Further increasing disorder, the system transitions to phase III, where
extrapolated values of $Q$ are compatible with $Q=0$, indicating that it is
topologically trivial. These results are also compatible with the zero-energy
LDOS shown in Fig. 4, revealing the existence of localized protected corner
states in the TOTAI phase II, and their absence in the trivial phase III. We
estimated $W_{c}^{\textnormal{III}}$ from the extrapolation to
$L\rightarrow\infty$ (analogously to Fig. 3(b)), of the crossing between the
energy gap and the mean level spacing (not shown) of the 50 states closer to
$E=0$, resulting in $W_{c}^{\textnormal{III}}=3.51(3)$. This estimation will
be compared to a different one based on localization properties in section
III.3.

(a)(b)(c)(d)(e)(f)(g)(h)

Figure 6:  Normalized localization length $\Lambda$ from the transfer matrix
method at (a,d) E = 0 versus $W$ for distinct $L$, and versus $L$ for distinct
energies at (b) $W=3.4$ and at (c) $W=5$. (e) LSR, (f) IPR and (g) fractal
dimension $D_{2}$ for $n$ eigenstates around zero-energy from exact
diagonalization, as a function of $W$. (h) $D_{2}$ versus E for distinct $W$
for $n=10$ eigenstates around $E$. Averages were taken over 200 realizations
of disorder. In (e) and (f), $L=20$. $D_{2}$ was computed by fitting using the
sizes $L\in\\{10,12,\dots,20\\}$ in (f) and $L\in\\{4,6,\dots,16\\}$ in (g).

###  III.2 Density of States

In Fig. 5 we show the DOS for different disorder strengths. We note that at
the topological transition from phase I to II, although the gap closes,
$\rho(0)$ is always zero [see Fig. 5(a)], behaving as $\rho(E)\sim E^{2}$
around $E=0$ [see Fig. 5(c)], as it would for a clean system with a Dirac cone
(which is the case of the clean 3D-BBH model in the topological transition
point). This was verified by performing a linear fit in a log-log plot of the
curves in Fig. 5(b) of positive $E$ values close to $E=0$ for
$W\in\\{2.5,2.6\\}$, which rendered slopes compatible with two (not shown). In
phase III, the energy gap closes again and $\rho(0)$ becomes finite. The DOS
starts to become peaked around $E=0$ for large $W$. Whether this finite DOS at
the Fermi level ($E=0$) is associated with a diffusive metal or an Anderson
insulator is discussed next.

###  III.3 Localization properties

In Fig. 6(a), we plot the normalized localization length $\Lambda$ along the
$z$-direction at $E=0$. The calculations of $\Lambda$ along other directions
yielded quantitatively identical results. We can see that $\Lambda$ decreases
with $L$ in phases I and II. This is because the system is gapped and the wave
function can therefore only propagate through (evanescent) localized modes at
$E=0$. At the topological phase transition, however, $\Lambda$ becomes
$L$-independent, as expected. In phase III, the system is gapless and has
extended states at $E=0$ since $\Lambda$ increases with $L$, as expected for a
diffusive metal. We also note that for energies where the DOS is finite, the
eigenstates are extended in phases I-III, as supported in Figs 6(b,c).

For large $W$, we see another phase transition at
$W_{c}^{\textnormal{IV}}=24(2)$ to a phase IV where $\Lambda$ again decreases
with $L$, Fig. 6(d). In this case, even though the system is gapless, the bulk
extended states become localized at $E=0$. In fact, localization occurs at all
energies and corresponds to the standard Anderson transition [27, 28, 29].

In order to make an additional independent estimation of the critical point
for the transition from phase II to III, we also analyzed the crossing points
between curves of adjacent $L$ in Fig. 6(a). Fitting the crossing points
analogously to what was done in 3(b), we extrapolate
$W_{c}^{\textnormal{III}}=3.56(3)$ in the thermodynamic limit. This is
compatible with the result obtained in section III.1 and their average is
presented in Tab. 1.

We now turn to the LSR analysis. In Fig. 6(e), we present the LSR for
eigenenergies around $E=0$, where we had to disregard some abnormally large
outlier spacings created due to finite-size effects (they correspond to
spacings between sets of degenerate states in the clean limit). The LSR in
phase III follows GOE statistics, completing the proof that phase III is a
diffusive metal. In phases I and II, where we access the statistics of the gap
edge, the states mostly follow the GOE ensemble for diffusive and extended
states. However, as we approach the transition point $W^{\mathrm{III}}_{c}$,
there is a sudden decrease in the LSR, especially at lower $n$ (closer to the
gap edge). To better understand this result, we calculated the IPR and the
fractal dimension, which we discuss next.

For the gapped phases, we computed the average IPR for eigenstates at the gap
edge, as for the LSR. In Fig. 6(f), we can see that the IPR is small in phases
I and III, which, in combination with the obtained fractal dimension
$D_{2}\approx 3$ in Fig. 6(g), indicates that the eigenstates closer to $E=0$
are extended. We also observe in Fig. 6(f) that the IPR becomes larger in
phase II, peaking close to the transition II $\rightarrow$ III. This is
concomitant with the fractal dimension results in Fig. 6(g), where it can be
seen that $D_{2}\approx 0$ close to the transition, suggesting the presence of
localized gap-edge states right before the gap closes. This correlates with
the sudden drop of the LSR. However, there are still some discrepancies
between the results for the LSR and fractal dimension (the LSR is still
significantly away from Poisson), that we attribute to strong finite-size
effects in phase II. Fig. 6(h) further shows that in phase III the states are
extended for any energy, while in phase II the states are only localized close
to $E=0$, at the gap edges. These localized states are likely related with
Lifshitz tails, whose exponentially suppressed DOS in the thermodinamic limit
justifies the strong finite size effects, especially for the LSR results.

###  III.4 Self-Consistent Born Approximation

Figure 7:  Effective renormalizaded intra-cell hopping amplitude
$\gamma^{\prime}$, computed through the SCBA, as a function of the clean
hopping amplitude $\gamma$ and the disorder strength $W$. The topological
transition curve at $\gamma^{\prime}=1$ is shown in red and the transition
numerically extracted from the topological invariant $Q$ is shown in the blue
points.

Disorder is introduced into the system in the form of added intra-cell hopping
amplitudes at each unit cell r, that is,

$V_{\textbf{r}}=\sum_{\alpha=1}^{12}V_{\textbf{r},\alpha}U_{\alpha}\,,$ (10)

where $\alpha(i,j)\in\\{1,\dots,12\\}$ is a bijection between the index of an
edge $\alpha$ and the indexes $i,j$ of the adjacent corners.
$V_{\textbf{r},\alpha(i,j)}=W\Delta_{r}^{ij}$ are the hopping strengths and

$\left[U_{\alpha(i,j)}\right]_{mn}=\frac{1}{\gamma}\left[\Gamma_{0}\right]_{mn}(\delta_{mi}\delta_{nj}+\delta_{mj}\delta_{ni})$
(11)

are the matrix elements of each separate intra-cell hopping term. Since the
disorder is uncorrelated,

$\langle V_{\textbf{r},\alpha}\rangle=0\ ,\ \langle
V_{\textbf{r},\alpha}V_{\textbf{r}^{\prime},\beta}\rangle=\frac{W^{2}}{12}\delta_{\textbf{r}\textbf{r}^{\prime}}\delta_{\alpha\beta}\
.$ (12)

Under the Self-Consistent Born approximation (SCBA)[50, 31, 36, 51], the
effective Bloch Hamiltonian at $E=0$ is
$H_{\textnormal{eff}}(\textbf{k})=H_{0}(\textbf{k})+\Sigma(E=0)$, where the
self-energy $\Sigma$ is computed self-consistently through the following
equation,

$\Sigma(E)=\frac{W^{2}}{12(2\pi)^{3}}\int_{BZ}d^{3}\textbf{k}\sum_{\alpha=1}^{12}U_{\alpha}GU_{\alpha}\,,$
(13)

where $G=\left[(E+i0^{+})\mathds{1}-H_{0}(\textbf{k})-\Sigma(E)\right]^{-1}$
is the Green’s function. Numerically, we find that
$\Sigma(0)=-\Gamma_{0}\sigma/\gamma$, $\sigma\in\mathbb{R}$. In the effective
Hamiltonian, this amounts to a normalization of the intra-cell hopping
strengths $\gamma\rightarrow\gamma^{\prime}=\gamma-\sigma$. Since the
effective model still corresponds to the 3D-BBH clean model, the topological
(trivial) phase occurs for $\gamma^{\prime}<1(>1)$. In Fig. 7, we observe that
the topological transition curve predicted by the SCBA agrees very well with
the one computed numerically from the topological invariant $Q$.

##  IV Conclusions

In summary, we have discovered the first example of a third-order topological
Anderson insulator, induced by chiral symmetry preserving disorder. The TOTAI
phase is characterized by a quantized quadrupole moment on the boundaries of
the 3D system, that corresponds to a quantized bulk octupole moment, and by
topologically protected localized corner states. Remarkably, the topological
transition to the TOTAI phase is captured with great accuracy by the self-
consistent Born approximation, up to very large disorder strengths.

Our findings can be tested experimentally in different metamaterials where
disorder can be tuned, such as mechanical metamaterials [19], electric
circuits [20, 37, 21] or photonic waveguides [23].

Finally, we note that in contrast to the disordered 2D BBH model [36], we have
not found a gapless HOTAI in 3D. This raises an interesting open question for
future research: do gapless TOTAIs exist?

##  V Acknowledgments

This work has been partially funded by the ERC Starting Grant 101042293
(HEPIQ) (H.L.). The authors MG and PR acknowledge partial support from
Fundação para a Ciência e Tecnologia (FCT-Portugal) through Grant No.
UID/CTM/04540/2019. EVC acknowledge partial support from FCT-Portugal through
Grant No. UIDB/04650/2020. MG acknowledges further support from FCT-Portugal
through the Grant SFRH/BD/145152/2019. We finally acknowledge the Tianhe-2JK
cluster at the Beijing Computational Science Research Center (CSRC), the
Bob|Macc supercomputer through computational project project
CPCA/A1/470243/2021 and the OBLIVION supercomputer, through projects
HPCUE/A1/468700/2021, 2022.15834.CPCA.A1 and 2022.15910.CPCA.A1 (based at the
High Performance Computing Center - University of Évora) funded by the ENGAGE
SKA Research Infrastructure (reference POCI-01-0145-FEDER-022217 - COMPETE
2020 and the Foundation for Science and Technology, Portugal) and by the
BigData@UE project (reference ALT20-03-0246-FEDER-000033 - FEDER and the
Alentejo 2020 Regional Operational Program. Computer assistance was provided
by CSRC’s, Bob|Macc’s and OBLIVION’s support teams.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Reviews of modern physics 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Reviews of Modern Physics 83, 1057 (2011).
  * Moore [2010] J. E. Moore, Nature 464, 194 (2010).
  * Shen [2012] S.-Q. Shen,  _Topological insulators_, Vol. 174 (Springer, 2012).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Reviews of Modern Physics 80, 1083 (2008).
  * Freedman  _et al._ [2003] M. Freedman, A. Kitaev, M. Larsen, and Z. Wang, Bulletin of the American Mathematical Society 40, 31 (2003).
  * Stern and Lindner [2013] A. Stern and N. H. Lindner, Science 339, 1179 (2013).
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. Ludwig, Physical Review B 78, 195125 (2008).
  * Fu [2011] L. Fu, Physical Review Letters 106, 106802 (2011).
  * Neupert and Schindler [2018] T. Neupert and F. Schindler, Topological Matter: Lectures from the Topological Matter School 2017 , 31 (2018).
  * Schindler  _et al._ [2018a] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Science advances 4, eaat0346 (2018a).
  * Khalaf [2018] E. Khalaf, Physical Review B 97, 205136 (2018).
  * Park  _et al._ [2019] M. J. Park, Y. Kim, G. Y. Cho, and S. Lee, Physical review letters 123, 216803 (2019).
  * Wang  _et al._ [2021a] B. Wang, X. Zhou, H. Lin, and A. Bansil, Physical Review B 104, L121108 (2021a).
  * Franca  _et al._ [2018] S. Franca, J. van den Brink, and I. Fulga, Physical Review B 98, 201114 (2018).
  * Chen  _et al._ [2020] R. Chen, C.-Z. Chen, J.-H. Gao, B. Zhou, and D.-H. Xu, Physical review letters 124, 036803 (2020).
  * Benalcazar  _et al._ [2017a] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Benalcazar  _et al._ [2017b] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Physical Review B 96, 245115 (2017b).
  * Serra-Garcia  _et al._ [2018] M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Imhof  _et al._ [2018] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nature Physics 14, 925 (2018).
  * Bao  _et al._ [2019] J. Bao, D. Zou, W. Zhang, W. He, H. Sun, and X. Zhang, Phys. Rev. B 100, 201406 (2019).
  * Peterson  _et al._ [2018] C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature 555, 346 (2018).
  * Dutt  _et al._ [2020] A. Dutt, M. Minkov, I. A. D. Williamson, and S. Fan, Light: Science & Applications 9, 131 (2020).
  * Schindler  _et al._ [2018b] F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Nature Physics 14, 918 (2018b).
  * Shumiya  _et al._ [2022] N. Shumiya, M. S. Hossain, J.-X. Yin, Z. Wang, M. Litskevich, C. Yoon, Y. Li, Y. Yang, Y.-X. Jiang, G. Cheng, _et al._ , Nature Materials 21, 1111 (2022).
  * Lee  _et al._ [2023] J. Lee, J. Kwon, E. Lee, J. Park, S. Cha, K. Watanabe, T. Taniguchi, M.-H. Jo, and H. Choi, Nature Communications 14, 1801 (2023).
  * Anderson [1958] P. W. Anderson, Physical review 109, 1492 (1958).
  * Lagendijk  _et al._ [2009] A. Lagendijk, B. Van Tiggelen, and D. S. Wiersma, Phys. today 62, 24 (2009).
  * Evers and Mirlin [2008] F. Evers and A. D. Mirlin, Reviews of Modern Physics 80, 1355 (2008).
  * Li  _et al._ [2009] J. Li, R.-L. Chu, J. K. Jain, and S.-Q. Shen, Physical review letters 102, 136806 (2009).
  * Groth  _et al._ [2009] C. Groth, M. Wimmer, A. Akhmerov, J. Tworzydło, and C. Beenakker, Physical review letters 103, 196805 (2009).
  * Meier  _et al._ [2018] E. J. Meier, F. A. An, A. Dauphin, M. Maffei, P. Massignan, T. L. Hughes, and B. Gadway, Science 362, 929 (2018).
  * Stützer  _et al._ [2018] S. Stützer, Y. Plotnik, Y. Lumer, P. Titum, N. H. Lindner, M. Segev, M. C. Rechtsman, and A. Szameit, Nature 560, 461 (2018).
  * Liu  _et al._ [2020] G.-G. Liu, Y. Yang, X. Ren, H. Xue, X. Lin, Y.-H. Hu, H.-x. Sun, B. Peng, P. Zhou, Y. Chong, and B. Zhang, Physical Review Letters 125, 133603 (2020).
  * Li  _et al._ [2020] C. A. Li, B. Fu, Z. A. Hu, J. Li, and S. Q. Shen, Physical Review Letters 125, 166801 (2020), 2008.00513 .
  * Yang  _et al._ [2021] Y.-B. Yang, K. Li, L.-M. Duan, and Y. Xu, Physical Review B 103, 085408 (2021).
  * Zhang  _et al._ [2021] W. Zhang, D. Zou, Q. Pei, W. He, J. Bao, H. Sun, and X. Zhang, Phys. Rev. Lett. 126, 146802 (2021).
  * Agarwala  _et al._ [2020] A. Agarwala, V. Juričić, and B. Roy, Physical Review Research 2, 012067 (2020).
  * Wang  _et al._ [2021b] J.-H. Wang, Y.-B. Yang, N. Dai, and Y. Xu, Physical Review Letters 126, 206404 (2021b).
  * Peng  _et al._ [2022] T. Peng, C.-B. Hua, R. Chen, Z.-R. Liu, H.-M. Huang, and B. Zhou, Physical Review B 106, 125310 (2022).
  * Wheeler  _et al._ [2019] W. A. Wheeler, L. K. Wagner, and T. L. Hughes, Physical Review B 100, 245135 (2019).
  * Kang  _et al._ [2019] B. Kang, K. Shiozaki, and G. Y. Cho, Physical Review B 100, 245134 (2019).
  * Resta [1998] R. Resta, Physical Review Letters 80, 1800 (1998).
  * Peng  _et al._ [2017] Y. Peng, Y. Bao, and F. von Oppen, Physical Review B 95, 235143 (2017).
  * Weiße  _et al._ [2006] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Reviews of modern physics 78, 275 (2006).
  * Hoffmann and Schreiber [2002] K.-H. Hoffmann and M. Schreiber,  _Computational Statistical Physics: From Billards to Monte-Carlo_ (Springer Science & Business Media, 2002).
  * Kramer and Schreiber [1996] B. Kramer and M. Schreiber, Computational Physics: Selected Methods Simple Exercises Serious Applications , 166 (1996).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Physical review b 75, 155111 (2007).
  * Wegner [1980] F. Wegner, Zeitschrift für Physik B Condensed Matter 36, 209 (1980).
  * Bruus and Flensberg [2004] H. Bruus and K. Flensberg,  _Many-body quantum theory in condensed matter physics: an introduction_ (OUP Oxford, 2004).
  * Mašek [1986] J. Mašek, Zeitschrift für Physik B Condensed Matter 64, 145 (1986).