# Data-Driven Fast Frequency Control using Inverter-Based Resources

E. Ekomwenrenren  Dept. Electrical & Comp. Eng.
University of Waterloo
Waterloo, ON, Canada  J. W. Simpson-Porco  Dept. Electrical & Comp. Eng.
University of Toronto
Toronto, ON, Canada  E. Farantatos, M. Patel, A. Haddadi, L. Zhu This work was
supported by EPRI Project #10009168: Wide-Area Hierarchical Frequency and
Voltage Control for Next Generation Transmission Grids, and by NSERC Discovery
Grant RGPIN-2017-04008. Grid Operations & Planning Division
Electric Power Research Institute
Palo Alto, CA, USA

###### Abstract

We develop and test a data-driven and area-based fast frequency control
scheme, which rapidly redispatches inverter-based resources to compensate for
local power imbalances within the bulk power system. The approach requires no
explicit system model information, relying only on historical measurement
sequences for the computation of control actions. Our technical approach fuses
developments in low-gain estimator design and data-driven control to provide a
model-free and practical solution for fast frequency control. Theoretical
results and extensive simulation scenarios on a three area system are provided
to support the approach.

###### Index Terms:

renewable energy, frequency control, transmission grid, smart grid

##  I Introduction

Driven by the need to decarbonize the existing power system generation
infrastructure [1], the transmission grid is increasingly being dominated by
inverter-based renewable energy resources (IBRs). Challenges arising from this
transition away from traditional generation include larger (and more frequent)
frequency deviations, faster frequency dynamics due to reduced system inertia,
and increased net load variability [1, 2]. As a result of this this
proliferation of IBRs in the grid, there is an increasing urgency to develop
new and faster frequency control methods.

To address the aforementioned challenges, in [3] the authors proposed a
hierarchical control scheme which coordinates IBRs to provide fast frequency
control. In this scheme, the bulk grid is partitioned into geographically
small local control areas (LCAs). Within each LCA, a local controller is
designed which processes (potentially, delayed) measurements from the LCA to
compute updated set-points for local IBRs (Figure 1). The local controller has
two key sub-blocks. The _disturbance estimator_ detects frequency events
(e.g., a large load disturbance and/or generation outage) by computing a real-
time estimate $\Delta\widehat{P}_{\rm u}$ of the net unmeasured active power
imbalance $\Delta P_{\rm u}$ in the LCA. The _power allocator_ then
continuously and optimally redispatches the IBRs to correct the imbalance,
subject to device limits. In situations where local IBR resources are
insufficient, a higher-level coordinating controller facilitates the optimal
transfer of additional power support from neighbouring LCAs. As our focus in
this paper will be on the LCA controller design, this coordinating control
layer will not be discussed further; see [3, Sec. III-B] for details.

Local Control Area Disturbance
---
Estimator
Power
---
Allocator
LCA
ControllerIBRs$z^{-\tau_{\mathrm{m}}}$$z^{-\tau_{\mathrm{m}}}$$z^{-\tau_{\mathrm{m}}}$$z^{-\tau_{\mathrm{c}}}$$\Delta\widehat{P}_{\rm
u}$$\Delta P_{\mathrm{tie}}$$\Delta P_{\mathrm{u}}$$\Delta f$$\Delta
P_{\mathrm{ibr}}^{\mathrm{c}}$$\Delta P_{\mathrm{ibr}}$ Figure 1: Block
diagram of area control structure for each LCA. Dashed lines denote sampled
signals.

The key to the fast operation of this scheme lies in disturbance estimator,
which must quickly produce an accurate estimate of the power imbalance within
the LCA based on the available measurements. In [3], the disturbance estimator
was designed in a model-based fashion, using a crude second-order lumped
dynamic LCA model [4]

$\begin{aligned}
2H\Delta{\dot{\omega}}&=-\tfrac{1}{R_{\mathrm{I}}}\Delta\omega+\Delta
P_{\mathrm{m}}-\Delta P_{\mathrm{u}}-\Delta P_{\mathrm{tie}}+\Delta
P_{\mathrm{ibr,tot}}^{\mathrm{c}}\\\ T_{\rm
R}\Delta{\dot{P}_{\mathrm{m}}}&=-\Delta P_{\mathrm{m}}-R_{\rm
g}^{-1}(\Delta{\omega}+T_{\mathrm{R}}F_{\mathrm{H}}\Delta{\dot{\omega}})\end{aligned}$

(1)

parameterized by several constants, such as total LCA inertia $H$, total IBR
and generator primary control gains $R_{\rm I}$ and $R_{\rm g}$, aggregated
turbine-governor time constant $T_{\rm R}$, and aggregated high-pressure
turbine fraction $F_{\rm H}$. This modelling decision was made for pragmatic
reasons, as more accurate high-order system models could be too burdensome to
build and maintain in practice [5]. Tests in [3] demonstrated that this design
approach can yield excellent closed-loop control when the model (1) reasonably
captures the LCA dynamics and when the lumped parameters in (1) are accurately
set. However, when either of these conditions fails, control performance
(e.g., post-disturbance settling time and overshoot) will indeed deteriorate.

Many modern approaches to frequency control, such as robust optimal control
[6], model predictive control [7, 8, 9, 10] and the coordinated dispatch of
IBRs [3], are of similar vein in that they require accurate model information
about the power system, which may difficult to obtain in practice [5]. In
contrast to model-based approaches, a model-free or data-driven approach would
allow for better adaptation to realistic system conditions, resulting in
faster disturbance estimation and rejection, and with reduced tuning of
parameters required for real-world implementation. Our objective in this paper
is to develop such a data-driven frequency controller.

In the literature, various researchers have tackled the problem of power
system frequency control using data-driven and learning-based approaches. The
most popular approach is based on reinforcement learning, wherein control
actions are undertaken to maximize some form of cumulative reward [11, 12, 13,
14]. Drawbacks of this approach are that it requires a significant number of
data samples, and can be particularly sensitive to the selection of hyper-
parameters [15]. Other researchers take a hybrid approach, fusing elements of
data-driven control with sequential system identification (ID) and control
[16, 17]. However, since a system identification step is needed to identify an
approximate model, the issue of non-trivial parameter tuning still remains.

With the goal of using a smaller number of data samples and circumventing the
need of a system identification step, we base our approach here for the design
of a frequency controller on recent data-driven predictive control approaches
proposed in the literature, which compute control actions without explicitly
identifying the system from data [18, 19, 20, 21, 22, 23]. These design
techniques are based on a fundamental result [24] in the behavioral approach
to systems and control theory [25], namely that a single recorded historical
trajectory of a controllable LTI system implicitly defines a model of the
system, as long as the input sequence is persistently exciting (i.e.,
sufficiently rich). Authors in [23] utilize this technique in a data-driven
model predictive control framework to regulate frequency in the power system
context. However, in their approach, the net-load demand in the system is
assumed available and used as a control input; this information may be
difficult to obtain precisely at all times, especially with the added net-load
variability and uncertainty in power supply introduced by renewables [26].

_Contributions:_ The main contributions of this paper are summarized as
follows. First, in Section III we combine ideas from dynamic estimator theory
and behavioural systems theory to design a data-based substitute for the
model-based disturbance estimator in [3]. Second, in Section IV-A, we leverage
the theory of Section III to propose a novel model-free frequency controller,
which achieves fast frequency regulation by rapidly redispatching inverter-
based resources to compensate for local power imbalances within the bulk power
system. The key advantage of our approach is that it does not rely on a
parametric system representation for computing control actions; the LCA system
model in Figure 1 is replaced by time-series data, which is sufficiently rich
in harmonic content to capture the dynamics of the system. This time-series
data is then directly used in our hierarchical control framework, without
passing through an explicit system identification step. Third and finally, in
Section IV-C, we extensively validate our approach via simulations on a
detailed nonlinear three-area power system presented in Section IV-B. Across
several scenarios — including load increases, heavy-renewable penetration,
generation trips, and three-phase faults — we illustrate how our data-driven
approach can provide fast and effective frequency control for the bulk grid.

_Notation:_ Given vectors or matrices $x_{1},\ldots,x_{N}$ which have equal
numbers of columns, $\mathrm{col}(x_{1},\ldots,x_{N})$ denotes their column
concatenation. For a matrix $A$ having full column rank, $A^{\dagger}=(AA^{\sf
T})^{-1}A^{\sf T}$ denotes its pseudoinverse.

##  II Control Preliminaries

Our approach relies on a technical foundation built from two topics in linear
control theory: disturbance estimator design, and behavioral systems theory.
Section II-A recalls some preliminaries on disturbance estimator design.
Section II-B provides additional insights and assumptions that simplify the
disturbance estimator problem introduced in Section II-A. Section II-C
introduces the basic facts from behavioral systems theory that enable data-
driven prediction based on historical system trajectories.

###  II-A Background on Disturbance Estimators

Consider the discrete-time linear time-invariant (LTI) model

$\displaystyle x(t+1)$ $\displaystyle=Ax(t)+Bu(t)+B_{d}d(t)$ (2)
$\displaystyle y(t)$ $\displaystyle=Cx(t)+Du(t)$

with time index $t\in\mathbb{Z}_{\geq 0}$, state $x(t)\in^{n}$, control input
$u(t)\in^{m}$, disturbance input $d(t)\in^{q}$, and output $y(t)\in^{p}$. We
assume that $d(t)$ is a _constant_ but unknown disturbance signal and that
$q\leq p$. The problem of _asymptotic disturbance estimation_ is to design a
causal system which processes $(u(t),y(t))$ to produce an estimate
$\hat{d}(t)$ of $d(t)$ at time $t$, which satisfies
$\lim_{t\to\infty}(\hat{d}(t)-d(t))=0$ irrespective of the initial conditions.

This problem admits a standard model-based solution via the design of a
Luenberger observer for the extended model

$\displaystyle\xi(t+1)$ $\displaystyle=\underbrace{\begin{bmatrix}A&B_{d}\\\
0&I_{q}\end{bmatrix}}_{\coloneqq\mathcal{A}}\xi(t)+\underbrace{\begin{bmatrix}B\\\
0\end{bmatrix}}_{\coloneqq\mathcal{B}}u(t)$ (3) $\displaystyle y(t)$
$\displaystyle=\underbrace{\begin{bmatrix}C&0\end{bmatrix}}_{\coloneqq\mathcal{C}}\xi(t)+\underbrace{D}_{\coloneqq\mathcal{D}}u(t).$

where $\xi(t)=\mathrm{col}(x(t),d(t))$ is the extended state vector. In (3),
the dynamic disturbance model $d(t+1)=d(t)$ with unknown initial condition is
equivalent to the constant disturbance in the model (2). A standard
predictive-type Luenberger observer for (3) can now be designed as

$\displaystyle\hat{\xi}(t+1)$
$\displaystyle=\mathcal{A}\hat{\xi}(t)+\mathcal{B}u(t)-\mathcal{L}(\hat{y}(t)-y(t))$
(4) $\displaystyle\hat{y}(t)$ $\displaystyle=\mathcal{C}\hat{\xi}(t)+Du(t)$

where $\mathcal{L}\in^{(n+q)\times p}$ is the estimator gain to be designed.
With $\tilde{\xi}(t)=\xi(t)-\hat{\xi}(t)$ denoting the estimation error, the
error dynamics are computed to be

$\tilde{\xi}(t+1)=(\mathcal{A}-\mathcal{L}C)\tilde{\xi}(t).$ (5)

It is a standard result that there exists $\mathcal{L}$ such that the error
dynamics (5) are exponentially stable if and only if $(A,C)$ is detectable and
if the plant (2) has no transmission zeros at $z=1$ on the $d\mapsto y$
channel (see, e.g., [27, Chp. 23]. Under these conditions, the estimator (4)
solves the asymptotic disturbance estimation problem, with the estimate
$\hat{d}(t)$ being obtained from the second block component of $\hat{\xi}(t)$.

###  II-B Low-Gain Disturbance Estimation

The design of the estimator gain $\mathcal{L}$ above requires the solution of
a standard multivariable stabilization problem for the detectable pair
$(\mathcal{A},\mathcal{C})$ of extended plant matrices. In particular then,
the design requires knowledge of the entire system matrix $A$, describing the
internal state dynamics of the model. Under additional assumptions, this
design requirement can be significantly relaxed and made more practical.

First, we will assume that $A$ is Schur stable, i.e., that the plant (2) is
internally exponentially stable. This assumption would be satisfied by, for
example, a linearized power system model describing frequency dynamics and
incorporating typical primary controllers. Under this assumption $(A,C)$ is
automatically detectable, and the transmission zero condition reduces to the
condition that $G_{d}(1)\in^{p\times q}$ have full column rank, where
$G_{d}(z)=C(zI-A)^{-1}B_{d}$ denotes the transfer matrix on the $d\to y$
channel. In this case, if one seeks an estimator gain of the form

$\mathcal{L}=\begin{bmatrix}0\\\ \varepsilon L\end{bmatrix}$

for some $L\in^{q\times p}$ and $\varepsilon>0$, then (4) reduces to

$\displaystyle\hat{x}(t+1)$ $\displaystyle=A\hat{x}(t)+Bu(t)+B_{d}\hat{d}(t)$
(6a) $\displaystyle\hat{y}(t)$ $\displaystyle=C\hat{x}(t)+Du(t)$ (6b)
$\displaystyle\hat{d}(t+1)$ $\displaystyle=\hat{d}(t)-\varepsilon
L(\hat{y}(t)-y(t)).$ (6c)

The estimator (6) admits the following simple interpretation. The first two
equations (6a)–(6b) perform an open-loop simulation of the plant, using the
control input $u(t)$ and the _estimated_ disturbance $\hat{d}(t)$ to produce
an estimate $\hat{y}(t)$ for the output. Given $\hat{y}(t)$ and the true
output measurement $y(t)$, the equation (6c) updates the disturbance estimate.
Crucial to our data-driven development to follow is the following observation:
the only purpose of (6a)–(6b) is to produce the output estimate $\hat{y}(t)$;
the state estimate $\hat{x}(t)$ is irrelevant for our purposes.

The following result, which will not be proven here, guarantees that the
design of such an estimator is possible.

######  Proposition II.1 (Low-Gain Disturbance Estimator Design)

Consider the disturbance estimator (6) for the plant (2). Assume that $A$ is
Schur stable, that $G_{d}(1)=C(I_{n}-A)^{-1}B_{d}$ has full column rank, and
set $L=G_{d}(1)^{\dagger}$. Then there exists $\varepsilon^{\star}>0$ such
that for all $\varepsilon\in(0,\varepsilon^{\star})$, the estimator (6) solves
the asymptotic disturbance estimation problem.

Proposition II.1 states that with the simple estimator gain selection
$L=G_{d}(1)^{\dagger}$, one can always obtain a stable estimator by starting
$\varepsilon>0$ at a small value and then tuning. The quantity $G_{d}(1)$ is
known as the DC gain, and can be obtained directly from recorded data (See
Section III for details). The above result also applies without changes to
modified estimator

$\displaystyle\hat{x}(t+1)$ $\displaystyle=A\hat{x}(t)+Bu(t)+B_{d}\hat{d}(t)$
(7a) $\displaystyle\hat{y}(t)$ $\displaystyle=C\hat{x}(t)+Du(t)$ (7b)
$\displaystyle\hat{d}(t)$ $\displaystyle=\hat{d}(t-1)-\varepsilon
L(\hat{y}(t)-y(t)).$ (7c)

which uses the most recent measurement $y(t)$ to compute $\hat{d}(t)$ as
opposed to (6c) which uses $y(t-1)$ to compute $\hat{d}(t)$.

###  II-C Primer on Behavioral Systems Theory

The model defined by (2) is a parametric representation of a LTI system. In
the setting of behavioral systems theory, a LTI system is instead interpreted
as defining an implicit constraint on the vector spaces of input and output
signals. Our treatment here will be very minimal; see [25] for a detailed
introduction.

Let $\sigma$ denote the signal shift operator defined by $(\sigma
x)(t)=x(t+1)$. The _behavior_ $\mathscr{B}$ of (2) is defined as the set of
all possible input-output sequences which are consistent with the model

$\displaystyle\mathscr{B}=$
$\displaystyle\left\\{(u,d,y)\in(^{m+q+p})^{\mathbb{Z}_{\geq
0}}\,\,:\,\,\exists x\in(^{n})^{\mathbb{Z}_{\geq 0}}\,\,\text{s.t.}\right.$
$\displaystyle\qquad\left.\sigma x=Ax+Bu+B_{d}d,\,\,y=Cx+Du\right\\}.$

The smallest possible state dimension consistent with $\mathscr{B}$ is called
the _order_ of the system, and is denoted by $n(\mathscr{B})$. The _lag_ of
$\mathscr{B}$, denoted by $\ell(\mathscr{B})$ is the smallest integer $\ell$
such that the matrix
$\mathcal{O}_{\ell}=\mathrm{col}(C,CA,\ldots,CA^{\ell-1})$ has rank
$n(\mathscr{B})$. As additional notation, we let $\mathscr{B}_{N}$ denote the
restriction of the behaviour to trajectories of length $N\in\mathbb{Z}_{\geq
1}$.

The following notion is essential. Let $T\in\mathbb{Z}_{\geq 1}$ and let
$z\in(^{m})^{T}$ be the length $T$ signal

$z=\mathrm{col}(z(1),\ldots,z(T)).$

We say $z$ is _persistently exciting of order $L$_ if the _Hankel matrix_

$\mathscr{H}_{L}(z)=\begin{bmatrix}z(1)&\cdots&z(T-L+1)\\\
\vdots&\ddots&\vdots\\\ z(L)&\cdots&z(T)\end{bmatrix}\in^{mL\times(T-L+1)}$

has full row rank. The idea is that if $z$ is persistently exciting, then it
is both sufficiently long and sufficiently rich in variation.

Suppose that, from the plant (2), we have collected $T\in\mathbb{Z}_{\geq 1}$
samples of input and output data

$\displaystyle u_{\rm d}$
$\displaystyle=\mathrm{col}(u(1),\ldots,u(T))\in(^{m})^{T}$ (8) $\displaystyle
d_{\rm d}$ $\displaystyle=\mathrm{col}(d(1),\ldots,d(T))\in(^{q})^{T}$
$\displaystyle y_{\rm d}$
$\displaystyle=\mathrm{col}(y(1),\ldots,y(T))\in(^{p})^{T}.$

A key result known as the _fundamental lemma_ [20] provides a characterization
of _all possible_ length $N$ system trajectories in terms of recorded data.

######  Lemma II.2 (Fundamental Lemma)

Assume that $\mathscr{B}$ is controllable and that $\mathrm{col}(u_{\rm
d},d_{\rm d})$ is persistently exciting of order $N+n(\mathscr{B})$. Then any
possible length $N$ trajectory $(u,d,y)$ of $\mathscr{B}$ can be represented
as

$\begin{bmatrix}\mathscr{H}_{N}(u_{\rm d})\\\ \mathscr{H}_{N}(d_{\rm d})\\\
\mathscr{H}_{N}(y_{\rm d})\end{bmatrix}g=\begin{bmatrix}u\\\ d\\\
y\end{bmatrix}$

for some vector $g\in^{T-N+1}$.

Put simply, _all possible_ system trajectories of length $N$ can be expressed
as a linear combination of the columns of the given matrix, which consists
purely of recorded historical data.

Lemma II.2 enables data-driven simulation of a LTI system [20]. To obtain a
unique response, a form of initilization is required. This is done through an
initial prefix trajectory of length $T_{\rm ini}$ is required for
initialization; the technical requirement is that $T_{\rm
ini}\geq\ell(\mathscr{B})$. Let $(u_{\rm p},d_{\rm p},y_{\rm p})$ be this
prefix trajectory. Let

$\begin{bmatrix}U_{\rm p}\\\ U_{\rm f}\end{bmatrix}=\mathscr{H}_{T_{\rm
ini}+N}(u_{\rm d}),\qquad$

be a corresponding partitioning of the Hankel matrices, with similar
partitionings for the data $d_{\rm d}$ and $u_{\rm d}$. Leveraging Lemma II.2,
and assuming that $u_{\rm d}$ is persistently exciting of order $T_{\rm
ini}+N+n(\mathscr{B})$, we have the following: given length $N$ inputs
$u=(u(T_{\rm ini}+1),\ldots,u(T_{\rm ini}+N)$ and $d=(d(T_{\rm
ini}+1),\ldots,d(T_{\rm ini}+N)$, one may uniquely compute the response
$y=(y(T_{\rm ini}+1),\ldots,y(T_{\rm ini}+N))$ by solving the system of linear
equations

$\begin{bmatrix}U_{\rm p}\\\ D_{\rm p}\\\ Y_{\rm p}\\\ U_{N}\\\ D_{N}\\\
Y_{N}\end{bmatrix}g=\begin{bmatrix}u_{\rm p}\\\ d_{\rm p}\\\ y_{\rm p}\\\ u\\\
d\\\ y\end{bmatrix}$ (9)

in the unknowns $g$ and $y$. Specifically, one solves the first five block
equations for $g$, then substitutes into the last equation to compute the
response $y$.

##  III Data-Driven Disturbance Estimators

We now combine the disturbance estimation ideas Section II-B with the
behavioral systems theory from Section II-C to develop a data-driven version
of the disturbance estimator (7).111A data-driven version of the estimator (6)
can be similarly obtained by modifying the results that follow. As in Section
II-B, we assume that $A$ is Schur stable and that $G_{d}(1)$ has full column
rank.

Recall that the idea behind the estimator (7) is to obtain a prediction
$\hat{y}(t)$ for $y(t)$ using (7a)–(7b), and then to use $\hat{y}(t)$ and the
current measurement $y(t)$ to update $\hat{d}(t)$ using (7c). Translating this
to a data-driven framework, our approach will be to use (9) in place of
(7a)–(7b) to generate the estimate $\hat{y}(t)$, and then to verbatim use
(7c). Suppose that we have collected historical data as in (8), and that the
inputs $\mathrm{col}(u_{\rm d},d_{\rm d})$ are persistently exciting of order
$T_{\rm ini}+1+n(\mathscr{B})$ with $T_{\rm ini}\geq\ell(\mathscr{B})$. At
time step $t$, we construct the vectors of recent past data

$\displaystyle\hat{y}_{\rm p}$ $\displaystyle=\mathrm{col}(\hat{y}(t-T_{\rm
ini}),\ldots,\hat{y}(t-1))$ $\displaystyle u_{\rm p}$
$\displaystyle=\mathrm{col}(u(t-T_{\rm ini}),\ldots,u(t-1))$
$\displaystyle\hat{d}_{\rm p}$ $\displaystyle=\mathrm{col}(\hat{d}(t-T_{\rm
ini}),\ldots,\hat{d}(t-1)).$

Following (9), we formulate and solve the system of equations

$\begin{bmatrix}U_{\rm p}\\\ D_{\rm p}\\\ Y_{\rm p}\\\ U_{\rm f}\\\ D_{\rm
f}\\\ Y_{\rm f}\end{bmatrix}g=\begin{bmatrix}u_{\rm p}\\\ \hat{d}_{\rm p}\\\
\hat{y}_{\rm p}\\\ u(t)\\\ \hat{d}(t)\\\
\hat{y}(t)\end{bmatrix}\quad\Rightarrow\quad\hat{y}(t)=\underbrace{Y_{\rm
f}\begin{bmatrix}U_{\rm p}\\\ D_{\rm p}\\\ Y_{\rm p}\\\ U_{\rm f}\\\ D_{\rm
f}\end{bmatrix}^{\dagger}}_{\coloneqq\mathcal{M}}\begin{bmatrix}u_{\rm p}\\\
\hat{d}_{\rm p}\\\ \hat{y}_{\rm p}\\\ u(t)\\\ \hat{d}(t)\\\ \end{bmatrix}.$
(10)

Note that (10) allows us to compute an estimate $\hat{y}(t)$ for $y(t)$ using
only recorded historical data, recently computed online variables $u_{\rm p}$,
$\hat{d}_{\rm p}$, $\hat{y}_{\rm p}$, and the current values $u(t)$, and
$\hat{d}(t)$ of the input and disturbance estimates. In fact, we do not even
require $\hat{d}(t)$ to evaluate the above, since due to strict causality of
(2) on the $d\mapsto y$ channel, the final column of the matrix $\mathcal{M}$
is identically zero. Our data-driven disturbance estimator is therefore

$\displaystyle\hat{y}(t)$ $\displaystyle=\mathcal{M}\mathrm{col}(u_{\rm
p},\hat{d}_{\rm p},\hat{y}_{\rm p},u(t),\mathbbold{0}_{q})$ (11a)
$\displaystyle\hat{d}(t)$ $\displaystyle=\hat{d}(t-1)-\varepsilon
L(\hat{y}(t)-y(t)).$ (11b)

In the ideal case where measurements are collected without noise, the results
produced by the data-driven estimator (11) are identical to the results
produced by the model-based estimator (7), since (11a) is a data-based
representation of (7a)–(7b). Therefore (11) is a data-based substitute for
(7).

The final issue to address concerns the tuning of the estimator gain $L$ in
(11). Proposition II.1 provides the tuning suggestion $L=G_{d}(1)^{\dagger}$.
While $G_{d}(1)$ could be obtained empirically from repeated step response
experiments, it can also be obtained directly from the exact same historical
data used to construct the matrix $\mathcal{M}$ in (11). The following result
is an adaptation of [28, Thm. 4.1].

######  Lemma III.1 (DC Gain From Trajectory Data)

Consider the previously defined historical data $(u_{\rm d},y_{\rm d},d_{\rm
d})$ and define

$\displaystyle y_{\rm d}^{\rm diff}$ $\displaystyle=(y_{\rm d}(2)-y_{\rm
d}(1),\ldots,y_{\rm d}(T)-y_{\rm d}(T-1))\in(^{p})^{T-1}$ $\displaystyle
u_{\rm d}^{\rm diff}$ $\displaystyle=(u_{\rm d}(2)-u_{\rm d}(1),\ldots,u_{\rm
d}(T)-u_{\rm d}(T-1))\in(^{m})^{T-1}$

with associated Hankel matrices $Y^{\rm
diff}=\mathscr{H}_{\ell(\mathscr{B})}(y_{\rm d}^{\rm diff})$ and $U^{\rm
diff}=\mathscr{H}_{\ell(\mathscr{B})}(u_{\rm d}^{\rm diff})$. Then

$G_{d}(1)=Y_{\rm f}\begin{bmatrix}Y^{\rm diff}\\\ U^{\rm diff}\\\ U_{\rm p}\\\
D_{\rm p}\end{bmatrix}^{\dagger}\begin{bmatrix}0\\\ 0\\\ 0\\\
I_{q}\end{bmatrix}.$

Combining Lemma III.1 with the tuning $L=G_{d}(1)^{\dagger}$, the disturbance
estimator (11) provides a _completely model-free_ solution to the asymptotic
disturbance estimation problem; the only required tuning is the single scalar
parameter $\varepsilon\in(0,1)$.

##  IV Data-Driven Fast Frequency Control using Inverter-Based Resources

We now detail the application of our data-driven disturbance estimation
methods to fast frequency control using IBRs. First, in Section IV-A we
describe how the theory outlined in Section II and Section III is adapted for
application to the area-based fast frequency control problem described in
Section I. Section IV-B and IV-C present the nonlinear three-area test system,
along with simulation results validating our approach on several scenarios
including load increases, generation trips, and three-phase faults.

###  IV-A Application of the Theory for Multi-Area Frequency Control

Consider a large interconnected power system, partitioned into many small
LCAs, each of which contains IBRs which are dispatchable within specified
limits. The model (2) will be taken as describing the linearized frequency
dynamics of an LCA around a dispatch point of interest. This model is
_unknown_ , and may have an arbitrarily high number of states, describing,
e.g., the electromechanical dynamics and control systems of synchronous
generators, wind turbines, load dynamics, and so forth.

Due to the small spatial scale of each LCA, the measured frequency is roughly
uniform within the LCA, even during transient conditions (modulo, e.g., intra-
area modes). Analogously, the effect of a power imbalance within the LCA on
the frequency is approximately independent of the specific location of the
imbalance within the LCA. Thus, to a good approximation, power imbalance
influences the frequency of each LCA in a lumped fashion. As a result of these
observations, for the model (2) we make the following selections for inputs
and outputs. The measurement $y(t)=\Delta f(t)\in$ will be a single local
measurement of frequency deviation222An average or weighted average of
frequency measurements from across the LCA may also be used.; this further
implies that $D=0$. The disturbance $d(t)=\Delta P_{\rm u}\in$ will model
aggregate _unmeasured_ generation-load imbalance within the LCA; this is the
disturbance we wish to quickly compensate through real-time redispatch of
IBRs. The control input $u(t)=\Delta P_{\mathrm{ibr,tot}}(t)$ will model the
sum of all IBR power set-points, again relative to scheduled dispatch values.
As disturbances and generation are effectively lumped, we further assume that
$B=B_{d}$, i.e., disturbance and control signals enter through the same
channel. All power flows to neighboring LCAs are assumed to be measured, and
thus the total tie flow $\Delta P_{\mathrm{tie}}(t)$ out of the area is
considered as a _measurable_ disturbance, and is lumped with the control
signal $u(t)$.

The procedure used for collecting historical data will be described in the
next section. Due to the assumptions that $B=B_{d}$ and $D=0$, in developing
the data-driven estimator of Section III, it suffices to collect historical
data $(u_{\rm d},y_{\rm d})$ in the absence of significant unmeasured
disturbances, i.e., during normal system conditions. In this case, the
estimator (11a) becomes

$\displaystyle\Delta\hat{f}(t)$ $\displaystyle=Y_{\rm f}\begin{bmatrix}U_{\rm
p}\\\ Y_{\rm p}\\\ U_{\rm f}\end{bmatrix}^{\dagger}\begin{bmatrix}\Delta
P_{\rm ibr,tot,p}-\Delta P_{\rm tie,p}+\Delta\hat{P}_{\rm u,p}\\\
\Delta\hat{f}_{\rm p}\\\ 0\\\ \end{bmatrix}$ (12a)
$\displaystyle\Delta\hat{P}_{\rm u}(t)$ $\displaystyle=\Delta\hat{P}_{\rm
u}(t-1)-\varepsilon\frac{1}{G_{d}(1)}(\Delta\hat{f}(t)-\Delta f(t)).$ (12b)

and the gain computation of Lemma III.1 reduces to

$G_{d}(1)=Y_{\rm f}\begin{bmatrix}Y^{\rm diff}\\\ U_{\rm p}\\\
\end{bmatrix}^{\dagger}\begin{bmatrix}0\\\ 1\end{bmatrix}\in.$

The algorithm (12) estimates the load unmeasured imbalance. To compensate the
imbalance, the total change in IBR power set-points required within the LCA is
updated as

$\Delta P_{\rm ibr,tot}(t+1)=\Delta\hat{P}_{\rm u}(t).$

Finally, to compute the set-points for individual IBRs within each LCA, the
total set-point change $\Delta P_{\rm ibr,tot}(t)$ is allocated amongst the
IBRs. This can be done via the optimal constrained allocation method presented
in [3], or alternatively, the total can be allocated to individual IBRs based
on participation factors and then saturated to respect the device power
limits; the former approach is used here.

###  IV-B Description of Test System and Data Collection

The three-LCA test system under consideration is shown in Figure 3, where each
of the individual areas is modified based on the IEEE 3-machine 9-bus system
[29]. In the modified model, two dispatchable IBRs are added to each LCA to
facilitate fast frequency control. For primary frequency support, each IBR
unit has a 5% droop curve on its respective base power. Doubly-fed induction
generator (DFIG) wind turbine systems replace two traditional synchronous
generators (SGs) in Area 1, and replace one SG each in Area 2 and Area 3. One
static var compensator (SVC) has been added to both Area 1 and Area 2, while
two SGs in Area 3 are replaced with two IBRs having the same ratings. The pre-
disturbance generation data and unit capacity limits for the system can be
found in Table I.

TABLE I: Generator and IBR Data.

Node | Gen. ID | Rating (MVA) | Dispatch (MW)
---|---|---|---
1 | WT1 | 142.2 | 72.24
1 | IBR2 | 50 | 15
3 | G1 | 192 | 126
3 | IBR1 | 50 | 25
5, 11, 17 | WT2, WT3, WT4 | 142.2 | 85
7 | G2 | 247.5 | 71.99
7 | IBR4 | 50 | 20
9 | G3 | 192 | 133
11 | IBR3 | 50 | 10
13 | G4 | 247.50 | 72.24
13 | IBR6 | 50 | 5
15 | G5 | 192 | 128
17 | IBR5 | 50 | 30

The next step in implementing the data-driven disturbance estimator (12) is
the one-time collection of measurements from the system in response to a
sufficiently rich input. To this end, we assume that the IBRs will be used to
excite the system for the data collection phase. As described in Section III,
the input signal should be persistently exciting of sufficient high order. For
the purposes of these tests, the set-point change

$\Delta P_{\rm ibr}(t)=\sin(12\pi t)+w(t)$ (13)

in units of MW was provided to each IBR in each LCA. The signal, plotted in
Figure 2, consists of a sinusoidal perturbation of 1 MW, with band-limited
white noise $w(t)$ with power spectral density (PSD) of 0.1. Note that the
amplitude of the perturbing IBR signal is relatively small compared to the
overall generation/demand in the system ($\approx 800$ MW).333Further
investigation into the design of practical input signals for data collection
is deferred to future work, but see [30, 31] for recent theoretical results.

Figure 2: Persistently exciting IBR set-point change for data collection
phase.

For each LCA, measurements of frequency deviation and net inter-area power
flow were collected for 10 seconds with a sampling period of 0.1s, leading to
a total of $T=101$ total historical data points for each LCA. The length of
recent past data used in (12) was $T_{\rm ini}=7$, which is a minimal value
obtained by starting from a large value and reducing until performance
degradation was first observed in simulation. The controller gain
$\varepsilon$ in (12) was set via tuning at $\varepsilon=0.1$.

Finally, to improve the numerical conditioning of the matrix inversion step in
(12), the matrix $\mathrm{col}(U_{\rm p},Y_{\rm p},U_{\rm f})$ is replaced by
a low-rank approximation thereof, obtained by computing the singular value
decomposition and retaining only the first three dominant singular values and
vectors [32]. The interpretation of this regularization step is that retaining
only a small number of singular values takes into account only the most
dominant sub-behaviours, and therefore removes the effects of unimportant fast
dynamics embedded in the data [19].

###  IV-C Case Studies

We consider five disturbance scenarios to illustrate the effectiveness of our
data-driven fast frequency controller design. The scenarios are:

  1. (i)

two large step load changes of different sizes in one area, the first being
small enough to compensate with local IBR resources, while the second requires
active power support from IBRs in adjacent areas;

  2. (ii)

two simultaneous disturbances, consisting of a step load change and a rate-
limited change in the wind power system; the test system in this case is
modified to have reduced inertia and increased wind turbine penetration;

  3. (iii)

the loss of a generator;

  4. (iv)

a symmetric three phase-to-ground fault;

  5. (v)

a step load change in a larger five-area test system.

In all scenarios, we consider the hierarchical fast frequency control
architecture proposed in [3], and compare the model-based disturbance
estimator of that work against the data-driven disturbance estimator presented
in this paper. a comparison of the model-based scheme and an AGC-type scheme
can be found in [3], where it was shown that our model-based redispatch can
significantly outperform a more traditional AGC-type implementation, and thus
we compare only against the model-based scheme. The performance is also
compared against a baseline case without any supplementary control scheme,
wherein frequency support is provided only through primary droop control
action of both generators and IBRs. All simulations are performed with
measurement and control signal delays of 300ms, representing worst-case wide-
area communication delays, and with white noise of standard deviation
$10^{-6}$ (resp. $2\times 10^{-2}$) added to the frequency (resp. inter-area
power flow) p.u. measurements.

Figure 3: Cyber-physical system illustrating frequency control approach [3].

#### Scenario #1: Step Load Increase

This scenario illustrates the performance of the control scheme during step
changes in loading. In the first case, a large load disturbance of 60 MW is
introduced at $t=2s$ at bus 8 in Area 2. The frequency response and
disturbance estimate of the system is shown in Figure 4. Both the data-driven
and model-based disturbance estimators localize the disturbance to Area 2, and
quickly redispatch the IBRs in Area 2 (Figure 5) to compensate, with minimal
transient response from the other areas. The frequency is quickly restored to
the nominal value, with minimal oscillation. Similarly, the measured net tie-
line deviation and IBR output power plots in Figure 5 show a fast, non-
oscillatory response. Due to the area-wise decentralized nature of the control
scheme, the disturbance estimate, tie-line deviation and IBR outputs in the
non-contingent areas return to their pre-disturbance values in steady-state.
The black dotted lines overlaid on the responses in Figure 4 show the response
when $\varepsilon$ is decreased by a factor of 10 to $\varepsilon=0.01$.
Tuning of the single scalar parameter $\varepsilon$ therefore allows for
adjustment of the transient response.

Figure 4: Frequency and disturbance estimate during a 60 MW load change at bus
8 in Area 2. Figure 5: Tie-line deviation and active power profiles during a
60 MW load change; dashed lines in the lower plots indicate the responses
under model-based estimation.

Figures 6 and 7 show the closed-loop response under a more severe load
disturbance of 130 MW, again at bus 8 in Area 2. The disturbance exceeds the
available spare IBR resources in the contingent area, and the inter-area
coordination scheme proposed in [3] is activated to facilitate support from
IBRs in neighboring areas. While the details of these coordination scheme are
not germane to our discussion here, the plots illustrate that the purely data-
driven approach presented here produces similar results to a model-based
approach. In summary, for both large and small load disturbances, our data-
driven estimator is able to quickly spatially localize and compensate for a
load disturbance, using absolutely no model information.

Figure 6: Frequency and disturbance estimate during a 130 MW load change at
bus 8 in Area 2. Figure 7: Tie-line deviation and active power profiles during
a 130 MW load change; dashed lines in the lower plots indicate the responses
under model-based estimation.

#### Scenario #2: Heavy Renewable Penetration

This scenario assesses the performance of the controller when one area of the
power system is fed entirely by renewable energy sources. For this scenario,
we have modified area three of the test system shown in Figure 3: in the
modified system, generators G4 and G5 are replaced by a non-dispatchable IBRs
with similar capacities and pre-disturbance dispatch points. A synchronous
condenser is connected at bus 15 in area 3 to regulate the voltage by
providing reactive power support.

For this scenario, a step disturbance of 40 MW is introduced at bus 14 in Area
3 at $t=2$s. Simultaneously, a rate limited change in the wind speed to the
wind turbine power system connected at bus 17 occurs. This is done to simulate
the net-load variability and uncertainty in power supply introduced by
renewables [26]. From the responses shown in Figures 8, 9, we observe that the
IBRs in Area 3 respond quickly to the disturbance, and additionally, slowly
ramp up their injections in response to the rate-limited decrease in the wind
power. Throughout this process, the frequency is maintained very close to the
nominal value.

From this test, we conclude that the data-driven controller functions well in
a low-inertia area, and even produces a closed-loop frequency response which
is slightly improved compared to the crude model-based design from [3].

Figure 8: Frequency response and disturbance estimates during combined step
load change and a rate-limited change in wind speed. Figure 9: Tie-line
deviation and active power profiles during combined step load change and a
rate-limited change in wind speed; dashed lines in the lower plots indicate
the responses under model-based estimation.

#### Scenario #3: Generator Trip

This scenario examines the performance of the controller during the loss of
generator G2 in Area 2 at $t=2$s. The response of the system to this loss is
plotted in Figure 10. Similar to Scenario #2 above, we can observe that while
the response of the controller under both data-driven and model-based
disturbance estimation is quite fast and satisfactory, the data-driven
controller also outperforms the model-based controller. This scenario
illustrates the robustness of the method, as the collected data used for the
design of the LCA controller was collected on the system _including_ the
inertia and primary response of generator G2. Despite the model mismatch
arising from the loss of the generator, the overall control response is
similar to that in Scenario #1.

Figure 10: Frequency and disturbance estimate during a loss of generator G2

#### Scenario #4: Three-Phase Fault

The scenario assesses the performance of the controller during a symmetrical
three-phase line-to-ground fault, which was introduced at bus 10 in Area 2 at
$t=2$s and cleared after $0.1$s. The response of the system is shown in
Figures 11 and 12. It can be observed that the frequency response of the
system with and without the controllers is very similar, with a small
transient response in the disturbance estimate computed by the estimators.
This indicates that the controller is able to effectively detect and respond
to frequency events. The performance of the data-driven estimator in this
scenario is satisfactory and is similar to the model-based estimator.

Figure 11: Frequency and disturbance estimate during a three-phase fault in
Area 2. Figure 12: Tie-line deviation and active power profiles during a
three-phrase fault in Area 2; dashed lines in the lower plots indicate the
responses under model-based estimation.

#### Scenario #5: Step loading change in larger test system

This scenario is introduced to test the performance of the data-driven
approach on the larger 5-LCA system shown in Figure 13, which is an IEEE
benchmark model with 68 buses and is extensively described in [33].

Following the same data collection approach from Section IV-B, the perturbing
set-point change for the IBRs in this system was as in (13), with the noise
term having a power spectral density of 1. For each LCA, measurements of
frequency deviation and net inter-area power flow were collected for 10
seconds with a sampling period of 25 ms. The selections of $T_{\rm ini}$ and
$\varepsilon$ were tuned for each area individually; in particular, it was
found that larger values are beneficial for performance in larger areas. For
instance, in the large NETS area the best tunings found were $T_{\rm ini}=119$
and $\varepsilon=0.3$, while in Area 4 the smaller values $T_{\rm ini}=7$ and
$\varepsilon=0.01$ were found to be more appropriate. The intuition behind
this is that larger areas are more dynamically complex, and require more data
for initialization within our data-driven approach.

Figure 13:  IEEE 68-bus test system [3].

For this case, we introduce a large step load change of 450 MW at bus 33 in
the NYPS area at $t=2$s. The frequency response of the system to this
disturbance is plotted in Figure 14, while Figure 15 shows the estimate of the
load disturbance with the data-driven and model-based approaches. Similar to
Scenarios #2 and #3 above, we can observe that while the response of the
controller under both data-driven and model-based disturbance estimation is
quite fast and satisfactory, the data-driven controller slightly outperforms
the model-based controller. This is because the data-driven controller
potentially captures more of the dynamics of the areas than the simplified
model used for the model-based control approach.

Figure 14: Frequency response to a load change of 450 MW at bus 33. Figure 15:
Disturbance estimate during a load change of 450 MW at bus 33.

##  V Conclusion

This paper has detailed the development of a data-driven area-based fast
frequency control scheme, which extends the previous work in [3] to be
completely model-free. The controller rapidly redispatches inverter-based
resources to compensate for local power imbalances within the bulk power
system. The design is based only on recorded historical data, requires no
model building or identification, and is tuned through adjustment of a single
scalar parameter, making it potentially attractive for practical
implementation. Simulation results on a nonlinear test system have been
presented to validate the approach. Future work will be primarily concerned
with minimizing the disruption to the system during the data collection phase
via the design of improved excitation signals for IBRs, along with extension
of the approach here to moving-horizon disturbance estimation.

## References

  * [1] IRENA, “Global energy transformation: A roadmap to 2050,” 2018.
  * [2] B. K. Poolla, S. Bolognani, and F. Dörfler, “Optimal placement of virtual inertia in power grids,” IEEE Trans. Autom. Control, vol. 62, no. 12, pp. 6209–6220, 2017.
  * [3] E. Ekomwenrenren, Z. Tang, J. W. Simpson-Porco, E. Farantatos, M. Patel, and H. Hooshyar, “Hierarchical coordinated fast frequency control using inverter-based resources,” IEEE Trans. Power Syst., vol. 36, no. 6, pp. 4992–5005, 2021.
  * [4] P. M. Anderson and M. Mirheydar, “A low-order system frequency response model,” IEEE Trans. Power Syst., vol. 5, no. 3, pp. 720–729, 1990.
  * [5] A. M. Prostejovsky, M. Marinelli, M. Rezkalla, M. H. Syed, and E. Guillo-Sansano, “Tuningless load frequency control through active engagement of distributed resources,” IEEE Transactions on Power Systems, vol. 33, no. 3, pp. 2929–2939, 2017.
  * [6] H. Bevrani, “Robust power system frequency control,” 2014.
  * [7] A. N. Venkat, I. A. Hiskens, J. B. Rawlings, and S. J. Wright, “Distributed mpc strategies with application to power system automatic generation control,” IEEE Trans. Control Syst. Tech., vol. 16, no. 6, pp. 1192–1206, 2008.
  * [8] A. M. Ersdal, L. Imsland, and K. Uhlen, “Model predictive load-frequency control,” IEEE Trans. Power Syst., vol. 31, no. 1, pp. 777–785, 2015.
  * [9] J. Köhler, M. A. Müller, N. Li, and F. Allgöwer, “Real time economic dispatch for power networks: A distributed economic model predictive control approach,” in Proc. IEEE CDC, (Melbourne, Australia), pp. 6340–6345, 2017.
  * [10] P. R. B. Monasterios and P. Trodden, “Low-complexity distributed predictive automatic generation control with guaranteed properties,” IEEE Trans. Smart Grid, vol. 8, no. 6, pp. 3045–3054, 2017.
  * [11] C. Chen, M. Cui, F. Li, S. Yin, and X. Wang, “Model-free emergency frequency control based on reinforcement learning,” IEEE Transactions on Industrial Informatics, vol. 17, no. 4, pp. 2336–2346, 2020.
  * [12] T. Yu, B. Zhou, K. W. Chan, L. Chen, and B. Yang, “Stochastic optimal relaxed automatic generation control in non-markov environment based on multi-step $q(\lambda)$ learning,” IEEE Transactions on Power Systems, vol. 26, no. 3, pp. 1272–1282, 2011.
  * [13] Z. Yan and Y. Xu, “Data-driven load frequency control for stochastic power systems: A deep reinforcement learning method with continuous action search,” IEEE Transactions on Power Systems, vol. 34, no. 2, pp. 1653–1656, 2018.
  * [14] Z. Yan and Y. Xu, “A multi-agent deep reinforcement learning method for cooperative load frequency control of a multi-area power system,” IEEE Transactions on Power Systems, vol. 35, no. 6, pp. 4599–4608, 2020.
  * [15] B. Recht, “A tour of reinforcement learning: The view from continuous control,” Annual Review of Control, Robotics, and Autonomous Systems, vol. 2, pp. 253–279, 2019.
  * [16] W. Liu, G. Geng, Q. Jiang, H. Fan, and J. Yu, “Model-free fast frequency control support with energy storage system,” IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 3078–3086, 2019.
  * [17] P. Hidalgo-Gonzalez, R. Henriquez-Auba, D. S. Callaway, and C. J. Tomlin, “Frequency regulation using data-driven controllers in power grids with variable inertia due to renewable energy,” in 2019 IEEE Power & Energy Society General Meeting (PESGM), pp. 1–5, IEEE, 2019.
  * [18] I. Markovsky and P. Rapisarda, “Data-driven simulation and control,” International Journal of Control, vol. 81, no. 12, pp. 1946–1959, 2008.
  * [19] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the deepc,” in 2019 18th European Control Conference (ECC), pp. 307–312, IEEE, 2019.
  * [20] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” IEEE Transactions on Automatic Control, vol. 65, no. 3, pp. 909–924, 2019.
  * [21] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Data-driven model predictive control with stability and robustness guarantees,” IEEE Transactions on Automatic Control, vol. 66, no. 4, pp. 1702–1717, 2020.
  * [22] H. J. Van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: a new perspective on data-driven analysis and control,” IEEE Transactions on Automatic Control, vol. 65, no. 11, pp. 4753–4768, 2020\.
  * [23] Y. Zhao, T. Liu, and D. J. Hill, “A data-enabled predictive control method for frequency regulation of power systems,” in 2021 IEEE PES Innovative Smart Grid Technologies Europe (ISGT Europe), pp. 01–06, IEEE, 2021.
  * [24] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. M. De Moor, “A note on persistency of excitation,” IFAC Syst & Control L, vol. 54, no. 4, pp. 325–329, 2005.
  * [25] J. C. Willems and J. W. Polderman, Introduction to mathematical systems theory: a behavioral approach, vol. 26.  Springer Science & Business Media, 1997.
  * [26] A. Ulbig, T. S. Borsche, and G. Andersson, “Impact of Low Rotational Inertia on Power System Stability and Operation,” in Proc. IFAC World C, vol. 47, (Cape Town, South Africa), pp. 7290–7297, Aug. 2014.
  * [27] W. S. Levine, The Control Systems Handbook: Control System Advanced Methods.  CRC press, 2018.
  * [28] G. Bianchin, M. Vaquero, J. Cortés, and E. Dall’Anese, “Data-driven synthesis of optimization-based controllers for regulation of unknown linear systems,” 2021.  https://arxiv.org/abs/2103.16067.
  * [29] A. Delavari, I. Kamwa, and P. Brunelle, “Simscape power systems benchmarks for education and research in power grid dynamics and control,” in IEEE Canadian Conf. on ECE, pp. 1–5, 2018.
  * [30] C. D. Persis and P. Tesi, “Designing experiments for data-driven control of nonlinear systems,” 2021.  https://arxiv.org/abs/2103.16509.
  * [31] H. J. van Waarde, “Beyond persistent excitation: Online experiment design for data-driven modeling and control,” IEEE Control Syst. Let., vol. 6, pp. 319–324, 2022.
  * [32] G. H. Golub and C. F. Van Loan, Matrix computations.  JHU press, 2013.
  * [33] B. Pal and B. Chaudhuri, Robust control in power systems.  Springer Science & Business Media, 2006.


# Wannier Functions Dually Localized in Space and Energy

Aaron Mahler<EMAIL_ADDRESS>Duke University, Department of Physics,
Durham, NC 27708  Jacob Z. Williams  Duke University, Department of Chemistry,
Durham, NC 27708  Neil Qiang Su  Department of Chemistry, Key Laboratory of
Advanced Energy Materials Chemistry (Ministry of Education) and Renewable
Energy Conversion and Storage Center (RECAST), Nankai University, Tianjin
300071, China  Duke University, Department of Chemistry, Durham, NC 27708
Weitao Yang  Duke University, Department of Chemistry, Durham, NC 27708  Duke
University, Department of Physics, Durham, NC 27708

###### Abstract

Wannier functions are real-space representations of Bloch orbitals that
provide a useful picture for chemical bonding and offer a localized
description of single-particle wave functions. There is a unitary freedom in
the construction of Wannier functions from Bloch orbitals, which can be chosen
to produce Wannier functions that have advantageous properties. A popular
choice for this freedom is the minimization of spatial variance, which leads
to maximally localized Wannier functions. We minimize a weighted sum of the
spatial and energy variances, yielding what we call dually localized Wannier
functions. Localization in energy results in Wannier functions that are
associated with a particular energy and allows for dually localized Wannier
functions to be comprised of occupied and unoccupied spaces simulatneously.
This can lead to Wannier functions that are fractionally occupied, which is a
key feature that is used in correcting the delocalization error in density
functional approximations. We show how this type of localization results in
bonding (antibonding) functions for the occupied (unoccupied) spaces around
the frontier energy of silicon in the diamond lattice and of molecular
ethylene. Dually localized Wannier functions therefore offer a relevant
description for chemical bonding and are well suited to orbital-dependent
methods that associate Wannier functions with specific energy ranges without
the need to consider the occupied and unoccupied spaces separately.

##  I Introduction

In the single-particle picture of electronic structure theory, the
eigenfunctions of the Hamiltonian may be delocalized in space. This is
especially true when the Hamiltonian is periodic and the solutions are Bloch
functions
$\psi_{n}^{\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot\mathbf{r}}u_{n}^{\mathbf{k}}(\mathbf{r})$,
where $u_{n}^{\mathbf{k}}(\mathbf{r})$ shares the periodicity of the
Hamiltonian. Bloch functions are indexed by $\mathbf{k}$, a point in
reciprocal space, and a band index corresponding to the $n^{\text{th}}$
eigenvalue for that $\mathbf{k}$-point, i.e.,
$H\psi_{n}^{\mathbf{k}}=\epsilon_{n}^{\mathbf{k}}\psi_{n}^{\mathbf{k}}$. As
can be seen from the plane wave $e^{i\mathbf{k}\cdot\mathbf{r}}$ in their
definition, the Bloch orbitals are delocalized in space while being perfectly
localized in energy, limiting their utility in describing spatially local
properties such as chemical bonds. Dual to the Bloch orbitals are the Wannier
functions, which are the Fourier transform of the Bloch orbitals over the
Brillouin zone [1]. They are localized in real space, making them useful for
evaluating position-dependent quantities such as the dipole moment [2, 3], as
a basis set for large-scale simulations [4, 5, 6, 7], interpolating band
structures [8], as a picture of chemical bonding [9, 10], and as localized
orbitals in orbital-dependent methods [11, 12, 13, 14].

If a finite number $N_{k}$ of $\mathbf{k}$-points in the first Brillouin zone
are sampled, Wannier functions take the form

$\lvert
w_{n}^{\mathbf{R}}\rangle=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\lvert\psi_{n}^{\mathbf{k}}\rangle,$
(1)

where $\mathbf{R}$ indexes unit cells in the unfolded supercell on which the
Wannier functions are periodic. Equal weighting of all $\mathbf{k}$-points is
equivalent to sampling on a Monkhorst-Pack mesh [15], which we assume
throughout the text. This $\mathbf{k}$-sampling procedure also yields Wannier
functions that are periodic on an unfolded supercell obeying Born–von Karman
boundary conditions, with the number of primitive unit cells in each real-
space lattice direction equal to the number of $\mathbf{k}$-points sampled in
the corresponding reciprocal lattice direction [16]. The Fourier transform
normalization convention in Eq. (1) also implies a Bloch orbital normalization
$\langle\psi_{n}^{\mathbf{k}}|\psi_{m}^{\mathbf{q}}\rangle=N_{k}\delta_{nm}\delta_{\mathbf{k}\mathbf{q}}$.
The inverse transform from Bloch orbital to Wannier function is then

$\lvert\psi_{n}^{\mathbf{k}}\rangle=\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}\lvert
w_{n}^{\mathbf{R}}\rangle.$ (2)

Since the Bloch orbitals are eigenfunctions, they have a phase freedom in
their formulation. This leads to a gauge freedom in the construction of the
Wannier functions:

$\begin{split}\lvert
w_{n}^{\mathbf{R}}\rangle&=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}e^{i\theta_{n}(\mathbf{k})}\lvert\psi_{n}^{\mathbf{k}}\rangle.\end{split}$
(3)

For obtaining localized Wannier functions, a natural choice for this gauge
freedom is to choose the phase factors $\theta_{n}(\mathbf{k})$ such that
$\lvert\psi_{n}^{\mathbf{k}}\rangle$ is as smooth as possible in its
dependence on $\mathbf{k}$. If it is analytic, the resulting Wannier function
is exponentially localized [17]. However, multiple bands may cross one another
in the Brillouin zone, yielding $\mathbf{k}$-points with a set of degenerate
Bloch functions. Any unitary combination $U^{\mathbf{k}}$ of the degenerate
bands remains a set of eigenfunctions of the Hamiltonian. These freedoms lead
to _generalized_ Wannier functions, which allows the bands to mix at a fixed
$\mathbf{k}$-point by an arbitrary unitary operator $U^{\mathbf{k}}$ [18],

$\begin{split}\lvert
w_{n}^{\mathbf{R}}\rangle&=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\sum_{m}U_{mn}^{\mathbf{k}}\lvert\psi_{m}^{\mathbf{k}}\rangle\\\
&=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\lvert\phi_{n}^{\mathbf{k}}\rangle.\end{split}$
(4)

We call the set $\\{\lvert\phi_{n}^{\mathbf{k}}\rangle\\}$ the _transformed_
Bloch orbitals. Thus, a generalized Wannier function can be comprised of a
combination of bands, and the combination may vary across the Brillouin zone.

###  I.1 Cost functions

The unitary freedom in the construction of the Wannier functions can be chosen
to satisfy a desired criterion. In the case of the maximally localized Wannier
functions (MLWFs) proposed by Marzari and Vanderbilt [9], the $U^{\mathbf{k}}$
minimize the spatial variance

$\begin{split}\Omega&=\sum_{n}\langle w_{n}|\Delta r^{2}|w_{n}\rangle,\\\
\end{split}$ (5)

where $\langle
f|x|g\rangle=\int_{\mathcal{D}}d\mathbf{r}\,f^{*}(\mathbf{r})x(\mathbf{r})g(\mathbf{r})$
for any operator $x$ and functions $f,g$. The domain of integration
$\mathcal{D}$ depends on the periodicity of $f$ and $g$. In the case of the
Wannier functions, $\mathcal{D}$ is the Born–von Karman unit cell [16]. We
write $\langle\Delta x^{2}\rangle=\langle
x^{2}\rangle-\lvert\langle{x}\rangle\rvert^{2}$ for the variance of $x$. Note
that the MLWF localization cost function of Eq. (5) is the bulk-system
extension of the Foster-Boys localization for finite systems [19]. The Wannier
functions are translationally symmetric on the primitive unit cell, so without
loss of generality we only include the Wannier functions indexed by the home
unit cell, $\mathbf{R}=\mathbf{0}$. Summing over all unit cells $\mathbf{R}$
in the Born–von Karman unit cell would only change the cost function by a
multiplicative constant. The MLWF formulation was originally applied to a
composite set of energy bands, which is a set of Bloch orbitals separated by
an energy gap from all other states at every point in the Brillouin zone. The
gradient of $\Omega$ was shown analytically for a composite set by Marzari and
Vanderbilt [9], and descent methods are applied to obtain the $U^{\mathbf{k}}$
which minimize $\Omega$. Applied to the occupied orbitals of gapped systems
with vanishing Chern numbers, MLWFs have been proven to be real and
exponentially localized in three and fewer dimensions [20, 17, 21, 22].

We now introduce the dual localization criteria, adding an energy variance
term to the MLWF cost function:

$\begin{split}F&=(1-\gamma)\sum_{n}\langle w_{n}|\Delta
r^{2}|w_{n}\rangle+\gamma\sum_{n}\langle w_{n}|\Delta h^{2}|w_{n}\rangle\\\
&=(1-\gamma)\Omega+\gamma\Xi,\end{split}$ (6)

where $h$ is the single-particle Hamiltonian, $\Omega$ denotes the spatial
spread, and $\Xi$ denotes the energy spread. The mixing term $\gamma$ can be
tuned to prioritize spatial or energy localization (see Sec. II.4). This cost
function was used in the localized orbital scaling correction (LOSC) [23],
which aims to systemically eliminate delocalization error using localized
orbitals [24]. In the case of finite systems, the localized orbitals based on
the optimization of the cost function, Eq. (6) are called orbitalets to
highlight its space and energy localization [24, 23]. The form of $F$ was
first suggested by Gygi et al. in [25]; however, the algorithm there proposed
is only applicable to $\Gamma$-sampled systems, where the Brillouin zone is
sampled only at the origin. In that case, the transformed Bloch orbitals are
equivalent to the Wannier functions, and varying the mixing term $\gamma$
progresses between MLWFs and Bloch orbitals. Additionally, the $\Gamma$-only
algorithm does not offer guidance on how to treat both occupied and unoccupied
spaces together. We show how using the cost function in Eq. (6) can produce
physically meaningful localizations when any number of virtual orbitals are
included in its domain. The resulting Wannier functions are localized in both
space and energy; we therefore refer to them as dually localized Wannier
functions (DLWFs).

###  I.2 Occupied and unoccupied spaces

The inclusion of energy localization allows occupied and unoccupied Bloch
orbitals to be considered together for Wannier function construction. When
using a cost function that solely minimizes spatial variance, only a set of
Bloch bands that are close together in energy are typically used to construct
MLWFs. This is because the basis of Hamiltonian eigenstates is complete and
the MLWF cost function allows Bloch bands to mix freely regardless of their
relative energy. Adding more bands to the construction of MLWFs will therefore
yield states arbitrarily localized in space, while losing virtually all
information about the energy content carried by the Bloch bands from which
they were constructed. On the other hand, if energy variance is included in
addition to spatial variance, then the mixing of Bloch orbitals that are far
apart in energy is suppressed. This means that a set of Bloch orbitals
spanning a large energy range can be considered and results in Wannier
functions that are localized in space as well as associated with a particular
energy.

When there is an energy gap between the occupied and unoccupied Bloch
orbitals, a set of Wannier functions can be comprised of the occupied Bloch
orbitals, and the projector onto the occupied space may be written in the
Bloch orbital or Wannier function basis as

$\begin{split}P_{\text{occ}}&=\frac{1}{N_{k}}\sum_{n\mathbf{k}}f_{n}^{\mathbf{k}}\lvert\psi_{n}^{\mathbf{k}}\rangle\\!\langle\psi_{n}^{\mathbf{k}}\rvert\\\
&=\sum_{n\mathbf{R}}f_{n}^{\mathbf{R}}\lvert
w_{n}^{\mathbf{R}}\rangle\\!\langle w_{n}^{\mathbf{R}}\rvert,\end{split}$ (7)

where $f_{n}^{\mathbf{k}},f_{n}^{\mathbf{R}}\in\\{0,1\\}$ are the occupations
of the Bloch orbitals and Wannier functions respectively. When the set of
Wannier functions includes both occupied and unoccupied bands, the Wannier
functions may not have integer occupancy: that is, $\langle
w_{n}^{\mathbf{R}}|P_{\text{occ}}|w_{n}^{\mathbf{R}}\rangle\not\in\\{0,1\\}$.
The fractional occupiations were used in the LOSC approach to correct the
major underestimation of the total energy in molecular ions at large bond
length, which is caused by the delocalization error in the functional
approximation. Fractional occupation also means the occupied projector
$P_{\text{occ}}$ is no longer diagonal in the transformed band index $n$. If
the transformed bands include all of the occupied Bloch orbitals, the occupied
projector can be written using the transformed Bloch orbitals as

$\begin{split}P_{\text{occ}}&=\frac{1}{N_{k}^{2}}\sum_{m\mathbf{q}}\sum_{n\mathbf{k}}\left(\langle\phi_{m}^{\mathbf{q}}|P_{\text{occ}}|\phi_{n}^{\mathbf{k}}\rangle\right)\lvert\phi_{m}^{\mathbf{q}}\rangle\\!\langle\phi_{n}^{\mathbf{k}}\rvert\\\
&=\frac{1}{N_{k}}\sum_{mn\mathbf{k}}\lambda_{mn}^{\mathbf{k}}\,\lvert\phi_{m}^{\mathbf{k}}\rangle\\!\langle\phi_{n}^{\mathbf{k}}\rvert,\end{split}$
(8)

where
$\lvert\phi_{n}^{\mathbf{k}}\rangle=\sum_{m}U_{mn}^{\mathbf{k}}\lvert\psi_{m}^{\mathbf{k}}\rangle$
is a transformed Bloch orbital. The $\lambda_{mn}^{\mathbf{k}}$ can be viewed
as local occupations of the transformed Bloch functions; observe that the
matrix $[\Lambda^{\mathbf{k}}]_{mn}=\lambda_{mn}^{\mathbf{k}}$ is not
necessarily diagonal in $n$ as it is when the Wannier functions are
constructed only from occupied Bloch bands. The trace of
$\Lambda^{\mathbf{k}}$ is $N_{f}^{\mathbf{k}}$, the number of Bloch states
below the Fermi energy at $\mathbf{k}$. For simplicity, we have omitted the
spin index, but these arguments also apply to open-shell systems, albeit for
each spin channel independently. As the unitary transform of a diagonal
projector, $\Lambda^{\mathbf{k}}$ is a Hermitian matrix with
$0\leq\lambda_{ii}^{\mathbf{k}}\leq 1$.

We can also write the occupied projector in terms of the Wannier functions as

$\begin{split}P_{\text{occ}}&=\sum_{m\mathbf{R}}\sum_{n\mathbf{T}}\left(\langle
w_{m}^{\mathbf{R}}|P_{\text{occ}}|w_{n}^{\mathbf{T}}\rangle\right)\lvert
w_{m}^{\mathbf{R}}\rangle\\!\langle w_{n}^{\mathbf{T}}\rvert\\\
&=\sum_{\mathbf{R}\mathbf{T}}\sum_{mn}\lambda_{mn}^{\mathbf{R}\mathbf{T}}\,\lvert
w_{m}^{\mathbf{R}}\rangle\\!\langle w_{n}^{\mathbf{T}}\rvert.\end{split}$ (9)

The matrix
$[\Lambda^{\mathbf{R}\mathbf{T}}]_{mn}=\lambda_{mn}^{\mathbf{R}\mathbf{T}}$ is
the element-wise discrete Fourier transform of the occupation matrix
$\Lambda^{\mathbf{k}}$ of transformed Bloch orbitals.
$\Lambda^{\mathbf{R}\mathbf{T}}$ can be viewed as the occupation matrix of the
Wannier functions; restricted to the home cell, it has unit cell trace
$\operatorname{Tr}_{c}[\Lambda^{(\mathbf{0}\mathbf{0})}]=\sum_{n}\Lambda^{(\mathbf{0}\mathbf{0})}_{nn}=N_{k}^{\text{-}1}\sum_{k}N_{f}^{\mathbf{k}}$.

##  II Methods

The procedure outlined by Marzari and Vanderbilt in [9] detailed how to find
the analytic gradient for minimizing $\Omega$ for a composite set of energy
bands. This was accomplished by splitting the cost function into two positive
definite quantities. One of these, $\Omega_{\text{I}}$, is invariant under the
Wannier functions’ gauge freedom, so minimizing $\Omega$ is equivalent to
minimizing only the gauge-dependent term $\widetilde{\Omega}$. For a set of
$N_{w}$ Wannier functions, the gauge-invariant term of the spatial cost is
given by

$\begin{split}\Omega_{\text{I}}&=\sum_{n}^{N_{w}}\Big{(}\langle
w_{n}^{\mathbf{0}}\rvert r^{2}\lvert
w_{n}^{\mathbf{0}}\rangle-\sum_{m}^{N_{w}}\sum_{\mathbf{R}}^{N_{k}}\left|\langle
w_{m}^{\mathbf{R}}\rvert\mathbf{r}\lvert
w_{n}^{\mathbf{0}}\rangle\right|^{2}\Big{)}.\end{split}$ (10)

The gauge-dependent term for the spatial cost is

$\begin{split}\widetilde{\Omega}&=\sum_{n\neq
m}^{N_{w}}\sum_{\mathbf{R}}^{N_{k}}\left|\langle
w_{m}^{\mathbf{R}}\rvert\mathbf{r}\lvert
w_{n}^{\mathbf{0}}\rangle\right|^{2}+\sum_{n}^{N_{w}}\sum_{\mathbf{R}\neq\mathbf{0}}^{N_{k}}\left|\langle
w_{n}^{\mathbf{R}}\rvert\mathbf{r}\lvert
w_{n}^{\mathbf{0}}\rangle\right|^{2}.\end{split}$ (11)

The valence bands of metals, and the conduction bands of most systems, cannot
be separated from the rest of the bands by an energy gap everywhere in the
Brillouin zone; such bands are said to form an entangled set. Souza et al.
developed a method called _disentanglement_ to extract a subset of interest
from a set of entangled bands [8]. Given $N_{b}$ Bloch bands, the
disentanglement procedure extracts $N_{w}\leq N_{b}$ composite bands, chosen
by minimizing $\Omega_{\text{I}}$. The disentangled bands represent the
smoothest possible $N_{w}$-dimensional subspace given the original $N_{b}$
bands, and are used to construct $N_{w}$ Wannier functions. To differentiate
between disentanglement and localization, we refer to the orbitals obtained by
minimizing $\Omega_{\text{I}}$ as disentangled Bloch orbitals, and to those
obtained by minimizing a cost function for a fixed number of (possibly
disentangled) Bloch bands as transformed Bloch orbitals.

When considering dual energy and spatial localization, we show two possible
ways to evaluate the energy of the Wannier functions needed for minimizing the
energy variance. First, we show how the energy variance can be evaluated using
the true Bloch orbital energies. Second, we demonstrate how using the
disentangled Bloch orbital energies creates a simpler and more readily
implementable representation. Finally, we detail the regime where these two
approaches agree.

###  II.1 Bloch orbital basis

Following the procedure for disentangling a set of Bloch orbitals
$\\{\lvert\psi_{n}^{\mathbf{k}}\rangle\\}$ outlined by [8] results in a set of
disentangled Bloch orbitals $\\{\lvert\varphi_{n}^{\mathbf{k}}\rangle\\}$. A
disentangled band can be represented in the Bloch orbital basis as

$\lvert\varphi_{n}^{\mathbf{k}}\rangle=\sum_{b}^{N_{b}}V_{nb}^{\mathbf{k}}\lvert\psi_{b}^{\mathbf{k}}\rangle,$
(12)

where $V^{\mathbf{k}}$ is an $N_{w}$$\times$$N_{b}$ matrix obeying the
relation $V^{\mathbf{k}}(V^{\mathbf{k}})^{\dagger}=\mathbb{1}$ ($\mathbb{1}$
is the $N_{w}$$\times$$N_{w}$ identity operator). We may then write a
transformed Bloch orbital for generalized Wannier function construction in the
basis of disentangled Bloch orbitals as

$\lvert\phi_{n}^{\mathbf{k}}\rangle=\sum_{a}^{N_{w}}U_{na}^{\mathbf{k}}\lvert\varphi_{b}^{\mathbf{k}}\rangle.$
(13)

A generalized Wannier function constructed from a set of disentangled Bloch
orbitals then takes the form

$\begin{split}\lvert
w_{m}^{\mathbf{R}}\rangle&=\frac{1}{N_{k}}\sum_{\mathbf{k}}^{N_{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\sum_{a}^{N_{w}}U_{ma}^{\mathbf{k}}\sum_{b}^{N_{b}}V_{ab}^{\mathbf{k}}\lvert\psi_{a}^{\mathbf{k}}\rangle\\\
&=\frac{1}{N_{k}}\sum_{\mathbf{k}}^{N_{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\sum_{b}^{N_{b}}L_{mb}^{\mathbf{k}}\lvert\psi_{b}^{\mathbf{k}}\rangle,\end{split}$
(14)

where $L_{mb}^{\mathbf{k}}=[U^{\mathbf{k}}V^{\mathbf{k}}]_{mb}$.

Using this formulation, we may break the energy cost function $\Xi$ into a
transformation-invariant term similar to $\Omega_{\text{I}}$, but where the
position operator $\mathbf{r}$ is replaced with $h$:

$\begin{split}\Xi_{\text{I}}&=\sum_{n}^{N_{w}}\Big{(}\langle
w_{n}^{\mathbf{0}}\rvert h^{2}\lvert
w_{n}^{\mathbf{0}}\rangle-\sum_{m}^{N_{w}}\sum_{\mathbf{R}}^{N_{k}}\left|\langle
w_{m}^{\mathbf{R}}\rvert h\lvert
w_{n}^{\mathbf{0}}\rangle\right|^{2}\Big{)}\\\
&=\frac{1}{N_{k}}\sum_{\mathbf{k}}^{N_{k}}\bigg{\\{}\operatorname{Tr}\Big{[}P^{\mathbf{k}}\big{(}H^{\mathbf{k}}\big{)}^{2}\Big{]}-\operatorname{Tr}\Big{[}\big{(}P^{\mathbf{k}}H^{\mathbf{k}}\big{)}^{2}\Big{]}\bigg{\\}},\end{split}$
(15)

where $P^{\mathbf{k}}=(V^{\mathbf{k}})^{\dagger}V^{\mathbf{k}}$ is the
projector from the space spanned by the $N_{b}$ Bloch orbitals to the
$N_{w}$-dimensional space of disentangled Bloch orbitals. If the commutator
$[P^{\mathbf{k}},H^{\mathbf{k}}]=0$, then by the cyclic property of the trace
$\Xi_{\text{I}}=0$. This happens when the disentangled Bloch orbitals are a
subset of the true Bloch orbitals, or for any mixing of degenerate Bloch
orbitals at a fixed $\mathbf{k}$-point.

###  II.2 Disentangled orbital basis

The energy spread cost function may also be evaluated using the energy of the
disentangled Bloch orbitals. In the procedure outlined in [8], the energy of
the disentangled subspace is diagonalized in the space spanned by the
disentangled Bloch orbitals. This means that

$\big{[}\widetilde{H}^{\mathbf{k}}\big{]}_{nm}=\langle\varphi_{n}^{\mathbf{k}}\rvert\tilde{h}\lvert\varphi_{m}^{\mathbf{q}}\rangle=N_{k}\delta_{\mathbf{k}\mathbf{q}}\delta_{nm}\tilde{\epsilon}_{n}^{\mathbf{k}},$
(16)

where
$\widetilde{H}^{\mathbf{k}}=V^{\mathbf{k}}H^{\mathbf{k}}(V^{\mathbf{k}})^{\dagger}$.
If we replace the full Hamiltonian $h$ by the disentangled Hamiltonian
$\tilde{h}$, then the projector $P^{\mathbf{k}}$ in Eq. (15) becomes the
identity matrix and $\Xi_{\text{I}}$ vanishes. Therefore, even if a
disentanglement leads to a nonzero $\Xi_{\text{I}}$, using $\tilde{h}$ instead
of the full $h$ in the definition of $\Xi$ yields $\Xi_{\text{I}}=0$. Using
the disentangled Hamiltonian also means the $\tilde{h}^{2}$ terms are
unitarily invariant:

$\begin{split}\sum_{n}^{N_{w}}\langle
w_{n}^{\mathbf{0}}\rvert\tilde{h}^{2}\lvert
w_{n}^{\mathbf{0}}\rangle=\frac{1}{N_{k}}\sum_{\mathbf{k}}\operatorname{Tr}\Big{[}(\widetilde{H}^{\mathbf{k}})^{2}\Big{]}.\end{split}$
(17)

This allows us to break the energy cost term into the squared average energy
term, $\Xi_{\text{SA}}=\sum_{n}\langle w_{n}\rvert\tilde{h}^{2}\lvert
w_{n}\rangle$, and an average energy squared term,
$\Xi_{\text{AS}}=\sum_{n}|\langle w_{n}\rvert\tilde{h}\lvert
w_{n}\rangle|^{2}$, such that the total energy spread cost is
$\Xi=\Xi_{\text{SA}}-\Xi_{\text{AS}}$. Similar to the spatial variance in the
$\Gamma$-sampled case [25], this means that minimizing the energy spread is
equivalent to maximizing the average energy squared $\Xi_{\text{AS}}$.

The difference between using the projected subspace Hamiltonian $\tilde{h}$
and the full Hamiltonian $h$ depends on the set of matrices
$\\{V^{\mathbf{k}}\\}$. If the disentangled multiset of eigenvalues is a
subset of the Bloch orbital eigenvalues at every $\mathbf{k}$-point,
$\\{\tilde{\epsilon}_{n}^{\mathbf{k}}\\}\subseteq\\{\epsilon_{n}^{\mathbf{k}}\\}$,
then $[\widetilde{H}^{\mathbf{k}},H^{\mathbf{k}}]=0$ and using $\tilde{h}$ is
equivalent to using $h$ in the energy cost function. A multiset is a set that
accounts for multiplicity, which is necessary to account for degenerate
eigenvalues. We find that we can use the existing disentanglement procedure to
create a disentangled band structure that is almost equivalent to the full
band structure in the energy range of interest. For example, we can find
Wannier functions around the Fermi energy by including a sufficiently large
number of unoccupied bands above the Fermi energy when applying the
disentanglement procedure. This produces a disentangled band structure that
only appreciably differs from the original band structure far away from the
Fermi energy. This is because the disentanglement procedure implemented in
wannier90 diagonalizes $\tilde{h}$ at the end of each iteration in the
minimization [8]. If the projected subspace fully spans the area of interest
around the Fermi energy, then the eigenvalues of $\tilde{h}$ will correspond
to those of $h$ near that energy value. We do note that it is possible to
include $\Xi_{\text{I}}$ in the disentanglement procedure, so that the
resulting disentangled eigenspectrum would more closely resemble the full
eigenspectrum for all energy values.

###  II.3 Energy cost gradient

Following the breakdown of the cost function into $\Xi_{\text{SA}}$ and
$\Xi_{\text{AS}}$ as outlined in the previous section, we now derive the
gradient of $\Xi_{\text{AS}}$ with respect to the unitary rotation
$U^{\mathbf{k}}$. We approximate $U^{\mathbf{k}}$ to first order using a small
anti-Hermitian matrix, $U^{\mathbf{k}}\approx\mathbb{1}+W^{\mathbf{k}}$, with
$(W^{\mathbf{k}})^{\dagger}=-W^{\mathbf{k}}$. Then we find the derivative with
respect to $U^{\mathbf{k}}$ using the matrix calculus identity

$\frac{d\,\operatorname{Re}\big{\\{}\operatorname{Tr}\big{[}M\,W\big{]}\big{\\}}}{dW}=\frac{1}{2}\big{(}M-M^{\dagger}\big{)}=\mathcal{A}[M].$
(18)

To see why this is useful, observe that we may write $\Xi_{\text{AS}}$ as

$\begin{split}\Xi_{\text{AS}}&=\sum_{n}\left\lvert\langle
w_{n}|\tilde{h}|w_{n}\rangle\right\rvert^{2}\\\
&=\sum_{n}\Big{\lvert}\frac{1}{N_{k}}\sum_{m,\mathbf{k}}|U_{mn}^{\mathbf{k}}|^{2}\tilde{\epsilon}_{m}^{\mathbf{k}}\Big{\rvert}^{2}\\\
&=\sum_{n}\Big{|}\frac{1}{N_{k}}\sum_{\mathbf{k}}B_{nn}^{\mathbf{k}}\Big{|}^{2},\end{split}$
(19)

where
$[B^{\mathbf{k}}]_{mn}=[U^{\mathbf{k}}\widetilde{H}^{\mathbf{k}}(U^{\mathbf{k}})^{\dagger}]_{mn}$
is the Hamiltonian in the transformed Bloch orbital basis. After transforming
the Bloch functions at each $\mathbf{k}$-point by
$U^{\mathbf{k}}\approx\mathbb{1}+W^{\mathbf{k}}$, the change in
$\Xi_{\text{AS}}$ to first order in $W^{\mathbf{k}}$ is

$\begin{split}d\Xi_{\text{AS}}&=\frac{4}{N_{k}^{2}}\sum_{n}\Big{(}\sum_{\mathbf{k}}\operatorname{Re}\big{\\{}[B^{\mathbf{k}}W^{\mathbf{k}}]_{nn}\big{\\}}\Big{)}\Big{(}\sum_{\mathbf{q}}B_{nn}^{\mathbf{q}}\Big{)}\\\
&=\frac{4}{N_{k}^{2}}\sum_{\mathbf{k}\mathbf{q}}\operatorname{Re}\big{\\{}\operatorname{Tr}[C^{(\mathbf{k},\mathbf{q})}W^{\mathbf{k}}]\big{\\}}\end{split}$
(20)

where
$C^{(\mathbf{k},\mathbf{q})}_{mn}=B^{\mathbf{k}}_{mn}B^{\mathbf{q}}_{mm}$.
This allows us to write the gradient in terms of $\mathcal{A}$, defined in Eq.
(18), as

$\begin{split}\frac{d\Xi_{\text{AS}}}{dW^{\mathbf{k}}}&=\frac{4}{N_{k}}\sum_{\mathbf{q}}\mathcal{A}[C^{(\mathbf{k},\mathbf{q})}]\\\
&=\frac{2}{N_{k}}\sum_{\mathbf{q}}\Big{(}C^{(\mathbf{k},\mathbf{q})}-\left(C^{(\mathbf{k},\mathbf{q})}\right)^{\dagger}\Big{)}.\end{split}$
(21)

In this form it appears the gradient at each $\mathbf{k}$-point depends on
every other $\mathbf{k}$-point, which would indicate $N_{k}^{2}$ scaling with
the number of $\mathbf{k}$-points. However, we can sum the matrix
$C^{(\mathbf{k},\mathbf{q})}$ over the Brillouin zone, removing the dependence
on $\mathbf{q}$ prior to evaluating the gradient. Thus

$\begin{split}\bigg{[}\frac{d\Xi_{\text{AS}}}{dW^{\mathbf{k}}}\bigg{]}_{mn}&=2B^{\mathbf{k}}_{mn}\left(E_{m}-E_{n}\right),\end{split}$
(22)

where $E_{m}=\langle w_{m}^{\mathbf{0}}\rvert\tilde{h}\lvert
w_{m}^{\mathbf{0}}\rangle$ is the average energy of a Wannier function. In
essence, these terms penalize the mixing of Bloch orbitals proportionally to
the difference in average energy of the Wannier functions they are used to
construct. This also means the computational cost of constructing an element
of $d\Xi_{\text{AS}}/dW^{\mathbf{k}}$ is independent of $N_{k}$.

###  II.4 Mixing parameter

We now remark on the mixing parameter $\gamma$ introduced in Eq. (6). Setting
$\gamma=0$ recovers the MLWF cost function $\Omega$. Setting $\gamma=1$
returns the Bloch orbitals in the case of $\Gamma$ sampling. However, when
$N_{k}>1$, setting $\gamma=1$ yields transformed Bloch orbitals that are
simply the Bloch orbitals ordered by their energy at each $\mathbf{k}$-point.
When there are band crossings in the Brillouin zone, the energy-ordered Bloch
orbitals will not be smooth in $\mathbf{k}$, giving poor spatial localization.
Choosing $\gamma$ strictly between 0 and 1 provides Wannier functions
localized in both space and energy.

In a related work for molecules, this cost function was used in an orbital-
dependent density functional method called the localized orbital scaling
correction (LOSC) [23]. In that work, $\gamma$ was chosen to minimize the
error in several experimentally realizable quantities in a test suite of
molecules. In Å and eV, the units of space and energy used in wannier90, this
value is $\gamma=0.47714$. For details on this value in relation to the value
published in [23], see the Supplemental Material.

When using localized orbitals in the context of LOSC [24, 23], we refer to
them as orbitalets, to highlight their compromise of spatial and energy
localization, following the concept of wavelets having compromise of spatial
and momentum localization [26]. A key difference between MLWFs and DLWFs is
that DLWFs naturally allow the inclusion of unoccupied orbitals in their
construction. Energy cutoffs must be enforced manually for MLWFs, because
including high-energy bands will result in unphysically localized Wannier
functions; indeed, including the full space of occupied and unoccupied Bloch
orbitals is expected to yield arbitrarily localized MLWFs, approaching Dirac
delta distributions as the number of unoccupied bands is increased. By
contrast, for suitable values of $\gamma$, the energy localization of DLWFs
induce convergence to stable frontier orbitals when enough higher-energy
orbitals are included. When $\gamma=0.47714$, we find that Bloch orbitals can
mix substantially in the construction of DLWFs at a given $\mathbf{k}$-point
as long as their energy difference is less than about $2\text{\,}\mathrm{eV}$.

###  II.5 Computational details

Following the MLWF formulation of Marzari and Vanderbilt [9], we compute the
gradient of the cost function $F$ at each $\mathbf{k}$-point. Given an initial
guess or disentanglement, either a conjugate gradient or steepest descent
algorithm is used to minimize $F$. We implement the energy-space localization
in a fork of the open-source wannier90 code [27, 28, 29]. In order to include
$\Xi$, some modifications to the descent algorithm were required. The steepest
descent portion of the algorithm was left unchanged, but we found that the
step size that produced the best localization was system-dependent. To account
for this, we sweep a range of step sizes in order to obtain the best minimum
for each localization. For steps when the conjugate-gradient descent and
parabolic line search were used, we found that the Polak–Ribiere coefficient
[30] provides better convergence than the default Fletcher–Reeves coefficient
[31, 32].

Because we allow the inclusion of unoccupied orbitals, we also use the
disentanglement procedure on the highest-energy conduction bands considered
for localization [8]. When including virtual bands, the choice of how many
bands and DLWFs to include is increased until the orbitals of interest,
typically the ones around the Fermi energy, are converged. Convergence of
high-lying virtual bands was found to be difficult; to sidestep this issue, we
implement a program option allowing the exclusion of specified bands from the
cost function convergence criterion. We consider the DLWFs converged if their
cost from Eq. (6) does not appreciably change with the addition of more
unoccoupied bands. The modified version of the code may be found at [33] and
more details on the functionality added to the wannier90 code are in the
Supplemental Material.

##  III Results

We obtain Bloch functions using the PBE density functional [34] in the plane-
wave basis, using the optimized norm-conserving Vanderbilt pseudopotentials
[35] generated from PseudoDojo [36]; calculations are performed using the
open-source Quantum ESPRESSO code suite [37, 38]. The Brillouin zone is always
sampled with Monkhorst-Pack meshes [15] centered at $\Gamma$, the origin of
reciprocal space. For silicon and ethylene, we use a wave function kinetic
energy cutoff $E_{\text{cut}}=60$ Ry; for copper, we set $E_{\text{cut}}=100$
Ry. In all cases $E_{\text{cut}}$ for the density is four times that used for
the wave functions. As mentioned in Sec. II, we use a modified version of the
wannier90 code [27, 28, 29] to obtain the DLWFs. We select the LOSC mixing
parameter $\gamma=0.47714$, unless stated otherwise. Since the DLWFs are not
always real, isoplots are of the DLWF densities, $|w_{n}(\mathbf{r})|^{2}$,
rather than of the DLWFs themselves. All lengths are reported in
$\mathrm{\SIUnitSymbolAngstrom}$ and all energies in $\mathrm{eV}$.

###  III.1 Silicon

First we consider the well-studied semiconductor silicon in the diamond
lattice. For the self-consistent calculation we use an 8$\times$8$\times$8
$\mathbf{k}$-mesh; for computational efficiency, we use a 4$\times$4$\times$4
$\mathbf{k}$-mesh for obtaining virtual Bloch orbitals and localization. We
use the experimental lattice parameter
$a_{\text{lat}}=$5.431\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ [39]. The gap
of Si is small, with the self-consistent calculation yielding
$0.71\text{\,}\mathrm{eV}$. It is an indirect gap, however, and the smallest
direct gap in the Brillouin zone is $2.56\text{\,}\mathrm{eV}$ at the $\Gamma$
point.

####  III.1.1 Occupied states

Constructing DLWFs from only the occupied bands of silicon yields Wannier
functions which are not degenerate in energy and have different shapes; this
contrasts with the MLWF procedure, which yields four degenerate Wannier
functions related by space-group symmetry operations. The lowest-energy DLWF
is tetrahedral, while the two highest-energy orbitalets are degenerate, each
having three lobes centered along a Si–Si bond. The MLWFs, on the other hand,
all have tight-binding character, with a single lobe centered along a bond.

####  III.1.2 Frontier States

We expect the valence bands of semiconductors to yield degenerate MLWFs that
approximate bonding molecular orbitals; localizing the same number of low-
lying conduction bands with the MLWF procedure often yields degenerate
functions of antibonding character. Since the conduction states are entangled
with higher-energy bands, Souza  _et al._ [8] disentangled 12 bands down to 8
when investigating silicon. Here, we use the same disentanglement with a
frozen disentanglement window at the Fermi energy, $6.23\text{\,}\mathrm{eV}$.
Using Eq. (6) for the cost function on this subspace, we obtain occupied
orbitals very similar to those found in Sec. III.1.1 above. The conduction
orbitals, however, form a fourfold degenerate anti-bonding set, qualitatively
equivalent to MLWFs when considering only the four lowest-lying conduction
bands.

We also use this system to explore the effects of varying values of $\gamma$
in the cost function. Allowing both valence and conduction bands to mix in the
MLWF construction yields eight Wannier functions degenerate in energy. As we
increase $\gamma$, we observe WFs more closely associated with specific energy
ranges and stepwise lifting of degeneracy (see Fig. 1), but spatial
localization is suppressed at the same time (see Fig. 2).

Figure 1: Degeneracy pattern of the average energy $\langle\tilde{h}\rangle$
of the DLWFs in silicon against $\gamma$ from Eq. (6). The color and line type
indicate the degeneracy of the DLWFs; we call $\lvert w_{m}\rangle$ and
$\lvert w_{n}\rangle$ degenerate if
$|\langle\tilde{h}\rangle_{n}-\langle\tilde{h}\rangle_{m}|<$0.15\text{\,}\mathrm{eV}$$.
At $\gamma=0$ (MLWFs), we obtain degenerate $sp^{3}$ hybrid orbitals, with
each MLWF having occupancy $\lambda_{nn}=0.5$ and $\langle h\rangle=6.3$ eV.
In the range $0.02<\gamma<0.05$, the DLWF degeneracy pattern rapidly changes
before splitting into two degenerate sets between $0.05\leq\gamma\leq 0.09$:
the occupied and unoccupied orbitals. The degeneracy of unoccupied DLWFs is
lifted at $\gamma=0.53$, while the first lifting of degeneracy for occupied
DLWFs occurs at $\gamma=0.1$. Figure 2: Spatial cost $\Omega$, energy cost
$\Xi$, and total cost $F$ against $\gamma$ in Eq. (6). Since the total cost is
a convex sum of $\Omega$ and $\Xi$, it always remains between their two
values. The sharp changes in $\Omega$ and $\Xi$ correspond to degeneracy
splitting in Fig. 1. Note the discontinuity in the vertical axis.

####  III.1.3 Converged frontier states

Next, we show convergence results when higher-energy virtual bands are
included. For this result, we find that the spatial variance of the highest-
energy occupied DLWF is well converged when 34 bands are disentangled to 30.
Notably, we find that the shapes of the resulting DLWFs are qualitatively the
same as those found when including only occupied orbitals (frontier orbitals),
as in Sec. III.1.1 (Sec. III.1.2). We note, however, that the spatial variance
$\Delta r^{2}$ of the lowest unoccupied DLWFs is larger than when only four
unoccupied bands are included. This is due to the fact that the
disentanglement in the latter case directly modifies the conduction bands,
smoothing them in $\mathbf{k}$-space and unphysically increasing the
localization of the resulting Wannier functions. Including more virtual
orbitals in the localization procedure means that disentanglement smooths only
the high-energy conduction bands, removing this effect. For isosurface plots
of the eight lowest-energy DLWF densities, see Fig. 3. For more data on
individual DLWFs, see the Supplemental Material.

Figure 3: DLWF density isoplots for silicon (isovalue $0.6$) using
$\gamma=0.47714$. They are labeled in ascending order of DLWF average energy
starting at the lowest energy. Note that even though this is the result of
using 30 disentangled Bloch orbitals per $\mathbf{k}$-point for construction,
the degeneracy pattern is the same as the frontier-only case show in Fig. 1
for the value of $\gamma=0.47714$.

###  III.2 Copper

We next investigate copper in an face-centered cubic lattice using an
experimental lattice parameter of
$3.614\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ [40]. For the self-consistent
calculation we use a 16$\times$16$\times$16 $bk$-mesh, while for the virtual
states and localization we use a 10$\times$10$\times$10 $\mathbf{k}$-mesh.
Since the unoccupied bands of copper are very steep, and energy localization
restricts mixing bands far apart in energy, we do not require many bands above
the Fermi level. Noting that there are 9.5 electrons per spin channel, we use
23 bands and disentangle to 17 bands with a frozen disentanglement window of
$40\text{\,}\mathrm{eV}$, well above the Fermi energy
($17.04\text{\,}\mathrm{eV}$). This results in a set of disentangled bands
that only differs from the parent DFA at the top of the energy window. Using
the $\gamma=0.47714$, we find that the fully occupied orbitals self-organize
into pure $s$, $p$, and $d$-type DLWFs (see Fig. 4). The partially occupied
DLWF associated with the frontier band has $d$-type character, but is
delocalized across multiple atoms. Information on individual DLWFs is provided
in the Supplemental Material.

Figure 4: DLWF density isoplots of copper at isovalue $0.6$ using
$\gamma=0.47714$. These are the six highest-energy occupied DLWFs labeled in
ascending order of average energy, starting with the fifth-lowest energy. The
four DLWFs of lowest energy (not shown) correspond to the occupied $3s$ and
$3p$ states (not shown). DLWF 10 is only partially occupied, as it is mainly
comprised of the frontier Bloch band. More atoms are shown for the isoplot of
DLWF 10 due to its delocalized nature.

###  III.3 Ethylene

To demonstrate the utility of dual localization in molecules, we simulate
ethylene with the geometry from [9]. We use a
$10\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ unit cell with the molecule
centered at the origin and sample only the $\Gamma$. Thus, disentanglement is
unnecessary, and we construct 36 DLWFs from 36 bands, 6 of which are occupied.
Setting $\gamma=0.47714$ recovers DLWFs closely related to molecular orbitals.
The lowest-energy DLWF corresponds to a $\sigma$ bond between the carbon
atoms, and the next four resemble combinations of C–C and C–H bonding
orbitals. The highest-energy occupied DLWF corresponds to a $\pi$ bonding
orbital, and the lowest unoccupied DLWF to a $\pi^{*}$ antibonding orbital.
Their isoplots are shown in Fig. 5, and the individual DLWF data are in the
Supplemental Material.

Figure 5: DLWF density isoplots for ethylene using $\gamma=0.47714$. DLWF 1
uses an isovalue of 160, DLWF 2–DLWF 7 use an isovalue of 60, and DLWF 8 uses
an isovalue of 30. The isosurface between DLWF 6 and DLWF 8 belongs to DLWF 8.

##  IV Discussion

We have shown that including energy localization in the construction of
Wannier functions results in DLWFs localized both in space and in energy. This
allows for bands widely separated in energy to be included in the Wannier
function construction without allowing the unphysical mixing of bands far
apart in energy; furthermore, the localization procedure naturally produces
orbitals associated with certain energies, organizing itself along the
Hamiltonian’s spectrum. DLWFs with occupation close to $1$ also retain the
behavior seen in MLWFs, appearing as tight-binding or atomic orbitals and
retaining the symmetry of the underlying system. Using the LOSC cost function,
$\gamma=0.47714$, the localization identifies atomic states such as the $d$
orbitals of copper (figure 4). The LOSC cost function also recovers molecular
orbitals such as the ethylene $\sigma$ and $\pi$ bonding orbitals in Fig. 5.
This shows that DLWFs have the potential to produce chemically relevant
frontier orbitals automatically, without the need for manually enforced energy
cutoffs or for separating the occupied and unoccupied spaces.

The fact that DLWFs offer a frontier orbital description localized in space
means that they can be employed for orbital-dependent DFT methods [11, 12, 13,
14, 23]. Many of these methods attempt to correct the well-known failure of
density functional approximations to give energy piecewise linear in
fractional occupation [41]. In order to correct this nonlinear behavior, many
of these corrections invoke localized frontier orbitals. A popular choice for
these orbitals is to use MLWFs constructed from the highest-energy composite
set of occupied Bloch bands. The fractional occupations in DLWFs are a key
features distingushing them from other Wannier functions, and play an
important role in the LOSC correction to density functional approximations
[24, 23]. As DLWFs are explicitly associated with an energy level, they offer
a localized charge description of frontier orbitals for both the occupied and
unoccupied spaces. This means they can also generate relevant Wannier
functions for both occupied and unoccupied frontier orbitals without having to
localize each separately.

When the Wannier function construction includes entangled bands, Damle and Lin
[42] have questioned whether sequential disentanglement and localization can
find the globally optimal solution as separate procedures. The unified
localization method proposed therein offers a possible solution to this
problem. For a composite set of bands, this method reduces to the selected
columns of the density matrix method [43, 44]; in that case, as with the
original MLWF construction, energy cutoffs can be enforced explicitly to limit
the energy spread of the Wannier functions. However, because the unified
method of Damle and coworkers relies on the exponential decay of off-diagonal
elements of the one-particle density matrix in space, it is unclear whether it
can be adapted to include energy localization in the construction of the
Wannier functions while allowing the inclusion of Bloch orbitals widely
separated in energy, as shown in this work.

Finally, we comment on the (lack of) reality of the DLWFs. When $\gamma=0$, so
that $F=\Omega$ in Eq. (6), Marzari and Vanderbilt conjectured in [9] that the
resulting maximally localized Wannier functions would be purely real (up to a
global phase). This was proven in [21] for Wannier functions constructed by
minimizing any functional symmetric under time reversal, provided that the
optimum is unique. Our cost function $F$ obeys time-reversal symmetry, and
numerical tests show that the DLWFs with an occupation of $1$ are real.
However, we observe substantial imaginary character in DLWFs that have an
occupation numerically distinguishable from $1$. Even though the arguments in
[21] do not differentiate between the occupied and unoccupied spaces, we
conjecture that the reality arguments do not apply to Wannier functions whose
occupation is far from 1. It is possible that the non-uniqueness of these
minima could account for their lack of reality.

##  V Acknowledgements

All of the authors acknowledge support from the National Science Foundation
(Grant No. CHE-1900338); A.M. from the Molecular Sciences Software Institute
Phase-II Software Fellowship; and J.Z.W. from the National Institutes of
Health (Grant No. 5R01GM061870).

## References

  * Wannier [1937] G. Wannier, Physical Review 52, 191 (1937).
  * Silvestrelli and Parrinello [1999a] P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. 82, 3308 (1999a), publisher: American Physical Society.
  * Silvestrelli and Parrinello [1999b] P. L. Silvestrelli and M. Parrinello, The Journal of Chemical Physics 111, 3572 (1999b), publisher: American Institute of PhysicsAIP.
  * Lee  _et al._ [2005] Y.-S. Lee, M. B. Nardelli, and N. Marzari, Phys. Rev. Lett. 95, 076804 (2005), publisher: American Physical Society.
  * Li and Marzari [2011] E. Y. Li and N. Marzari, ACS Nano 5, 9726 (2011), publisher: American Chemical Society.
  * Shelley  _et al._ [2011] M. Shelley, N. Poilvert, A. A. Mostofi, and N. Marzari, Computer Physics Communications 182, 2174 (2011).
  * Goedecker [1999] S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999).
  * Souza  _et al._ [2001] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001), publisher: American Physical Society.
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997), publisher: American Physical Society.
  * Resta [2000] R. Resta, Berry’s phase and geometric quantum distance (2000), publisher: EPFL. Lausanne.
  * Stengel and Spaldin [2008] M. Stengel and N. A. Spaldin, Phys. Rev. B 77, 155106 (2008), publisher: American Physical Society.
  * Borghi  _et al._ [2014] G. Borghi, A. Ferretti, N. L. Nguyen, I. Dabo, and N. Marzari, Phys. Rev. B 90, 075135 (2014), publisher: American Physical Society.
  * Ma and Wang [2016] J. Ma and L.-W. Wang, Sci Rep 6, 24924 (2016).
  * Wing  _et al._ [2020] D. Wing, G. Ohad, J. B. Haber, M. R. Filip, S. E. Gant, J. B. Neaton, and L. Kronik, arXiv:2012.03278 [cond-mat] (2020), arXiv: 2012.03278.
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Ashcroft and Mermin [2003] N. Ashcroft and D. Mermin,  _Solid State Physics_ (Thomson Press, 2003).
  * Cloizeaux [1964] J. D. Cloizeaux, Phys. Rev. 135, A685 (1964).
  * Cloizeaux [1963] J. D. Cloizeaux, Phys. Rev. 129, 554 (1963), publisher: American Physical Society.
  * Foster and Boys [1960] J. M. Foster and S. F. Boys, Rev. Mod. Phys. 32, 300 (1960), publisher: American Physical Society.
  * Kohn [1959] W. Kohn, Phys. Rev. 115, 809 (1959).
  * Brouder  _et al._ [2007] C. Brouder, G. Panati, M. Calandra, C. Mourougane, and N. Marzari, Phys. Rev. Lett. 98, 046402 (2007), publisher: American Physical Society.
  * Panati and Pisante [2013] G. Panati and A. Pisante, Commun. Math. Phys. 322, 835 (2013).
  * Su  _et al._ [2020] N. Q. Su, A. Mahler, and W. Yang, J. Phys. Chem. Lett. 11, 1528 (2020).
  * Li  _et al._ [2018] C. Li, X. Zheng, N. Q. Su, and W. Yang, National Science Review 5, 203 (2018).
  * Gygi  _et al._ [2003] F. Gygi, J.-L. Fattebert, and E. Schwegler, Computer Physics Communications 155, 1 (2003).
  * Daubechies [1992] I. Daubechies,  _Ten Lectures on Wavelets_, CBMS-NSF Regional Conference Series in Applied Mathematics (Society for Industrial and Applied Mathematics, 1992).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 178, 685 (2008).
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, J. Phys.: Condens. Matter 32, 165902 (2020).
  * Polak and Ribiere [1969] E. Polak and G. Ribiere, R.I.R.O. 3, 35 (1969).
  * Fletcher and Reeves [1964] R. Fletcher and C. M. Reeves, The Computer Journal 7, 149 (1964).
  * Press  _et al._ [1992] W. Press, S. Taukolsky, W. Vetterling, and B. Flannery,  _Numerical Recipes in Fotran, The Art of Scientific Computing_ , 2nd ed. (Cambridge University Press, 1992).
  * [33] github.com/mtesseracted/wannier90/tree/costSE.
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Hamann [2013] D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * van Setten  _et al._ [2018] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, Computer Physics Communications 226, 39 (2018).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009), publisher: IOP Publishing.
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. d. Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys.: Condens. Matter 29, 465901 (2017), publisher: IOP Publishing.
  * Hom  _et al._ [1975] T. Hom, W. Kiszenik, and B. Post, J Appl Crystallogr 8, 457 (1975).
  * Rutt  _et al._ [2006] O. J. Rutt, G. R. Williams, and S. J. Clarke, Chem. Commun. , 2869 (2006).
  * Perdew  _et al._ [1982] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, PHYSICAL REVIEW LETTERS 49, 4 (1982).
  * Damle and Lin [2018] A. Damle and L. Lin, Multiscale Model. Simul. 16, 1392 (2018).
  * Damle  _et al._ [2015] A. Damle, L. Lin, and L. Ying, J. Chem. Theory Comput. 11, 1463 (2015).
  * Damle  _et al._ [2017] A. Damle, L. Lin, and L. Ying, Journal of Computational Physics 334, 1 (2017).