# A Review of the GIC Blocker Placement Problem

Arthur K. Barnes 1*  —, Adam Mate 1*  —, and Russell Bent 1*  — Manuscript
submitted: Oct. 09, 2023. Current version: Mar. 20, 2024. 1 The authors are
with the Advanced Network Science Initiative at Los Alamos National
Laboratory, Los Alamos, NM 87545 USA.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Approved for public release; distribution is unlimited.

###### Abstract

Space weather poses a tremendous threat to power systems: geomagnetic
disturbances could result in widespread disruptions and long-duration
blackouts, including severe damage to system components. To mitigate their
impacts, a handful of strategies exist, with the most promising being the
deployment of transformer neutral blocking devices. The high cost of these
devices, however, precludes their installation at all substations; this
motivates the development of effective solutions for the cost-effective
placement of such devices. While the current state-of-the-art in blocker
placement methods is insufficient to be applied to real-sized power grids,
ongoing research continues to increase the size of networks for which the
placement problem remains tractable. Along these lines, the contributions of
this paper are two fold: first, a comprehensive overview of the current state-
of-the-art in blocker placement methods is provided; and second, a complete
optimization formulation – implemented and benchmarked in an open-source
software – for the blocker placement problem is presented.

###### Index Terms:

geomagnetic disturbance, power systems resilience, power systems optimization,
mitigation strategy, transformer neutral blocking device.

##  I Introduction

The bulk electric system (BES) faces risks on many different time-scales –
both short- and long-term events, varying based on the amount of warning time
present before the threat and the time expected to restore service after the
threat – which can lead to major system disruptions [1]. From a grid
resilience perspective, the most consequential and hardest to predict risks
are called high-impact low-frequency events (HILFs) [2]. HILFs are a challenge
for protecting the BES because they may cause widespread, long duration
blackouts – potentially involving transient instability issues, cascading
outages – or severe damage to grid components, such as large power
transformers, which are challenging and time-consuming to repair [3, 4].

The North American Electric Reliability Corporation identifies geomagnetic
disturbances (GMDs) as HILFs [2]: these solar-driven incidents disrupt the
Earth’s magnetic field and drive geomagnetically induced currents (GICs) in
the conductive infrastructure, and may be predicted a few hours or few days
before the threat. Since they are rare, there is a lack of empirical data on
how the BES will respond to such threats; GMDs change the physics of grid
operation, so power-flow formulations must be modified [5]. Therefore,
accurately simulating their impact and developing cost-effective mitigation
strategies that can operate on their time-scales is a complex research
challenge.

###  I-A Approaches to Mitigate GMD Impacts

The research community has developed approaches over the years to help plan
for and mitigate the potential impacts of GMDs on the BES. Considered
mitigation strategies include load-shedding [6], transmission line switching
[7, 8, 9], generation switching (re-dispatch) [6], shunt/series-capacitor
switching [10], and addition of transformer neutral (dc-current) blocking
devices [11].

Line switching reduces the total transmission capacity of the power grid, and
likely requires generation re-dispatch, customer load-shedding, or both to
achieve a new generation-demand balance. Generation switching, removing
generation from the grid, breaks GIC flow via disconnecting generator step-up
transformers. Transformer neutral blocker placement – the installation of GIC
dc-current blocking devices – consists of injecting a shunt capacitor in
series with the grounding point of transformer neutrals; this does not
interfere with the power-carrying capacity of the grid, but may interfere with
ground-fault relaying, which in turn requires additional switching devices
leading to increased capital costs.

This paper focuses on the placement of transformer neutral blocking devices to
mitigate GMD impacts. Given the high cost of these devices, sparse placement
is necessary; it is an unrealistic expectation to install such devices at all
transformer neutrals in all substations. In addition, the literature has
suggested that the number of transformers that experience high GICs during GMD
events is also sparse [12]. Finally, the placement itself is a challenging
problem because installing a blocker at one transformer can increase GICs at
other transformers throughout the power grid; therefore, simple approaches,
such as greedy methods, are likely unsuitable in this case. These three
factors motivate the need for optimization-based placement methods to guide
the choices of where to place blocking devices.

TABLE I: Summary of GIC Blocker Placement Methods Problem Formulation | Objective | Solver | Busses | Runtime | CPU | Authors
---|---|---|---|---|---|---
dc GIC flow | min. $I_{eff}$ | GA |  100 | – | – | [13]
dc GIC flow | min. $I_{eff}$ | GA |  110 | – | – | [11]
time-extended dc GIC flow | avg./peak $q_{loss}$ min. over time | SA | – | – | – | [14]
SDP, relaxed placement | min. total $q_{loss}$ or max. $q_{loss}$ | – | 20 | – | – | [15]
SOCP, relaxed placement, lin. $q_{loss}$ | – | – | – | – | – | [16]
lin. sensitivity factors, line switching | – | greedy | – | 0.29 s | Core i7 | [7]
scen.-based lin. $q_{loss}$, ac power flow | min. blocker count | GA |  1,000 | – | – | [17]
– | min. blocker count | surrogate opt. | 118 | – | – | [18]
MINLP, ac power flow | – | GA | – | 6 h | Core i7-2600 | [19]
MINLP, ac power flow | – | branch-and-cut | 1,875 | 2:47 h | Core i5 | [20]
MINLP, no ac power flow | – | SA | 1,875 | 48 h | 3.10GHz Core i5 | [14]
MINLP, harm. ac power flow, xfmr heating | – | – | 118 | 25 h | 12 cores | [21]

###  I-B Blocker Placement for GMD Mitigation

GICs are quasi-dc currents – produced by the low-frequency electric fields
that GMDs induce on the Earth’s surface – that appear in the conductive
infrastructure and flow into the high-voltage network through the neutrals of
transformers [22]. They are dependent on system characteristics (geographical
location of substations, transformer parameters), geomagnetic source fields
(amplitude, frequency content, spatial characteristics), and the Earth
conductivity structure (substation grounding resistance, influence on
geoelectric fields) [23, 24].

Placement of blocking devices on transformer neutrals can break the flow of
GICs and block them from entering/exiting the electrical grid at grounded
points. However, dc-current blocking devices do not prevent induced currents
on transmission lines themselves, and only drive the “effective” GICs
($I_{eff}$) to zero for auto-transformers and transformers with a single
grounded winding (typically delta-grounded-wye generator step-up
transformers); these devices have limited benefits for transformers with
multiple grounded windings, such as grounded-wye-grounded-wye or three-winding
transformers. Additionally, even with complete dc-blocker placement, it is not
guaranteed that $I_{eff}$ will be zero across the entire system, nor that the
placement problem will be feasible [24]; alternate mitigation methods, such as
series capacitor placement or line switching, should also be explored.

The literature on blocking device placement can generally be categorized using
two criteria: the problem formulation, and the method used to solve the
problem. Given the computational complexity associated with modeling the
physics of power flows, existing research has considered a wide range of
problem formulations, such as complete non-linear formulations [19, 20, 21,
14, 18], mathematical relaxations [15, 16, 7] and linearized approximations
[13, 11, 14, 17]. Some work has focused only on formulating the equivalent dc
network for calculating GIC flows [13, 11], others on the dc network with
reactive power loss ($q_{loss}$) calculations [15, 14] for coupling with ac
power flows, and finally, on formulations that feature the fully coupled dc
and ac networks [19, 16, 21, 7, 17]. Other potential additions include
transformer hot-spot constraints [21], multi-time-series extension [14, 4],
and various objectives; often the objective is to minimize total or maximum
$I_{eff}$ [13, 11], minimize total or maximum $q_{loss}$ [14], minimize
maximum transformer hot-spot temperature, minimize load-shedding, or minimize
dc-blocker placement cost. Solution methods may rely on linear solvers, local
solvers, or genetic algorithms/simulated annealing [13, 11].

While research have tried to tackle the scalability of this problem, current
state-of-the-art is limited to test cases with fewer than $2\times 10^{3}$
busses, far less than the size of real-world BES networks, such as the U.S.
Eastern Interconnect with $6\times 10^{4}$ busses. The rest of this paper
presents a complete optimization formulation for the GIC dc-blocker placement
problem – which includes features of practical network models such as line
charging or switched shunts – and highlights challenges in obtaining a
globally optimum solution in reasonable time. A reference implementation for
this formulation is provided in PowerModelsGMD.jl111https://github.com/lanl-
ansi/PowerModelsGMD.jl [4], which is benchmarked on a set of publicly
available test cases.

##  II Problem Formulation

This section describes the developed optimization-based problem formulation
for the strategic placement of blocking devices, based on full ac power flow
equations. Relative to conventional ac transmission network problems – such as
load flow, optimal power flow, optimal line switching, minimum load-shedding –
a positive-sequence representation of the problem is insufficient; the blocker
placement problem requires the addition of a separate quasi-dc representation
of the transmission network, which is coupled to the ac network through a set
of constraints associated with transformers. The quasi-dc representation
itself is fairly similar to the zero-sequence representation of the ac
network, with two major distinctions: inductances are neglected and assumed to
be short circuits, while capacitances are neglected and assumed to be open
circuits. It is an architectural choice for implementation whether the dc
network is to be modeled as a completely separate network, or added to the ac
network as a set of special dc busses (nodes) and lines (branches or edges).

In contrast to traditional ac problems, additional data is required for the
blocker placement problem: ground electric field magnitude and direction,
which are assumed to be known; geographical locations of transmission line
endpoints, or the full path of lines for increased fidelity; resistance to
remote earth of substation grounding grids; and transformer core
configurations. Furthermore, it is necessary to explicitly model generator
step-up transformers, which are frequently not included in published
transmission network datasets.

The end-to-end workflow of the GIC blocker placement problem is as follows:

  1. 1.

generation of the dc network given the ac network;

  2. 2.

line coupling calculations to determine the induced voltages along
transmission lines;

  3. 3.

calculations of the “effective” GICs ($I_{eff}$) on transformers;

  4. 4.

calculations of the transformer reactive power losses ($q_{loss}$); and

  5. 5.

solving the ac power flow problem.

Depending on the exact problem to be solved, some of these steps may be
coupled within the same optimization problem rather than being solved
sequentially.

###  II-A Objective Function

A number of objective functions for guiding the choice of GIC blocker
placement are reasonable, including:

  * •

minimizing the sum-squared of the $I_{eff}$ (subject to constraints on the
number of blockers placed);

  * •

minimizing the total cost of placing GIC dc-blockers (subject to power flow
constraints and a limit on the amount of allowed load-shedding); or

  * •

minimizing the amount of load-shedding (subject to constraints on the number
of blockers placed).

There is a clear trade-off between power system resilience (as expressed by
the total amount of load served during a GIC event) and the total number of
blockers placed throughout the system. Additional objective functions may also
be considered: the total of the weighted sum of $I_{eff}$ ($I_{eff}$ is a
weighted sum of the transformer winding currents, which is proportional to the
core magnetic flux density); the total of the weighted sum of the cost of
placing GIC dc-current blockers; or the total cost of operating generation
(however, as HILF events are considered, fuel cost of generation is arguably
not a concern).

The cost associated with operating generation, typical for optimal power flow
problems under nominal conditions, is as follows:

$\allowdisplaybreaks\small\mathrm{min.}\smashoperator[]{\sum_{i\in\mathcal{G}}^{}}\kappa_{i}^{0}+\kappa_{i}^{1}p_{i}^{g}+\kappa_{i}^{2}(p_{i}^{g})^{2}$
(1)

where $\mathcal{G}$ is the set of generators in the ac network; $p_{i}^{g}$ is
the power output of generator $i$; and
$\kappa_{i}^{0},\,\kappa_{i}^{1},\,\kappa_{i}^{2}$ are the cost coefficients
of generator $i$.

####  II-A1 Minimize cost of sum-squared GIC

The cost associated with minimizing GIC throughout the system is:

$\allowdisplaybreaks\small\mathrm{min.}\smashoperator[]{\sum_{i\in\mathcal{E}^{d}}^{}}|I_{i}|^{2}$
(2)

where $\mathcal{E}^{d}$ is the set of directed edges in the dc network,
oriented from dc node $i$ to dc node $j$; and $I_{i}$ is the $I_{eff}$ on edge
$i$ of that network.

####  II-A2 Minimize cost of blocker placement

The cost associated with placing GIC dc-current blockers is:

$\allowdisplaybreaks\small\mathrm{min.}\smashoperator[]{\sum_{i\in\mathcal{N}^{db}}^{}}\kappa_{i}^{b}z_{i}^{b}$
(3)

where $\mathcal{N}^{db}$ is the set of nodes in the dc network that correspond
to transformer neutrals and are candidates for placed blockers;
$\kappa_{i}^{b}$ is the cost of placing a blocker at node $i$; and
$z_{i}^{b}\in\\{0,1\\}$ is the binary placement variable for the blocker in
node $i$.

####  II-A3 Minimize amount of load shed

The cost associated with load-shedding is:

$\allowdisplaybreaks\small\mathrm{min.}\smashoperator[]{\sum_{i\in\mathcal{D}}^{}}\mathinner{\\!\left\lvert
p_{i}^{d}\right\rvert}\kappa_{i}^{d}z_{i}^{d}$ (4)

where $\mathcal{D}$ is the set of loads in the ac network; $p_{i}^{d}$ is the
dispatched power of load $i$; $\kappa_{i}^{d}$ is the cost of shedding load
$i$; and $z_{i}^{d}\in\\{0,1\\}$ is the load-shedding variable corresponding
to load $i$.

###  II-B AC Power Flow Constraints

The ac power flow constraints consist primarily of the set of equations that
represent the physics of power flow on an ac transmission network, and the set
of inequalities that represent the operational limits of the network.

####  II-B1 Slack bus constraint

To provide symmetry breaking, $\theta_{i}$ voltage angles are constrained to
be zero:

for $\forall\;i\in\mathcal{N}^{as}$

$\allowdisplaybreaks\footnotesize\theta_{i}=0$ (5)

where $\mathcal{N}^{as}$ is a set of slack nodes in the ac network.

####  II-B2 Nodal power balance constraint

for $\forall\;i\in{\mathcal{N}^{a}}$

$\displaystyle\smashoperator[]{\sum_{j\in{\mathcal{E}_{i}^{a+}}}^{}}p_{j}^{+}-\smashoperator[]{\sum_{j\in{\mathcal{E}_{i}^{a-}}}^{}}p_{j}^{-}=\sum_{j\in\mathcal{G}_{i}}p_{j}^{g}z_{j}^{g}-\sum_{j\in\mathcal{D}_{i}}p_{j}^{d}z_{j}^{d}-v_{i}^{2}\sum_{j\in\mathcal{S}_{i}}g_{j}^{s}z_{j}^{s}$
(6a)
$\displaystyle\smashoperator[]{\sum_{j\in{\mathcal{E}_{i}^{a+}}}^{}}q_{j}^{+}-\smashoperator[]{\sum_{j\in{\mathcal{E}_{i}^{a-}}}^{}}q_{j}^{-}=\sum_{j\in\mathcal{G}_{i}}q_{j}^{g}z_{j}^{g}-\sum_{j\in\mathcal{D}_{i}}q_{j}^{d}z_{j}^{d}-v_{i}^{2}\sum_{j\in\mathcal{S}_{i}}b_{j}^{s}z_{j}^{s}-\sum_{j\in\mathcal{L}_{i}}q_{j}^{loss}$
(6b)

where $\mathcal{N}^{a}$ is the set of nodes in the ac network;
$\mathcal{E}_{i}^{+}$ is the set of directed edges connected to node $i$ that
are oriented towards node $i$ in the ac network, while $\mathcal{E}_{i}^{-}$
is the set of directed edges connected to node $i$ that are oriented against
node $i$; $p_{j}^{+}$ and $q_{j}^{+}$ are the real and reactive power flowing
out the end of edge $j$, respectively, while $p_{j}^{-}$ and $q_{j}^{-}$ are
the real and reactive power flowing into the end of edge $j$, respectively;
$\mathcal{G}_{i}$ is the set of generators in the ac network connected to node
$i$, $p_{j}^{g}$ is the real power output of generator $j$ and $q_{j}^{g}$ is
the reactive power output of generator $j$, and $z_{j}^{g}\in\\{0,1\\}$ is the
binary on-off variable associated with generator $j$; $\mathcal{D}_{i}$ is the
set of loads in the ac network connected to node $i$, $p_{j}^{d}$ is the real
power demand of load $j$ and $q_{j}^{d}$ is the reactive power demand of load
$j$, and $z_{j}^{d}\in\\{0,1\\}$ is the binary on-off variable associated with
load $j$; $v_{i}$ is the voltage magnitude at node $i$, while $q_{j}^{loss}$
is the reactive power losses associated with transformer primary windings
connected to node $i$; $\mathcal{S}_{i}$ is the set of shunts in the ac
network connected to node $i$, $g_{j}^{s}$ is the real portion of the
admittance of shunt $j$, while $b_{j}^{s}$ is the reactive portion of the
admittance of shunt $j$, and $z_{j}^{s}\in\\{0,1\\}$ is the binary on-off
variable associated with shunt $j$; and $\mathcal{L}_{i}$ is the set of
$q_{loss}$ pseudo-loads – associated with the GMD-induced reactive power loss
for transformers – in the ac network connected to node $i$, while
$q_{j}^{loss}$ is the GMD-induced reactive power loss associated with
transformer $j$ connected to node $i$.

####  II-B3 Branch flow constraint

Constraints for edge $i$, connecting nodes $j$ and $k$:

for $\forall\;i\in{\mathcal{E}^{a}}$

$\displaystyle p_{i}^{-}$
$\displaystyle=-g_{i}v_{j}v_{k}\cos(\theta_{j}-\theta_{k})+b_{i}v_{j}v_{k}\sin(\theta_{j}-\theta_{k})+g_{i}v_{j}^{2}$
(7a) $\displaystyle q_{i}^{-}$
$\displaystyle=-b_{i}v_{j}v_{k}\cos(\theta_{j}-\theta_{k})+g_{i}v_{j}v_{k}\sin(\theta_{j}-\theta_{k})+\left(b_{i}+\frac{1}{2}b_{i}^{sh}\right)v_{j}^{2}$
(7b) $\displaystyle p_{i}^{+}$
$\displaystyle=g_{i}v_{j}v_{k}\cos(\theta_{j}-\theta_{k})-b_{i}v_{j}v_{k}\sin(\theta_{j}-\theta_{k})-g_{i}v_{k}^{2}$
(7c) $\displaystyle q_{i}^{+}$
$\displaystyle=b_{i}v_{j}v_{k}\cos(\theta_{j}-\theta_{k})-g_{i}v_{j}v_{k}\sin(\theta_{j}-\theta_{k})-\left(b_{i}+\frac{1}{2}b_{i}^{sh}\right)v_{k}^{2}$
(7d)

where $\mathcal{E}^{a}$ is the set of directed edges in the ac network;
impedances are represented in rectangular coordinates, where $g_{i}$ and
$b_{i}$ are the real and reactive portions of the series admittance of edge
$i$, respectively, while $b_{i}^{sh}$ is the shunt admittance of edge $i$
resulting from line capacitance; and voltages are represented in polar
coordinates, where $v_{j}$ is the magnitude of the voltage at node $j$, and
$\theta_{j}$ is the angle of the voltage at node $j$.

###  II-C AC Operational Limit Constraints

Voltage limits at each node are:

for $\forall i\in{\mathcal{N}^{a}}$

$\allowdisplaybreaks\footnotesize\underline{v}_{i}\leq
v_{i}\leq\overline{v}_{i}$ (8)

where the lower and upper bounds of the voltage magnitude $v_{i}$ at node
$i\in\mathcal{N}^{a}$ are set by the range
$[\underline{v}_{i},\,\overline{v}_{i}]$.

Thermal limits of edges in both directions are:

for $\forall i\in{\mathcal{E}^{a}}$

$\displaystyle(p_{i}^{-})^{2}+(q_{i}^{-})^{2}\leq\overline{s}_{i}^{2}$ (9a)
$\displaystyle(p_{i}^{+})^{2}+(q_{i}^{+})^{2}\leq\overline{s}_{i}^{2}$ (9b)

where $\overline{s}_{i}$ is the thermal limit of edge $i$.

Voltage stability limits on the phase angle difference on edge $i$, between
nodes $j$ and $k$, are enforced:

for $\forall i\in{\mathcal{E}^{a}}$

$\allowdisplaybreaks\footnotesize\underline{\theta}_{i}\leq\theta_{j}-\theta_{k}\leq\overline{\theta}_{i}$
(10)

where $\underline{\theta}_{i}$ and $\overline{\theta}_{i}$ represent the upper
and lower bounds on the phase angle difference.

The capacity of power generation:

for $\forall i\in\mathcal{G}$

$\displaystyle\underline{p}_{i}^{g}$ $\displaystyle\leq
p_{i}^{g}\leq\overline{p}_{i}^{g}$ (11a) $\displaystyle\underline{q}_{i}^{g}$
$\displaystyle\leq q_{i}^{g}\leq\overline{q}_{i}^{g}$ (11b)

where $\mathcal{G}$ is a set of generators in the ac network connected to node
$i\in\mathcal{N}^{a}$; $\underline{p}_{i}^{g}$ and $\overline{p}_{i}^{g}$
represent the lower and upper bounds on real power for generator $i$; and
$\underline{q}_{i}^{g}$ and $\overline{q}_{i}^{g}$ represent the lower and
upper bounds on reactive power for generator $i$.

###  II-D Quasi-DC Voltage and Current Constraints

To formulate the nodal current balance constraint, an edge in the dc circuit
is linked to an edge in the ac circuit:

for $\forall i\in{\mathcal{N}^{d}}$

$\allowdisplaybreaks\footnotesize\smashoperator[]{\sum_{j\in\mathcal{E}^{d+}_{i}}^{}}I_{j}-\smashoperator[]{\sum_{j\in\mathcal{E}^{d-}_{i}}^{}}I_{j}=(1-z_{i}^{b})a_{i}^{s}V_{i}$
(12)

where $\mathcal{N}^{d}$ is the set of nodes in the dc network;
$\mathcal{E}_{i}^{d+}$ is the set of directed edges connected to node $i$ in
the dc network that are oriented towards node $i$, while
$\mathcal{E}_{i}^{d-}$ is the set of directed edges connected to node $i$ that
are oriented against node $i$; $I_{j}$ is the GIC flow through dc branch $j$;
$z_{i}^{b}\in\\{0,1\\}$ is the binary placement variable for the blocker for
node $i$, which when placed forces GIC flow into a transformer neutral to 0;
$V_{i}$ is the quasi-dc voltage at node $i$; and $a_{i}^{s}$ is the resistance
to remote earth at $i$. This constraint introduces a bilinear term with both
$z_{i}^{b}$ and $V_{i}$.

Ohm’s law constraint:

for $\forall i\in\mathcal{E}^{d}$

$\allowdisplaybreaks\footnotesize
I_{i}=z_{l}a_{i}(V_{j}-V_{k}+\mathbf{E}_{i})$ (13)

where $l=A_{i}$ is the ac edge associated with dc edge $i$; $z_{l}$ is the
status of ac edge $l$; and $a_{i}$ is the conductance of dc edge $i$.

GIC on a transmission line is determined by the induced voltage source
$\mathbf{E}_{i}$, an unknown real-valued quantity that depends on the
displacement of the line in addition to the magnitude and orientation of the
ground electric field along that displacement. While this developed
formulation does not address how to handle this uncertainty, existing
literature offers solutions; for example, scenario-based approaches or robust
optimization approaches could be included [25, 9].

The notion of $I_{eff}$ is to calculate a weighted sum of currents through
transformer windings that is proportional to the transformer core flux, and
therefore to the GIC-induced $q_{loss}$ of the transformer. This varies for
different transformer configurations [26]:

for $\forall i\in{\mathcal{E}^{a}}$

$\allowdisplaybreaks\footnotesize\widetilde{I}_{i}=\begin{cases}I_{i^{H}}&\text{if
$i\in\mathcal{E}^{\Delta}$}\vspace{.1cm}\\\
{\frac{\alpha_{i}I_{i^{H}}+I_{i^{L}}}{\alpha_{i}}}&\text{if
$i\in\mathcal{E}^{y}$}\vspace{.1cm}\\\
{\frac{\alpha_{i}I_{i^{S}}+I_{i^{C}}}{\alpha_{i}+1}}&\text{if
$i\in\mathcal{E}^{\infty}$}\vspace{.1cm}\\\
I_{i^{H}}+\frac{I_{i^{L}}}{\alpha_{i}}+\frac{I_{i^{T}}}{\beta_{i}}&\text{if
$i\in\mathcal{E}^{3}$}\vspace{.1cm}\\\ 0&\text{otherwise}\end{cases}$ (14)

where $\mathcal{E}^{\Delta}$ is the set of delta-delta connected transformers,
$\mathcal{E}^{y}$ is the set of grounded-wye-delta connected transformers,
$\mathcal{E}^{\infty}$ is the set of autotransformers, while $\mathcal{E}^{3}$
is the set of three-winding transformers;
$\alpha_{i}=V_{i^{H}}^{b}/V_{i^{L}}^{b}$ is the transformer turns ratio for
the primary and secondary windings, while
$\beta_{i}=V_{i^{H}}^{b}/V_{i^{3}}^{b}$ is the transformer turns ratio for the
primary and tertiary windings for the case of three-winding transformers;
$V_{i^{H}}^{b}$, $V_{i^{L}}^{b}$, and $V_{i^{L}}^{3}$ are the nominal voltages
for the high-side node, low-side node, and tertiary node corresponding to
transformer $i$, respectively; $\widetilde{I}_{i}$ is the $I_{eff}$ current on
the dc edge; $I_{i^{H}}$ is the high-side winding current and $I_{i^{L}}$ is
the low-side for two- or three-winding transformers; $I_{i^{S}}$ is the series
winding current and $I_{i^{C}}$ is the common winding current for auto-
transformers; and
$i^{H}=H_{i}^{d},\,i^{L}=L_{i}^{d},\,i^{S}=S_{i}^{d},\,i^{C}=C_{i}^{d},\,i^{3}=T_{i}^{d}$
are set index variables for transformer terminal nodes in the dc network.

The magnitude of $I_{eff}$:

for $\forall i\in\mathcal{E}^{\tau}$

$\allowdisplaybreaks\footnotesize\overline{I}_{i}=\left|{\widetilde{I}_{i}}\right|$
(15)

Last, the $q_{loss}$ increase for a transformer resulting from GIC-induced
half-cycle saturation:

for $\forall i\in\mathcal{E}^{\tau}$

$\allowdisplaybreaks\footnotesize\small
q^{loss}_{j}=\sqrt{\frac{2}{3}}\cdot\frac{S_{i}^{b}}{V_{j}^{b}}\cdot
K_{i}\overline{I}_{i}^{d}v_{j}$ (16)

where $j=H_{i}$ is the high-side node associated with transformer $i$;
$q_{j}^{loss}$ is the GIC-induced reactive power loss at the high-side node
$j$; $S_{i}^{b}$ is the base power of transformer $i$; $V_{j}^{b}$ is the
rated voltage at the high-side node $j$; $K_{i}$ is the GIC reactive power
loss coefficient of transformer $i$; $\overline{I}_{i}^{d}$ is the absolute
value of $I_{eff}$; and $v_{j}$ is the high-side node voltage. This constraint
also introduces a bilinear term in $\overline{I}_{i}^{d}$ and $v_{j}$. While
the ac node voltage has relatively tight bounds, typically in the range
$[0.9,1.1]$, this is not the case for $I_{eff}$, which is typically in the
range $[0,1000]$ A.

###  II-E Blocker Placement Constraints

To handle the trade-off between the total amount of load served (load lost)
during a severe contingency and the total number (cost) of GIC dc-blockers
placed throughout the system, additional constraints are introduced.

####  II-E1 Cost of dc-current blockers

To restrict the total cost of blocker placement under a specified threshold
budget:

$\allowdisplaybreaks\small\smashoperator[]{\sum_{i\in\mathcal{N}_{b}^{d}}^{}}\kappa_{i}^{b}z_{i}^{b}\leq\overline{c}_{b}$
(17)

where $\kappa_{i}^{b}$ is the cost of placing a GIC dc-blocker at node $i$;
and $\overline{c}_{b}$ is the total budget for blocker placement.

This constraint is intended to be used in conjunction with a minimum load-
shedding objective.

####  II-E2 Number of dc-blockers

To restrict the total number of GIC dc-blockers placed in the system under a
specified threshold:

$\allowdisplaybreaks\small\smashoperator[]{\sum_{i\in\mathcal{N}^{db}}^{}}z_{i}^{b}=n^{b}$
(18)

where $n^{b}$ is the total number of blockers to be placed.

####  II-E3 Cost of load-shedding

To restrict the total cost of load-shedding under a specified threshold, given
that the system has loads with varying priorities:

$\allowdisplaybreaks\small\smashoperator[]{\sum_{i\in\mathcal{D}}^{}}\mathinner{\\!\left\lvert
p_{i}^{d}\right\rvert}\kappa_{i}^{d}z_{i}^{d}\leq\overline{c}_{d}$ (19)

where $\overline{c}_{d}$ is the maximum allowable cost of load-shedding.

####  II-E4 Fraction of load served

In the case, when the cost of blocker placement (Subsubsection II-E1) rather
than load-shedding (Subsubsection II-E3) is the quantity to be minimized, a
supporting constraint may be the fraction of load served; this must be greater
than a specified threshold:

$\allowdisplaybreaks\footnotesize\small\smashoperator[]{\sum_{i\in\mathcal{D}}^{}}\mathinner{\\!\left\lvert
p_{i}^{d}\right\rvert}z_{i}^{d}\leq\rho_{d}\smashoperator[]{\sum_{i\in\mathcal{D}}^{}}\mathinner{\\!\left\lvert
p_{i}^{d}\right\rvert}$ (20)

where $\rho_{d}$ is the minimum fraction of load that must be served.

##  III Reference Implementation and Benchmarking

The above introduced problem formulation was implemented in PowerModelsGMD.jl
v5.0, an open-source Julia tool within the InfrastructureModels software
ecosystem [27]. For the optimal placement of GIC dc-current blocking devices,
the selected objective was to minimize the number of blockers placed –
therefore, minimizing the cost of blocker placement, as described in
Subsubsection II-E1 – subject to Constraint (17); additionally, the fraction
of load served was used as supporting constraint – load-shedding upper bounded
was set to 15% of total load, with all loads weighted equally, as described in
Subsubsection II-E4 – subject to Constraint (20). This implementation is
presented as a reference and is not intended for practical use on realistic
network sizes; the intention is to quantify the computational difficulty of
the blocker placement problem before heuristics are applied to scale to larger
network sizes.

The implementation was benchmarked on a computer running Windows 10 with an
Intel Xeon Silver 4112 CPU with 2.60 GHz and 64 GB RAM, using Juniper [28] as
a solver, and applied to three publicly available test cases of varying sizes
(from 4 to 150 busses) [29, 30, 31]. Results are summarized in Table II:
listing the optimal number of blockers placed (versus the total number of
possible locations where blockers can be placed), the percentage of load met,
the unit cost of the placed blocker(s), and the required computation time.

TABLE II: Summary of GIC Blocker Placement Methods Case | Busses | Blockers | Load Met | Cost | Runtime
---|---|---|---|---|---
|  |  |  |  | (CPU Seconds)
B4GIC | 4 | 1/2 | 100% | 1.0 | 12.9
EPRI21 | 20 | 1/8 | 81.8% | 1.0 | 9.2
UIUC150 | 169 | 1/98 | 100% | 1.0 | 1759

##  IV Conclusion

The threat of GMDs pose a risk to the BES, and addressing the challenges is
difficult because of the lack of empirical data. This means that simulation
and optimization methods are important for analysis and impact mitigation.

This paper presented a complete base problem formulation for the optimal
placement of transformer neutral (dc-current) blocking devices, an effective
technology to reduce GICs without reducing the power-carrying capacity of the
grid. Although this formulation is not tractable – as it is an NP-hard problem
with mixed-integer nonlinear, nonconvex constraints – and is formulated as an
ac polar formulation – other formulations, such as semidefinite or second-
order cone programming [32], may be used as well – this is an important first
step towards the development of a scalable, optimization-based method to guide
the placement of blocking devices.

The outlined problem formulation was implemented in PowerModelsGMD.jl and
benchmarked on three test cases of varying sizes, providing comparison for
future work. As the benchmarking experiments show, obtaining a globally
optimum solution in reasonable time is increasingly difficult as the network
size increases. Relaxations and approximations – for the GIC flow and
$q_{loss}$ coupling constraints (including the absolute value, the $q_{loss}$
constraint, and the dc current nodal balance constraints) that all include
bilinear terms – need to be explored an applied to improve the computational
performance and scale the base problem to realistic network sizes. This is the
subject of ongoing research work.

## References

  * [1] National Academies of Sciences, Engineering, and Medicine.  Enhancing the Resilience of the Nation’s Electricity System.  The National Academies Press. Washington DC, 2017.
  * [2] NERC and U.S. Department of Energy.  High-Impact, Low-Frequency Event Risk to the North American Bulk Power System.  Technical report, North American Electric Reliability Corporation, Jun. 2010.
  * [3] A. Mate et al.  Relaxation Based Modeling of GMD Induced Cascading Failures in PowerModelsGMD.jl.  In Proceedings of the 2021 IEEE/PES 53rd North American Power Symposium, 2021.
  * [4] A. Mate et al.  Analyzing and Mitigating the Impacts of GMD and EMP Events on the Electrical Grid with PowerModelsGMD.jl, Jan. 2021.  https://arxiv.org/abs/2101.05042.
  * [5] G. P. Juvekar et al.  GIC-Inclusive State Estimator for Power System Awareness During Geomagnetic Disturbance Events.  IEEE Transactions on Power Systems, 36(4):2966–2974, Jul. 2021\.
  * [6] T. J. Overbye et al.  Towards Developing Implementable High Altitude Electromagnetic Pulse E3 Mitigation Strategies for Large-Scale Electric Grids.  In Proceedings of the 2022 IEEE Texas Power and Energy Conference, 2022.
  * [7] M. Kazerooni et al.  Mitigation of Geomagnetically Induced Currents Using Corrective Line Switching.  IEEE Transactions on Power Systems, 33(3):2563–2571, May 2018\.
  * [8] M. Lu et al.  Optimal Transmission Line Switching Under Geomagnetic Disturbances.  IEEE Transactions on Power Systems, 33(3):2539–2550, May 2018\.
  * [9] M. Ryu et al.  Mitigating the Impacts of Uncertain Geomagnetic Disturbances on Electric Grids: A Distributionally Robust Optimization Approach.  IEEE Transactions on Power Systems, 37(6):4258–4269, Nov. 2022\.
  * [10] R. Zhang et al.  Machine Learning-Aided Enhancement of Power Grid Resilience to Electromagnetic Pulse Strikes.  In Proceedings of the 2022 North American Power Symposium, 2022\.
  * [11] P.-H. Yang et al.  Optimal Placement of Grounding Small Resistance in Neutral Point for Restraining Voltage Fluctuation in Power Grid Caused by Geomagnetic Storm.  IET Generation, Transmission & Distribution, 13(8):1456–1465, Apr. 2019.
  * [12] R. Horton.  Magnetohydrodynamic Electromagnetic Pulse Assessment of the Continental U.S. Electric Grid.  Technical Report 3002009001, Electric Power Research Institute, 2017.
  * [13] X. Ning et al.  Research on Optimal Placement for GIC Mitigation with Blocking Device.  In IOP Conference Series: Materials Science and Engineering, 2019\.
  * [14] Y. Liang et al.  Optimal Blocking Device Placement for Geomagnetic Disturbance Mitigation.  IEEE Transactions on Power Delivery, 34(6):2219–2231, Dec. 2019\.
  * [15] H. Zhu et al.  Blocking Device Placement for Mitigating the Effects of Geomagnetically Induced Currents.  IEEE Transactions on Power Systems, 30(4):2081–2089, Jul. 2014\.
  * [16] H. Zhu et al.  Blocking Device Placement for Mitigating the Effects of Geomagnetically Induced Currents.  In Proceedings of the 2016 IEEE Power and Energy Society General Meeting, 2016.
  * [17] W. Shengquan et al.  Optimal Blocking Devices Placement for Geomagnetic Disturbance Mitigation Based on Sensitivity of Induced Geoelectric Fields.  IEEE Access, 10:132814–132821, Nov. 2022.
  * [18] A. Lamichhane.  Optimal Mitigation of Geomagnetically Induced Current Effects in Power Systems Considering Transformer Thermal Limits.  Masters Thesis, York University, Dec. 2022.
  * [19] A. H. Etemadi et al.  Optimal Placement of GIC Blocking Devices for Geomagnetic Disturbance Mitigation.  IEEE Transactions on Power Systems, 29(6):2753–2762, Mar. 2014\.
  * [20] Y. Liang et al.  Optimal Blocker Placement for Mitigating the Effects of Geomagnetic Induced Currents Using Branch and Cut Algorithm.  In Proceedings of the 2015 North American Power Symposium, 2015\.
  * [21] A. Rezaei-Zare et al.  Optimal Placement of GIC Blocking Devices Considering Equipment Thermal Limits and Power System Operation Constraints.  IEEE Transactions on Power Delivery, 33(1):200–208, Feb. 2018\.
  * [22] R. Pirjola.  Geomagnetically Induced Currents During Magnetic Storms.  IEEE Plasma Sci., 28(6):1867–1873, Dec. 2000.
  * [23] K. S. Shetye et al.  Modeling and Analysis of GMD Effects on Power Systems: An Overview of the Impact on Large-Scale Power Systems.  IEEE Electrific. Mag., 3(4):13–21, Dec. 2015.
  * [24] D. H. Boteler et al.  Modeling Geomagnetically Induced Currents.  Space Weather, 15(1):258–276, Dec. 2016.
  * [25] A. Barnes et al.  Resilient Design of Large-Scale Distribution Feeders with Networked Microgrids.  Electric Power Systems Research, 171(0):150–157, Jun. 2019.
  * [26] T. J. Overbye et al.  Integration of Geomagnetic Disturbance Modeling into the Power Flow.  In Proceedings of the 2012 North American Power Symposium, 2012\.
  * [27] R. Bent et al.  InfrastructureModels: Composable Multi-Infrastructure Optimization in Julia.  2023\.  Under Review.
  * [28] O. Kröger et al.  Juniper: An Open-Source Nonlinear Branch-and-Bound Solver in Julia.  In Proceedings of the Integration of Constraint Programming, Artificial Intelligence, and Operations Research, pages 377–386, 2018.
  * [29] PowerWorld Corporation.  Quick-Start Guides.  https://www.powerworld.com/training/quick-start-guides, May 2023\.
  * [30] R. Horton et al.  A Test Case for the Calculation of Geomagnetically Induced Currents.  IEEE Transactions on Power Delivery, 27(4):2368–2373, Oct. 2012\.
  * [31] A. Birchfield et al.  Statistical Considerations in the Creation of Realistic Synthetic Power Grids for Geomagnetic Disturbance Studies.  IEEE Transactions on Power Systems, 32(2):1502–1510, Mar. 2017\.
  * [32] C. Coffrin et al.  The QC Relaxation: Theoretical and Computational Results on Optimal Power Flow, Jul. 2015.  http://arxiv.org/abs/1502.07847.