1

Optimal Transmission Line Switching under
Geomagnetic Disturbances

arXiv:1701.01469v4 [cs.SY] 5 Oct 2017

Mowen Lu∗ , Harsha Nagarajan† , Emre Yamangil† , Russell Bent‡ , Scott Backhaus‡ , Arthur Barnes‡

Abstract—In recent years, there have been increasing concerns
about how geomagnetic disturbances (GMDs) impact electrical
power systems. Geomagnetically-induced currents (GICs) can
saturate transformers, induce hot spot heating and increase reactive power losses. These effects can potentially cause catastrophic
damage to transformers and severely impact the ability of a
power system to deliver power. To address this problem, we
develop a model of GIC impacts to power systems that includes
1) GIC thermal capacity of transformers as a function of normal
Alternating Current (AC) and 2) reactive power losses as a
function of GIC. We use this model to derive an optimization
problem that protects power systems from GIC impacts through
line switching, generator redispatch, and load shedding. We
employ state-of-the-art convex relaxations of AC power flow
equations to lower bound the objective. We demonstrate the
approach on a modified RTS96 system and the UIUC 150-bus
system and show that line switching is an effective means to
mitigate GIC impacts. We also provide a sensitivity analysis of
optimal switching decisions with respect to GMD direction.
Index Terms—GMD, transmission line switching, convex relaxations, ACOPF, GIC.

I. I NTRODUCTION
OLAR flares and coronal mass ejections drive geomagnetic disturbances (GMD) that lead to changes in
the Earth’s magnetic field, which then create geo-electric
fields. These low-frequency geo-electric fields induce quasiDC currents, also known as Geomagnetically-Induced Currents (GICs), in grounded sections of power system networks
[1]–[3]. The GIC are superimposed on the usual alternating
currents (AC) and bias the AC such that the maximum currents
are increased. In many power system components, this bias is
not a major concern, however, in transformers, this bias can
lead to half-cycle saturation of transformers and the loss of
magnetic flux to regions outside of the transformer core. The
energy stored in the stray flux increases the reactive power
consumption of the transformer, which can affect system voltages. The stray flux also drives eddy currents that can cause
excessive transformer heating leading to reduced transformer
life or, potentially, immediate damage [4].
The potential impacts of GMD to transformers in the bulk
electric power system have led the United States government to
increase the understanding of and mitigate the impacts of such
events [5], [6]. To mitigate the potential risks introduced by
GIC to power systems, the electric power industry has actively
improved GIC modeling and GIC monitoring [7]–[12]. These
models have been used to conduct risk analysis [10], [13]

S

∗ Graduate Student, Department of Industrial Engineering, Clemson University, SC, USA
† Center for Nonlinear Studies, Los Alamos National Laboratory, NM, USA
‡ Information Systems and Modeling (A-1 division), Los Alamos National
Laboratory, NM, USA

that investigate the sensitivity of transformer reactive power
losses due to GIC and concluded that risk and risk mitigation
warrants further study.
One focus in the recent literature has been on mitigating
the effect of GIC on transformer reactive power consumption
and subsequent drops in system voltages and potential voltage
collapse. One approach to mitigation is the installation of
DC-current blocking devices to keep the GIC from entering
through transformer neutrals [14], however, these devices
are expensive, with costs for a single unit close to $500K
[15]–[17]. In an attempt to minimize the projected cost of
mitigation, optimization-based methods have been developed
to guide the siting of these blocking devices. Instead of
performing a full power systems analysis that includes the
AC, GIC and full AC power flow equations, these papers
have primarily focused on minimizing induced reactive losses
independent of the normal AC currents. The intuition of these
surrogate models is that small amounts of reactive losses imply
small voltage impacts and, presumably, a secure power system.
Beyond voltage effects, the literature on risk mitigation
associated with transformer heating is relatively sparse. Existing studies focus on assessing transformer susceptibility
to GIC effects [18] and formulating the thermal response of
transformer cores to different levels of GMDs [19]. However,
this approach was strictly a screening study and did not
recommend methods for mitigation.
The work discussed above is a very important start, but it
leaves a number of open questions, which we address in this
manuscript. First, the installation of blocking devices is very
expensive and cost may pose a barrier to adoption. Instead,
we focus on developing a GIC-aware optimal power flow
(OPF) model that uses existing controls such as generator
dispatch, load shedding, and line switching to mitigate the
risks of GIC impacts. Second, we incorporate the AC physics
of power flow into the GIC-aware OPF because these physics
play an important role in the impacts associated with GIC. For
example, while minimizing reactive losses may imply small
voltage problems across the whole system, these models focus
on total losses and can miss relatively large voltage problems
in a small part of a system. More importantly, models of
hotspot thermal heating inherently depend on both GIC and
AC.
The setting considered in this manuscript is very challenging. It combines transformer reactive losses, transformer
heating, and full AC power flow into an optimization-based
operational mitigation setting with line switching. By itself,
optimal transmission line switching (OTS) with AC power
flow physics is a mathematically challenging problem that

2

includes nonlinearities, nonconvexities and discrete variables.
Existing solution methodologies designed for OTS heavily rely
on tight convex relaxations and advanced discrete optimization
techniques. In recent literature, various convex relaxations and
disjunctive representations have been developed. These include
second-order-conic (SOC) relaxations [20], quadratically constrained (QC) relaxations [21] and Semi-definite programming
relaxations [22]. In the context of transmission expansion
planning applications, the QC relaxations have been effective
[23], [24] and we use this model here. Despite these recent
advances in optimization methods for OTS, global methods
still cannot scale to systems with 500 nodes.
The main contributions of this paper are the formulation
and initial algorithmic solution approaches to an operational
decision support tool that incorporates:
1) A model of transformer heating as a response to AC and
GIC-induced DC,
2) A realistic, coupled model of convex, relaxed AC power
flows with GIC effects and an algorithm to recover good
feasible solutions quickly, and
3) An optimization problem that protects the system from
reactive losses and thermal heating induced by GIC.
II. GIC MODELING AND ACOTS FORMULATION
N OMENCLATURE
Sets
N a , N d , N o set of nodes in the AC and DC circuit, respectively, where N o = N a ∩ N d
g
a
N ⊆ N set of nodes with exactly one generator
I ⊆ Nd
set of substation neutrals
E a , E d , E set of edges in the AC and DC circuit, respectively,
where E = E a ∪ E d
o
a
E ⊆E
set of transmission lines
E g ⊆ E a set of edges eij such that either i or j ∈ N g
E τ ⊆ E d set of transformer edges used to model the high
voltage primary windings of GSU transformers
and the common windings of autotransformers in
the DC circuit.
Ei+ ⊆ E
set of outgoing edges connected to AC/DC node i
Ei− ⊆ E
set of incoming edges connected to AC/DC node
i
Ei
set of all edges connected to AC/DC node i, where
Ei = Ei+ ∪ Ei−
τ
τ
Ei ⊆ E
set of DC edges used to compute Qloss
(as dei
scribed later) for node i
Parameters
c0i , c1i , c2i generation cost coefficients of generator i ∈ N g
ηe0 , ηe1 , ηe2 coefficients of the thermal limit curve of transformer line e ∈ E τ
µ
cost of load shedding
am
admittance of the grounding line at bus m ∈ I, 0
if bus m 6∈ I
ae
DC admittance of edge e ∈ E d
Je
induced current by GMD on line e ∈ E d
re , xe
resistance and reactance of line e ∈ E a
ge , be
conductance and susceptance of line e ∈ E a
gi , bi
shunt conductance and susceptance at bus i ∈ N a

dpi , dqi
bce
se
θ
θM
a
Ie
ke
V i, V i
gpi , gpi
gq i , gq i
φ
Vd
LN ,LE

real and reactive power demand at bus i ∈ N a
line charging susceptance of line e ∈ E a
apparent power limit on line e ∈ E a
phase angle difference limit
Big-M parameter given by |E a |θ
AC current flow limit on line e ∈ E a
loss factor of transformer line e ∈ E τ
AC voltage limits at bus i ∈ N a
real power generation limits at generator i ∈ G
reactive power generation limits at generator i ∈ G
the angle of the geo-electric field relative to east
GMD induced voltage source
the north and east components of the displacement
of each transmission line, respectively
EN ,EE
strength of the north and east geo-electric field,
respectively
Binary Variables
ze
1 if line e ∈ E a is switched on; 0 otherwise
Continuous Variables
θi
phase angle at bus i ∈ N a
Vi
voltage magnitude at bus i ∈ N a
d
Vi
induced DC voltage magnitude at bus i ∈ N d
le
AC magnitude squared on line e ∈ E a
d
GIC flow on transformer line e ∈ E τ
Ie
a
Iee
AC magnitude on line e ∈ E a
the effective GIC on transformer line e ∈ E τ
Ieed
Qloss
GIC-induced
reactive power loss at bus i ∈ N a
i
pij , qij
real and reactive power flow on line eij ∈ E a , as
measured at node i
fip , fiq
real and reactive power generated at bus i ∈ N a
lip , liq
real and reactive power shed at bus i ∈ N a
Each edge, eij ∈ E, is given an arbitrary orientation
from bus i to bus j. We omit the ij subscript when the
−
orientation is not relevant. For e ∈ E d , we use notation →
e
to denote the associated AC edge of e. This is a one-to-one
mapping for transmission lines and a many-to-one mapping
for transformers (discussed later).

A. GIC Modeling
Je calculation The computation of transformer hot spot
heating and GIC-induced reactive power losses depends on
the induced current sources (Je ) on each power line e ∈ E d
in the network, which itself depends on the strength and
direction of the geo-electric field associated with the GMD.
These relationships are modeled in Eq.(1)
Je = ae V d = ae

I

~ · d~l,
E

(1)

~ is the geo-electric field at the location of the transwhere, E
mission line, and d~l is the incremental line segment length,
including direction [12]. In practice, the actual geo-electric
field varies with time and geographical locations. Using a
common assumption that the north and east components of the
geo-electric field are constant in the geographical area of the

3

transmission line [11], [12], [16]1 , Je is calculated as (superand sub-scripts indicating edges are omitted):
J = aV d = a(EN LN + EE LE ) = a|E|(sin(φ)LN + cos(φ)LE ), (2)

where LN , LE , EN , EE and φ are as described in the nomenclature (see Appendix I of [12]). Given their short length,
generally Je = 0 for transformers, i.e. e ∈ E τ .
Transformer modeling The two most common transformers
in electrical transmission systems subject to GIC are network
transformers and generator step-up (GSU) transformers. Network transformers are generally located relatively far from
generators and transform voltage between different sections of
the transmission system. In contrast, GSUs connect the output
terminals of generators to the transmission network. Many
IEEE transmission reliability test networks explicitly model
network transformers, but generally do not model GSUs.
However, GSUs and the neutral leg ground points they provide
are critical when modeling GICs and methods to mitigate the
impact of GICs.
In this manuscript, we modify the IEEE RTS test network
by adding a GSU transformer between each generator and
its injection bus (see Fig. 1(a)). Consistent with common
engineering practice, we assume that each GSU is grounded
on its high voltage side that connects to the transmission
network. We also model the switching of the circuit breaker
between the high side of the GSU and the transmission
network using a binary variable that allows the GSU to be
isolated from the network and the quasi-DC GIC to protect the
GSU. This switching is performed if the generator output is
zero. Although the IEEE test networks include network transformers, transformer type and grounding data are typically
not provided. In this manuscript, we assume that all network
transformers are auto-transformers and each transformer has a
single neutral ground on the low voltage side.
Figure 1 includes examples of both GSU and network auto
transformer modeling. Figure 1(a) shows a four-bus section
of the transmission system with a single network transformer
(Tjk ) and two GSU transformers (Tia , Tib ) independently connecting two generators (Gai , Gbi ) to the same injection bus
i. In the simplified AC network of Fig. 1(b), bus ia and ib
model output terminals of generator Gai and Gbi , respectively.
Each GSU transformer Tia (Tja ) is reduced to a (single)
series impedance iia (iib ) with a circuit breaker. The network
transformer Tjk is reduced to a (single) series impedance (jk)
with a circuit breaker. Under this transformation, the number
of buses and lines of the AC network grow to |N o | + |N g |
and |E o | + |N g |, respectively, where |N o | and |E o | model the
original set of buses and edges in the network. Fig. 1(c) shows
an equivalent single-phase DC circuit of the example system
in nodal form. In this figure, m and n model the neutral point
of substation A and B, respectively. RTia and RTia denote
the resistance of the primary HV winding of Tia and Tib ,
respectively.2 RCTjk and RSTjk represent the resistance of
the common and series windings of Tjk , respectively. For
grounded GSU transformers, the effective GIC flows through
1 Our model does not depend on this assumption. It only depends on J as
e
an input parameter.
2 R corresponds to the inverse DC admittance, i.e. a = 1 ,
R

the primary HV winding. For example, in Fig. 1(c), the
effective GIC of GSU transformer Tib is IeTd b which is the
i
GIC flow from bus i to m on Tib [13], [25]. For an autotransformer, the effective GIC is derived from the GIC flows
through both the series and common windings as shown in
Fig. 1(c), i.e.,
αA IS + IC
αIH + IL
=
Ied =
α
αA + 1
where α is the turns ratio and αA = αS /αC = α − 1 (Eq.(14)
and Eq.(15) in [25]). In this manuscript, we assume the turn
ratios of all auto-transformers are one. As a result, Ied = |IH +
IL | = |IC |, i.e., the effective GIC is the GIC flow through the
common winding.3 Thus, in the four-bus network, the effective
GIC of autotransformer Tjk is the GIC flow on line (k, n) of
d
.
Fig. 1(c), i.e., Iekn
SUB A

Tia

SUB B

Gai

i

ia

j
k

Tib

Gai

s

Tjk

i
i

Gbi

j

k

s

b

Gbi

(a) 4-bus system
Jij

i
i

m

Jks

IS

j IH
Rij

RT a

(b) Equivalent AC network

RT b
i

3RgA

s

k
IL Rks

RSTjk

IC

n

RCTjk
3RgB

(c) Equivalent DC network
Fig. 1: Schematic of the GSU and network transformer modeling used here.
Picture (c) illustrates the effective GIC calculations. For the autotransformer,
Tjk , j is the high voltage (HV) bus and k is the low voltage (LV) bus. Rg
denotes the substation grounding resistance. IH (IL ) and IS (IC ) are the GIC
flows through HV (LV) winding and series (common) winding, respectively.

GIC-Effects During GMDs, the quasi-DC GICs may flow
through transformers with grounded neutral legs. This quasiDC current combines with the normal operating AC current
creating half-cycle saturation and loss of magnetic flux from
the transformer core and leads to several undesirable effects.
The two effects that we consider are eddy current-driven
transformer heating and excess reactive power consumption
from the excess magnetic energy stored in the stray magnetic
flux. Both of these effects are challenging to model from first
principles, and even if such models existed, they would be
too complex to include in the OTS formulation considered
here. Instead, we use a combination of manufacturer test and
specification data and simplified models.
For eddy current-driven transformer heating, we use GIC
capability curves (e.g. see Fig. 2) that may be based on either
3 The

model remains convex for any constant turns ratio by substituting
for |IC |.

αIH +IL
α

4

AC (p.u.)

4.5
ACmax
4.0

− Vi Vj be sin(θj − θi )) ∀eij ∈ E a \ E g
bc
qji = ze (−(be + e )Vj2 + Vi Vj be cos(θj − θi )
2
− Vi Vj ge sin(θj − θi )) ∀eij ∈ E a \ E g
bc
pij + pji = ze re (le + bce qij + ( e )2 Vi2 ) ∀eij ∈ E a \ E g
2
bc
c
qij + qji = ze (xe (le + be qij + ( e )2 Vi2 )
2
bc
− e (Vi2 + Vj2 )) ∀eij ∈ E a \ E g
2
pij + pji = 0, qij + qji = 0 ∀eij ∈ E g

sampled points
constructed quadratic curve

3.5
180oC (2 min.)
3.0
2.5

AC feasible
region

50% ACmax
2.0
1.5
0.20

0.25

0.30

0.35 0.40 0.45
GIC (p.u.)

0.50

0.55

Fig. 2: Fitted curve for thermal GIC capability of a transformer. Here, we used
180◦ as the maximum allowed temperature of transformers for short-term (2
minutes) peak GIC pulses and assumed that a transformer cannot be loaded
to greater than 100% of its MVA limit. The figure shows the coefficients
of the constructed quadratic function (curve), ηe0 , ηe1 and ηe2 , fitted to the
collected 11 (blue) points of the GIC thermal capacity measurements. The
feasible region of the transformer load current is the area under the curve and
is expressed as constraint (3w).

manufacturer acceptance test data or on electromagnetic and
thermal modeling of the transformer design. These curves
provide an upper bound on a feasible operating range in
the space of AC loading and GIC. The upper bound is also
a function of the duration of the combined AC and GIC
loading (typically given for 30 minute and 2 minute durations).
The sampled points (blue) in Fig. 2 are sampled from a
transformer manufacturer’s 2-minute duration curve [26]. Over
a reasonable operating range, these points are well represented
by the best-fit quadratic (red) curve with the feasible operating
region lying below and to the left of the curve.
Excess reactive power losses due to GIC has been studied
in the literature [10], [13], [15], [16]. We adopt the simplified
model in [13] which is shown in Eq.(5). These reactive losses
create voltage sags that can adversely impact system operation.
The previous work has focused on minimizing these losses to
improve system safety. In this manuscript, we explicitly model
the AC power equations (voltage magnitudes) so that we can
enforce voltage limits directly.
B. ACOTS with GIC constraints
A complete ACOTS model with topology reconfiguration
that accounts for GIC-induced transformer thermal heating and
transformer reactive power heating is formulated below.
min

X

c2i (fip )2 + c1i fip + ze (c0i ) +

i∈G,e∈Ei

X

µ(lip + liq )

AC power flow equations
X
X
pij +
pij = fip + lip − dpi − Vi2 gi
eij ∈Ei+

eji ∈Ei−

X

X

qij +

eij ∈Ei+

(3a)

i∈N

∀i ∈ N a

qij = fiq + liq − dqi + Vi2 bi − Qloss
i

(3b)

∀i ∈ N a

eji ∈Ei−

(3c)
− Vi Vj ge cos(θi − θj )
− Vi Vj be sin(θi − θj )) ∀eij ∈ E a \ E g
bc
= ze (−(be + e )Vi2 + Vi Vj be cos(θi − θj )
2
− Vi Vj ge sin(θi − θj )) ∀eij ∈ E a \ E g

pij =

qij

ze (ge Vi2

pji =

ze (ge Vj2

− Vi Vj ge cos(θj − θi )

(3d)

(3e)

p2ij

2
+ qij
= le Vi2
a
2
le = (Iee )
∀e ∈

∀eij ∈ E

a

2
2
p2ij + qij
≤ ze s2e , p2ji + qji
≤ ze s2e
a
0 ≤ Ieea ≤ ze I e ∀e ∈ E a

∀i ∈ N

(3g)
(3h)

(3i)
(3j)
(3k)

a

E
Operational limit constraints

V i ≤ Vi ≤ V i

(3f)

(3l)
∀eij ∈ E a

a

(3m)
(3n)
(3o)

|θi − θj | ≤ ze θ + (1 − ze )θM ∀eij ∈ E a \ E g
ze gpi ≤ fip ≤ ze gpi ∀i ∈ G, e ∈ Ei

(3p)
(3q)

ze gq i ≤ fiq ≤ ze gq i

(3r)

∀i ∈ G, e ∈ Ei

GIC effects on transformers
X
X
d
−
−
z→
z→
e Je −
e Je = −am Vm
−
e∈Em

+
e∈Em

−

X

d
ze ae (Vm
− Vnd ) +

X

d
ze ae (Vnd − Vm
) ∀m ∈ N d

d
d
−
Ied = z→
∀emn ∈ E τ
e ae (Vm − Vn )
d
d
d
d
Iee ≥ Ie , Iee ≥ −Ie
∀e ∈ E τ
a
0 ≤ Ieed ≤ max 2I ê ∀e ∈ E τ
a
−
Ie→
e ≤

Qloss
i

∀ê∈E a
0
ηe + ηe1 Ieed

=

(3s)

−
enm ∈Em

+
emn ∈Em

X

+ ηe2 (Ieed )2 ∀e ∈ E τ
ke Vi Ieed ∀i ∈ N a

(3t)
(3u)
(3v)
(3w)
(3x)

e∈Eiτ

ze ∈ {0, 1}

∀e ∈ E a

(3y)

The objective function (3a) minimizes total generator dispatch costs and load shedding costs. Constraints (3b) – (3r)
describe system constraints for the buses and branches in the
AC circuit. Constraints (3b) and (3c) represent the nodal real
and reactive power balance, including the increased reactive
power losses (demand) due to GICs. Constraints (3d) through
(3g) model the AC power flow on each transmission line with
on-off variables ze . The flow on any line is forced to zero if the
line is switched off. Constraints (3h) through (3j) model power
loss equations associated with AC power flow. In constraint
(3j), fictitious lines between output terminals of generators and
their injection buses are modeled as transportation edges (i.e.,
|pij | = | − pji | = fip , |qij | = | − qji | = fiq ∀eij ∈ E g ).
Nonconvex constraint (3k) evaluates current magnitude lij , an
auxiliary variable introduced to bound the squared AC current
flow magnitude in constraint (3l). Constraints (3m) through
(3r) describe the operational limits of the grid; constraint
(3m) models operational thermal limits of lines in both directions. Constraint (3o) limits the voltage magnitude at buses.
Constraint (3p) applies appropriate bounds on phase angle
difference between two buses when the line exists. Constraints
(3q) and (3r) model the availability and capacity of power
generation. A generator is offline if its line is switched off.

5

The DC circuit and the effects associated with the GMD
are formulated in constraints (3s)-(3w). Recall that we link
an edge, e ∈ E d in the DC circuit to an edge in the AC
−
circuit with →
e . Also recall that the HV primary winding of
GSU transformer eij ∈ E a is modeled by introducing a node
and edge in the DC circuit (node m in Fig. 1(c)). Similarly,
the common winding of autotransformer eij ∈ E a is modeled
in the DC circuit by introducing additional nodes and edges
(see Fig. 1(c)). By using these notations, constraints (3s) and
(3t) calculate the GIC flow on each DC line by applying
Kirchhoff’s current law. The GIC on a line is determined by
the induced current source and the quasi-dc voltage difference
between two buses [12]. GIC flow is forced to 0 by ze when e
is switched off. Since the value of Ied can be negative, decision
variables Ieed are introduced to model the magnitude (absolute
value) of GIC flows (i.e., Ieed = |Ied |). Instead of introducing
additional discrete variables, constraint (3u) is used to model
and relax the magnitude of Ied . Constraint (3v) denotes the
maximum allowed value of GIC flowing through transformers.
We assume this limit is twice the upper bound of AC flows
in the network. Constraint (3w) guarantees that the hot spot
temperature of transformers due to the combination of AC
and GICs is below the thermal limits for peak GIC. Constraint
(3x) computes the reactive power load due to GIC transformer
saturation [1], [10], [16], [25] by using the effective GIC on
the primary winding in GSU transformers and the common
winding in autotransformers (Eiτ ). The couplings between AC
power flows and GIC occur in constraints (3c), (3w), and (3x).

on/off variable ze ∈ {0, 1}, a disjunctive quadratic relaxation
is used to convexify the nonlinear function ze cos(θij ) by inR
troducing a new variable cs
b ij ∈ hze cos(θij )i , as formulated
in (6).
1 − cos(θ) 2
(θij + (ze − 1)(θu )2 )
(θ)2
b ij ≤ ze
ze cos(θ) ≤ cs
cs
b ij ≤ ze −

The ACOTS with GIC constraints is a mixed-integer, nonconvex optimization problem that is generally computationally
very difficult to solve. We adopt the convex relaxations developed by [21] and later show that the results obtained with
the relaxation is (empirically) tight. We now discuss the key
features of the relaxations extended to the problem with GIC.
Handling bilinear terms Given any two variables xi , xj ∈
R, the McCormick relaxation is used to linearize the bilinear
MC
product xi xj by introducing a new variable xc
.
ij ∈ hxi , xj i
The feasible region of xc
ij is defined by inequalities (4). Note
that the MC relaxation is exact if one variable is binary.

sbij ≤ cos(θ̄/2)θij + ze (sin(θ̄/2) − θ̄/2 cos(θ̄/2))
+ (1 − ze )(cos(θ̄/2)θ

M

− (1 − ze )(cos(θ̄/2)θ

M

(7c)

Based on the above relaxations, we replace the non-convex
constraints in (3d), (3e), (3f) and (3g) with equations (8):
pij = ge zc
v ij − ge w
ccij − be ws
c ij ∀eij ∈ E a \ E g
c
b
qij = −(be + e )c
z v ij + be w
ccij − gij ws
c ij ∀eij ∈ E a \ E g
2

cs
b ij ∈ hze cos(θij )iR , sbij ∈ hze sin(θij )iR
M C−q

, zc
v ij ∈ hzij , vbi i

MC

,

(8b)

(9a)
(9b)

w
bij ∈ hVi , Vj iM C ,

(9c)

MC

w
ccij ∈ hw
bij , cs
b ij i

, ws
c ij ∈ hw
bij , sbij i

MC

(9d)

Other nonconvex constraints Further, non-convex constraints
(3k) and (3l) are relaxed to a convex, rotated second-order
conic constraint by using the introduced lifted variable vbi (for
(3k)) as follows:
2
p2ij + qij
≤ le vbi ∀eij ∈ E a ,
a
le ≥ (Ieea )2 , le ≤ (I e )Ieea ∀e ∈ E a

(10a)
(10b)

Convex relaxations of the reformulated thermal heating limit
constraint (3w) and excess reactive power losses equation (3x)
are stated here in equations (11):

(4a)

a
0
1 ed
2 bd
−
Ie→
∀e ∈ E τ
e ≤ ηe + ηe Ie + ηe Ie
X
d
Qloss
=
ke VcI e ∀i ∈ N a
i

xi ≤ xi ≤ xi , xj ≤ xj ≤ xj

(4e)

M C−q
Ibed ∈ hIeed i

On/off trigonometric terms In constraints (3d), (3e), (3f) and
(3g), if the line eij is switched off, [21] suggests the following
procedure to deactivate the associated trigonometric terms:
Given the phase angle difference variable θij = θi − θj and

(8a)

where, the new variables zv
d
ccij and ws
c ij , admit feasible
ij , w
regions as given in equations (9):

(4b)
(4c)
(4d)

(5a)
(5b)
(5c)

(7b)

+ 1)

ze sin(−θ) ≤ sbij ≤ ze sin(θ)

xc
ij ≥ xi xj + xj xi − xi xj
xc
ij ≤ xi xj + xj xi − xi xj
xc
ij ≤ xi xj + xj xi − xi xj

xbi ≥ x2i
xbi ≤ (xi + xi )xi − xi xi
xi ≤ xi ≤ xi

(7a)

+ 1)

sbij ≥ cos(θ̄/2)θij − ze (sin(θ̄/2) − θ̄/2 cos(θ̄/2))

xc
ij ≥ xi xj + xj xi − xi xj

Quadratic terms Given a variable xi ∈ R, a second-order
conic relaxation can be applied to convexify the quadratic term
M C−q
x2i by introducing a new variable xbi ∈ hxi i
, as defined
in equation (5).

(6b)

Similarly, for ze sin(θij ), a disjunctive polyhedral relaxation
R
is applied by introducing a new variable sbij ∈ hze sin(θij )i ,
as described in equation (7).

vbi ∈ hVi i

C. Convex Relaxations

(6a)

(11a)
(11b)

e∈Eiτ

∀e ∈ E τ

d
MC
VcI e ∈ hVi , Ieed i

(11c)
a

∀i ∈ E , ∀e ∈

Eiτ

(11d)

Convex relaxations of the GIC injection constraints (3s) and
(3t) are described in equation (12).
X

−
z→
e Je −

X

ae zc
v de +

+
e∈Em

Ied = ae zc
v de
zc
v de

∈

d
−
z→
e Je = −am Vm

−
e∈Em

+
e∈Em

−

X

X

ae zc
v de ∀m ∈ N d

(12a)

−
e∈Em

∀e ∈ E τ

d
−
hz→
e , (Vm

−

MC
Vnd )i

(12b)
∀emn ∈ E

d

(12c)

6

III. C ASE S TUDY
In this section, we analyze the performance and sensitivity
of a power system when exposed to varying strengths of geoelectric fields induced by GMDs. We use a modified version of
single area IEEE RTS-96 system [27]. Its size is comparable to
previous work [16] that considered minimization of the quasistatic GICs and not a full AC-OPF with topology control.
The derived and modified parameters of IEEE RTS-96 are
presented in Table I–II. We arbitrarily placed the system
in western Pennsylvania to give the model a geographic
orientation. We assume the cost of shedding load is twice
the cost of the most expensive generator. We performed all
computations using the high performance computing resources
at Los Alamos National Laboratory with Intel(R) Xeon(R)
CPU E5-2660 v3 @2.60GHz processor and 120 GB memory
installed. All cases were solved using CPLEX 12.7.0 (default
options). Knitro 10.2.1 (default options) was used as the local
solver. JuMP was used as an algebraic modeling language [28].
For reference, the peak geo-electric field during the HydroQuebec event of 1989 was 2 V/km (≈3.2 V/mile) [29], [30].
References [31] and [32] suggests that 100-year GMDs could
cause 5 V/km (≈9 V/mile) and 13 V/km (≈21 V/mile), respectively, at some high-latitude locations. In our case studies, we
consider middle ground, but still extreme, geo-electric fields
of 12 V/mile and 14 V/mile. We also study the directionality
of the event by considering field directions between 0◦ and
360◦ spaced by 10◦ .
To analyze the benefits for GIC mitigation of generator
dispatch and load shedding and the combined effects of those
two controls plus topology control, we studied three cases. To
describe these cases, we define z∗x and c∗x to be the optimal
topology (line on/off) decisions and objective (minimum total
costs), respectively, for case x. The models are defined below
as (Mx ). The solutions of z∗x and c∗x are obtained from the
convex relaxations of Mx described in Section II-C. The cases
we consider are:
1) C1: The ACOTS model neglecting GIC effects (c∗o , z∗o ):
Mo := Min{(3a): (3b)-(3r); (3y); Qloss = 0}
2) C2: The ACOTS with GIC effects ( c∗gmd , z∗gmd ):
Mgmd := Min{(3a): (3b)-(3y)}
3) C3: The ACOPF (fixed z = z∗o ) with GIC effects (c∗f ,
z∗o ):
Mf := Min{(3a): (3b)-(3y); z = z∗o }
4) C4: The ACOPF (fixed z = 1) with GIC effects (c∗1 , 1):
M1 := Min{(3a): (3b)-(3y); z = 1}
Case C1 defines the topology z∗o and evaluates the objective
that results from neglecting the effects of GICs. Case C3
evaluates the new cost c∗f that results from mitigating GIC
with generation dispatch on the topology of C1. Case C4 is
similar to case C3, but all lines are closed. Case C2 considers
the effects of both generation dispatch and topology control.
All results in this section are based on the convex relaxation
(except for Section III-G which evaluates the quality of the
relaxation by recovering feasible solutions to the original nonconvex formulation).
c∗o

TABLE I: Transformer data. In this manuscript, all network transformers are
auto-transformers. The transformer winding resistance and k are estimated
based on the test cases provided in [11], [12].
Name

Type

Resistance
W1
(Ohm)

Bus
No.

Resistance
W2
(Ohm)

Bus
No.

Line
No.

k
(p.u.)

A
A
A
A
A
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G
G

Auto
Auto
Auto
Auto
Auto
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU
GSU

0.12
0.12
0.12
0.12
0.12
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3
0.3

3
9
9
10
10
1
1
1
1
2
2
2
2
7
7
7
13
13
13
14
15
15
15
15
15
15
16
18
21
22
22
22
22
22
22
23
23
23

0.18
0.18
0.18
0.18
0.18
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A

24
11
12
11
12
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

7
14
15
16
17
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8

1
2
3
4
5
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

A. GIC Modeling Validation
To validate the GIC modeling in Eq.(3s)–(3v), we tested
our model on the 6-bus system given in Appendix II of
[12] and compared our solution with the results provided in
this reference. Table III displays the GIC flows obtained by
solving Eq.(3s)–(3v) with fixed ze = 1 (in [12] the GIC flows
are calculated without line switching options). This solution
matches the results found in the reference (see Eq.(B.5)–
Eq.(B.10)).
B. Case C1: Potential Damage by GICs
Under normal circumstances without GMDs, line switching
decisions are determined by economic dispatch. More specifically, the optimal system topology is obtained by solving an
ACOTS model without the GIC-effects constraints (Case C1).
Figure 3(a) shows the optimal normal topology, z∗o , where
some generators are not injecting real or reactive power. For
example, generators 16 through 20 are shut down at node 15,
and their GSU transformers are disconnected from the network
using the circuit breakers. Referring to Fig. 1, we note that
this action does not significantly affect the topology of the
AC network, which is only affected by switching transmission
lines. This action removes GSU transformer ground points
from the DC network topology over which the GICs flow.

7

TABLE II: Power system model parameters. The nominal line length parameters of a single area of RTS-96 [27] are used to perform an approximate
geospatial layout of the power system nodes. (b) The substation grounding
resistance GR is estimated from typical values of grounding resistance of
substations provided in [33]. (c) The original line parameters reo and xoe are
scaled by the ratio kle of the new to original line lengths.
(a) Transmission line data
Line
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

From
Bus

To
Bus

Length
(miles)

1
1
1
2
2
3
3
4
5
6
7
8
8
9
9
10
10
11
11
12
12
13
14
15
15
15
15
16
16
17
17
18
18
19
19
20
20
21

2
3
5
4
6
9
24
9
10
10
8
9
10
11
12
11
12
13
14
13
23
23
16
16
21
21
24
17
19
18
22
21
21
20
20
23
23
22

3.98
53.15
22.78
33.16
44.49
33.56
0.00
26.89
23.38
19.96
16.04
43.51
43.51
0.00
0.00
0.00
0.00
35.95
33.98
35.95
70.48
57.39
27.36
12.18
35.44
35.44
38.43
18.77
18.57
10.75
72.84
17.96
17.96
29.97
29.97
15.59
15.59
51.83

(b) Substation data
Name

Latitude

Longitude

GR
(Ohm)

SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB
SUB

40.44
40.44
40.90
40.70
40.70
41.08
40.50
40.53
41.03
41.22
41.48
41.45
41.63
41.86
42.01
41.77
42.01
41.95
42.41
42.02

-78.80
-78.73
-79.61
-79.26
-79.07
-78.61
-78.20
-78.50
-78.99
-78.35
-79.26
-79.71
-79.75
-79.94
-79.86
-79.45
-78.95
-79.52
-78.73
-78.65

0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

all orientations of the GMD. For example, when the event
is oriented at 10◦ , GSU transformers 21, 22 and 23 are
overheated. As the event is shifted to 80◦ , then network
transformer 1 is the transformer at risk. These results provide a
baseline to evaluate alternative operating paradigms that ensure
system security.

(a) Optimal topology z∗o in case C1
(c) Other parameters
µ
a
Ie

$ 1000 /MW (or MVar)

re
xe

Te / min{V i , V j }
(βe )reo
(βe )xo
e

θ

30◦

12 V/mile
14 V/mile

170
160

22, 23

150
140
130
120

23

110
100
90

1

φo 80
70

21

60
50
40

21, 22

21

21, 22

30

TABLE III: The test case in [12], Appendix II
Variable Name

Transformer

GIC flow (amps)

I12
Is
I34
IT 1
Ic
IT 3

T2 series
T1
T2 common
T3

-627.02
-763.26
-763.26
627.02
136.24
-763.26

Case C1 assumes that generation and system topology are
optimized for cost while neglecting the impact of GICs. This
impact is calculated using Eqs. (3s) through (3u) to evaluate
the feasibility of thermal limit constraint (3v). Figure 3(b)
shows how many GSU and network transformers would be
overheated under C1 depending on the direction and strength
of the GMD. Figure 3(b), presents results for c∗gmd from 0 to
180◦ because the strength of the geo-electric field is uniform,
and the effects do not depend on field direction. For example,
when the electric field is 12 V/mile, the GSU transformer
23 (at node 18) is overheated when the event is oriented
between 100◦ and 170◦ . When the strength is increased to
14 V/mile, one or more transformers are overheated at almost

20

21, 22, 23
21, 22, 23

22

10

22, 23

0

0

1

2

3

4

Number of overheated transformers

(b) Potential damage of transformers by GICs
Fig. 3: Evaluation of the power system in Table I–II for case C1. (a) The grey
nodes are loads. The blue lines indicate network transformers. The green
and magenta nodes indicate GSU transformers. Transformer IDs are listed
next to the node. The solution to case C1 does not allow reconfiguration
of the topology via the network transformers, however, if a generator is not
injecting real or reactive power into the network, its GSU is disconnected
using the circuit breakers in Fig. 1. Generation nodes with disconnected GSU
transformers are magenta and their (network or GSU) transformer IDs are
marked as magenta as well. (b) The case C1 solution is tested by applying
geo-electric fields of strength 12 V/mile and 14 V/mile for all directions. The
label above each bar indicates IDs of overheated transformers. Red and green
are used to label network transformers and GSU transformers, respectively.

C. Case C2: GIC Mitigation via ACOTS
Using case C2, the cost benefits of simultaneous controlling
generation dispatch and network topology to mitigate GIC
effects are evaluated.
1) Cost Analysis: For geo-electric field strengths of 12
V/mile and 14 V/mile, case C2 is solved for orientations of

8

the field from 0 to 360◦ , which results in a total cost c∗gmd
∗
(see Fig. 4) and topology ygmd
(discussed later). Figure 4
∗
only presents results for cgmd from 0 to 180◦ because of the
symmetry discussed above.
The results in Fig. 4 show that the directions of the
geo-electric field are not all equivalent because the cost of
mitigation c∗gmd varies with direction. The most costly GMDs
occur when the event is oriented between 20◦ and 140◦ . The
increase in cost between 12 V/mile and 14 V/mile is primarily
due to changes in generator dispatch and is not significant.
For example, the difference in cost between the 12V and 14V
per mile case is 1.40% when the GMD is oriented at 60◦ .
Moreover, the dispatch cost is smaller when GIC effects are
neglected (Case C1). However, the transformer thermal limit
constraints are violated when GIC effects are applied to the
network (as seen in Fig. 3(b)). Thus, there is an implicit higher
cost associated with replacing the damaged equipment and
unexpected load shed when the transformer fails.
12 V/mi
14 V/mi
No GIC

Total Cost ($)

56298

the details of the orientation and the difficulty in making
accurate predictions of geo-electric field direction suggest that
the ACOTS formulation should be extended to a stochastic
formulation over the field direction in future work.
The results displayed in Fig. 6 for different geo-electric field
orientations suggests similar conclusions. At a fixed 14 V/mile
in Fig. 6, the optimal topology solutions switch off several
long transmission lines oriented along the geo-electric field,
but some transmission lines still display significant sensitivity
to orientation.
Comparing Fig. 6 (14 V/mile) with Fig. 5 (12 V/mile)
shows that some topology solutions at low field strength
persist to higher field strength, however significantly more
transmission lines are switched off to avoid large GIC in
the network and in GSU and network transformers. The
properties of the topology solutions for different geo-electric
field strengths again suggests that the ACOTS solution should
be extended to a stochastic or robust formulation over field
strength. Finally, we note that while the solution adjusts the
topology, it does not create islands—a mitigation strategy that
is sometimes suggested. However, islands could form in larger,
more complex networks.

55998
55698
55398
55098
0

20

50

80
φo

110

140

170

Fig. 4: The total cost c∗gmd for case C2 for different geo-electric field
orientations and strengths.

2) Topology Control Analysis: In Fig. 4, the topology of the
network varies with the strength and direction of the GMD.
In the remainder of this section, we focus on GMD events
oriented between 80◦ and 110◦ because these directions are
most sensitive to GMD. Figures 5 and 6 display the network
topology for geo-electric field strengths of 12 V/mile and 14
V/mile, respectively.4
In Fig. 5, in the 80◦ geo-electric field case, only one
transmission line is switched off (1,5). In the 90◦ geo-electric
field case, three transmission lines are switched off. Two
are intuitively long lines oriented along the geo-electric field
and incur larger GICs. Switching these lines off removes a
significant vulnerability. The third line is nearly perpendicular
to the geo-electric field and is also switched off. This counterintuitive topology control is being used to reroute power flow
away from other, more susceptible transmission lines.
In the 100◦ orientation case, some of the lines in the
90◦ case remain in the solution and some disappear. All
transmission lines are switched on when the event orientation
is at 110◦ . Thus, the sensitivity of the topology solution to
4 We note that there are multiple generators located at buses 1, 2 and 15.
Generators 1 and 2 (at bus 1) have the same cost and capacity, as do generators
5 and 6 at bus 2 and generators 16 through 20 at bus 15. Thus, there are
equivalent dispatch solutions.

(a) 12 V/mile, 80◦

(b) 12 V/mile, 90◦

(c) 12 V/mile, 100◦

(d) 12 V/mile, 110◦

Fig. 5: Topology solutions for case C2 at 12 V/mile strength and orientations
from 80◦ –110◦ . Switched off lines are colored magenta and the IDs of unused
generators are labeled beside their connected substations.

D. Case C2 versus Case C3: Cost Benefits of Topology
Optimization
The inclusion of topology control into the ACOTS formulation increases the complexity of the problem, but it also
provides significant cost savings over a less complex ACOPF.
The cost savings is evaluated by comparing case C2 (where
topology control is allowed) with case C3 (where the topology
is fixed to that found in case C1). Figure 7 displays the
percentage cost savings of C2 (ACOTS) over C3 (ACOPF) for
field strengths between 12 and 14 V/mile and field directions

9

TABLE IV: Percentage of the total cost in cases C2 and C3 due to generator
dispatch and load shedding. For 12 V/mile and 14 V/mile strengths, the
average, minimum and maximum percentage of total cost is computed over
the geo-electric field orientation from 0◦ –180◦ .
Dispatch Cost(%)
Strength

(a) 14 V/mile, 80◦

Avg.

Min.

Max.

Avg.

Min.

Max.

12 V/mile

C2
C3

100.0
97.8

100.0
90.8

100.0
100.0

0.0
2.1

0.0
0.0

0.0
9.2

14 V/mile

C2
C3

100.0
86.1

100.0
66.46

100.0
100.0

0.0
13.9

0.0
0.0

0.0
33.5

(b) 14 V/mile, 90◦

(a) The total cost c∗f in Case C3.
(c) 14 V/mile, 100◦

(d) 14 V/mile, 110◦

Fig. 6: Same as Fig 5 but for a geo-electric field strength of 14 V/mile.

between 0◦ –180◦ . Under the most severe GMD conditions
explored, the benefit of topology control is as much as 54%.
Table IV further breaks down the cost savings of case C2
over case C3 into generator dispatch costs and load shedding
costs. For the 14 V/mile field strength case, the topology
control in case C2 enables nearly all of the load to be served.
In contrast, the fixed topology in C3 results in load shedding
costs of 13.9% on average and 33.54% in the worst case.

45

48%

36
40%

27

32%

18

24%

9

Cost saving (%)

54

56%

8%

80

φo

110

140

170

12.0

Elec

50

tric

12.6
20

field

13.2

0%

(V/
m

13.8

ile)

0

16%

0

Load shedding Cost(%)

Case

Fig. 7: Combined savings from generator dispatch and load shedding costs
enabled by the optimal topology z∗gmd found by ACOTS relative to the
dispatch and load shedding cost incurred by the ACOPF of case C3 with
the topology fixed to z∗0 .

E. Case C3 versus Case C4: Performance of Network Reconfiguration
The results displayed in Figure 8 compare the topology of
case C3 (where topology is fixed to case C1) with case C4
(where all lines and generators). Similar to the results in Fig.

(b) The total cost c∗1 in Case C4.

Fig. 8: Cost comparison of Case C3 and C4 for different geo-electric field
orientations and strengths.

8(a), the total cost, c∗1 , varies with event direction. Ignoring
GIC effects also induces a lower operating cost (Fig. 8(b)).
However, like the results of Fig. 3(b), we also observed that the
transformer thermal limit constraints (3w) are violated when
GIC effects are neglected. For example, if the field strength
is 14 V/mile and has an orientation of 20◦ , GSU transformers
22 and 23 overheat. Comparing Fig. 8(a) (C3) with Fig. 8(b)
(C4) shows that the optimal topology control found in case
C1 induces a higher cost than case C4 for orientations through
40◦ to 130◦ under the field of 14 V/mile. This is due to forced
disconnect of generators in Case C1 (e.g., generator 14 at bus
13) which could be dispatched more effectively, when no line
can be switched off, to mitigate GIC effects.
F. Computational Analysis
1) Computational Speed: Table V summarizes the computational time properties of ACOTS with GIC subject to convex
relaxations (from section II-C) by solving them to optimality.
We observed that the times are higher at larger geo-electric
field strength, likely because of the increased complexity due
to larger number of possible topology changes. Though the
computations are time-intensive (in Table V), we observed that
by terminating the solver at larger optimality gaps (say 5%),
a solution was obtained within 320 seconds, which was 2-3
orders of magnitude quicker than solving to optimality.
2) Scalability to Larger Network: The computational time
required to solve small test cases implies that a key limitation
for practical deployment of this model is scalability. This is
not surprising as solution methods for OTS suffer similar
limitations. However, in the case of GMD mitigation, highquality solutions that are close to optimality are often sufficient. On UIUC 150-bus system, figure 9 shows the feasible
solution costs when the optimization is terminated with a 5%
optimality gap. Like are other results, the cost of dispatch

10

TABLE V: Computational time for geo-electric field strengths of 12 V/mile
and 14 V/mile on the RTS system. The average, minimum, maximum and
standard deviation of solving time are presented over the geo-electric field
orientation from 0◦ –360◦ . Results without parentheses present elapsed time
for solving to optimality. Solutions displayed in parentheses denote solving
time when the optimization is terminated with a 5% optimality gap.

Of∗ eas := min O(fip , lip , liq )
s.t.

O(fip , lip , liq )

≤

∗
Olb
(1

(13a)
+ δ),

Constraints (3b) − (3x),
Wall Time (sec)
Strength(V/mile)

Avg.
167.1
(10.6)
1249.9
(56.0)

12
14

Min.
31.5
(10.2)
32.6
(10.3)

Max.
1066.9
(11.0)
6676.3
(314.1)

264.1
(0.3)
1880.9
(85.8)

Total Cost ($)

varies with direction and is higher than the dispatch cost when
GIC is ignored. In Table VI, we provide the computational
times required to obtain results displayed in Fig. 9. These
results suggest that it is practical to use heuristic methods
in conjunction with state-of-the-art convex relaxations to find
high-quality solutions to larger scale systems on the time
scales required for GMD mitigation efforts.
12 V/mile
No GIC

353352
351852
350352
348852
347352
345852
344352
342852
341352
339852

ze =

Std. dev.

ze∗

(13b)
(13c)

a

∀e ∈ E .

(13d)

O(fip , lip , liq )

where
represents the objective function (3a).
Constraint (13b) specifies that the objective function cost must
∗
be within a small percentage δ of the lower bound Olb
, where
δ is a specified parameter (3% for testing). Thus, formulation
(13) guarantees a feasible solution (close to global optimum)
if there exists one for the specified δ. A similar approach has
been shown to be effective for OPF-based problems in [34].
The properties of the optimality gap of the test case (single
area IEEE RTS-96 system) from these studies, summarized in
Table VII, suggest that the relaxed solution is always within
3% of the optimal solution, indicating that the relaxation is
empirically tight to the original MINLP. It is also noteworthy
to mention that the convergence time and quality of the local
solver (Knitro 10.2.1) solutions were tremendously improved
by solving the formulation (13). For instance, we observed
gaps up to 70% (instead of 3%) by solving formulation (13)
without the objective-cost constraint in (13b).
TABLE VII: Optimality gaps for the RTS-96 system between the lower bound
∗ ) and the feasible solution (O ∗
(Olb
f eas ) recovered for the original nonconvex
model. Values shown are evaluated over various geo-electric field orientations.

0

20

50

80

φo

110

140

Optimality Gap (%)

170

Fig. 9: Feasible solution costs of UIUC 150-bus system for geo-electric field
strength of 12 V/mile in Case 2.
TABLE VI: Computational time of the UIUC 150-bus system. The average,
minimum, maximum and standard deviation of solving time are presented
over the geo-electric field orientation from 0◦ –360◦ .

Strength (V/mile)

Avg.

Min.

Max.

Std. Dev.

12
14

0.9
1.8

0.01
0.01

2.99
3.00

1.4
1.5

Given the optimality gaps observed here, it is worth noting
that the smaller cost fluctuations (up to 3% in Figure 4) could
be due to the relaxations.

Wall Time (sec)
Strength(V/mile)
12
14

Avg.

Min.

Max.

Std. dev.

996.3
1162.3

648.9
508.1

1586.6
2896.1

290.4
750.2

G. Recovering AC Feasible Solutions
The ACOTS with GMD constraints formulated in this
manuscript is solved applying several hierarchical convex
relaxations of the AC power flow physics. Since the solutions
obtained may not necessarily lie in the original nonconvex
feasible region, we present a simple methodology to test
the quality of the relaxed solutions and obtain AC feasible
solutions.
In [21], for OPF-based problems, it is empirically shown
that the lower bounds obtained from convex quadratic relaxations are mostly close to globally optimal objective values.
Thus, we exploit this fact and apply an objective-cost-based
constraint to the problem of ACOTS with GMD. Let the
optimal value of formulation (3) with convex relaxations be
∗
Olb
and let ze∗ ∀e ∈ E a be the respective optimal topology.
For the fixed topology ze∗ , we solve the following original
nonlinear, nonconvex program (without integer variables):

IV. C ONCLUSIONS
We formulated a detailed topology control optimization
model to mitigate the impacts of GMD on electrical transmission systems. The mathematical formulation minimizes the total generation dispatch and load shedding subject to nonconvex
AC power flow physics, effects of geomagnetically-induced
currents on transformer heating and transformer reactive power
consumption. Further, we leveraged recently developed convex
relaxation approaches to handle the nonlinearites due to AC
transmission switching and GIC constraints, which we subsequently observed to provide near global optimum solutions.
While this paper has made contributions in showing that
switching can mitigate the impacts of GMD events, there
remain a number of important future directions. For example,
new algorithms are needed to solve larger problems. Here, the
ACOTS with GMD is naturally posed as a 2-stage program
with topology decisions in the master problem. Thus decomposition algorithms, like Benders’, are a natural direction to
consider. Second, based on our empirical observations, convex
relaxation solutions are often tight, thus local search techniques, like meta-heuristics and state-of-the-art global search

11

methods, could yield high quality solutions quickly [35], [36].
In addition, there are a number of modeling enhancements
that need to be considered. For example, integration of N-1
security (contingency) constraints are important to increase the
resiliency of transmission systems under GMD extreme events.
Moreover, capturing other effects of GMD on transformers,
the modeling of time-extended variations in geo-electric field
strengths will be important. Further, there is often uncertainty
in predictions of direction and strength of the GMD event,
thus it will be important to development methods that produce
solutions that are robust to errors in predictions. Finally, this
paper modeled GSU and auto transformers, but there are other
types of transformers like GWye-GWye Auto and Delta-wye
that will need to be modeled.
ACKNOWLEDGMENT
The work was funded by the Center for Nonlinear Studies at
LANL, the Defense Threat Reduction Agency project Advancing Knowledge of Networks for Understanding Robustness,
and the LANL laboratory directed research and development
project Impacts of Extreme Space Weather Events on Power
Grid Infrastructure: Physics-Based Modeling of Geomagnetically Induced Currents (GICs) During Carrington-Class Geomagnetic Storms. It was carried out under the auspices of
the NNSA of the U.S. DOE at LANL under Contract No.
DE-AC52-06NA25396.
R EFERENCES
[1] V. D. Albertson, J. M. Thorson Jr, R. E. Clayton, and S. C. Tripathy,
“Solar-induced-currents in power systems: cause and effects,” Power
Apparatus and Systems, IEEE Trans. on, no. 2, pp. 471–477, 1973.
[2] A. VD, J. THORSON, and S. MISKE, “Effects of geomagnetic storms
on electrical power-systems,” IEEE TRANS. ON POWER APPARATUS
AND SYSTEMS, no. 4, pp. 1031–1044, 1974.
[3] V. Albertson, B. Bozoki, W. Feero, J. Kappenman, E. Larsen, D. Nordell,
J. Ponder, F. Prabhakara, K. Thompson, and R. Walling, “Geomagnetic
disturbance effects on power systems,” IEEE Trans. on power delivery,
vol. 8, no. 3, pp. 1206–1216, 1993.
[4] “Effects of geomagnetic disturbances on the bulk power system - technical report,” North American Electric Reliability Corporation, 2013.
[5] “Executive order–coordinating efforts to prepare the nation for space
weather events,” White House Office of the Press Secretory, 2016.
[6] Federal Energy Regulatory Commission, “156 FERC 61,215: Reliability
Standard for Transmission System Planned Performance for Geomagnetic Disturbance Events,” Tech. Rep., 2016.
[7] I. A. Erinmez, J. G. Kappenman, and W. A. Radasky, “Management of
the geomagnetically induced current risks on the national grid company’s
electric power transmission system,” Journal of Atmospheric and SolarTerrestrial Physics, vol. 64, no. 5, pp. 743–756, 2002.
[8] P. Cannon, M. Angling, L. Barclay, C. Curry, C. Dyer, R. Edwards,
G. Greene, M. Hapgood, R. B. Horne, D. Jackson et al., Extreme space
weather: impacts on engineered systems and infrastructure. Royal
Academy of Engineering, 2013.
[9] Q. Qiu, J. A. Fleeman, and D. R. Ball, “Geomagnetic disturbance:
A comprehensive approach by american electric power to address the
impacts.” Electrification Magazine, IEEE, vol. 3, no. 4, pp. 22–33, 2015.
[10] T. J. Overbye, K. S. Shetye, T. R. Hutchins, Q. Qiu, and J. D. Weber,
“Power grid sensitivity analysis of geomagnetically induced currents,”
Power Systems, IEEE Trans. on, vol. 28, no. 4, pp. 4821–4828, 2013.
[11] R. Horton, D. Boteler, T. J. Overbye, R. Pirjola, and R. C. Dugan, “A
test case for the calculation of geomagnetically induced currents,” Power
Delivery, IEEE Trans. on, vol. 27, no. 4, pp. 2368–2373, 2012.
[12] “Computing geomagnetically-induced current in the bulk-power system
- application guide,” North American Electric Reliability Corporation,
2013.

[13] T. J. Overbye, T. R. Hutchins, K. Shetye, J. Weber, and S. Dahman,
“Integration of geomagnetic disturbance modeling into the power flow:
A methodology for large-scale system studies,” in North American
Power Symposium (NAPS), 2012. IEEE, 2012, pp. 1–7.
[14] L. Bolduc, M. Granger, G. Pare, J. Saintonge, and L. Brophy, “Development of a DC current-blocking device for transformer neutrals,” Power
Delivery, IEEE Trans on, vol. 20, no. 1, pp. 163–168, 2005.
[15] Y. Liang, H. Zhu, and D. Chen, “Optimal blocker placement for
mitigating the effects of geomagnetic induced currents using branch
and cut algorithm,” in North American Power Symposium (NAPS), 2015.
IEEE, 2015, pp. 1–6.
[16] H. Zhu and T. J. Overbye, “Blocking device placement for mitigating
the effects of geomagnetically induced currents,” Power Systems, IEEE
Trans. on, vol. 30, no. 4, pp. 2081–2089, 2015.
[17] B. Kovan and F. de León, “Mitigation of geomagnetically induced
currents by neutral switching,” Power Delivery, IEEE Trans. on, vol. 30,
no. 4, pp. 1999–2006, 2015.
[18] R. Girgis, K. Vedante, and G. Burden, “A process for evaluating the
degree of susceptibility of a fleet of power transformers to effects of
gic,” in T&D Conference and Exposition, 2014 IEEE PES. IEEE,
2014, pp. 1–5.
[19] L. Marti, A. Rezaei-Zare, and A. Narang, “Simulation of transformer
hotspot heating due to geomagnetically induced currents,” Power Delivery, IEEE Trans. on, vol. 28, no. 1, pp. 320–327, 2013.
[20] B. Kocuk, S. S. Dey, and X. Sun, “New formulation and strong MISOCP
relaxations for AC optimal transmission switching problem,” IEEE
Trans. on Power Systems, 2017.
[21] H. Hijazi, C. Coffrin, and P. Van Hentenryck, “Convex quadratic
relaxations for mixed-integer nonlinear programs in power systems,”
Mathematical Programming Computation, pp. 1–47, 2014.
[22] X. Bai, H. Wei, K. Fujisawa, and Y. Wang, “Semidefinite programming
for optimal power flow problems,” International Journal of Electrical
Power & Energy Systems, vol. 30, no. 6, pp. 383–392, 2008.
[23] H. Nagarajan, R. Bent, P. Van Hentenryck, S. Backhaus, and E. Yamangil, “Resilient transmission grid design: AC relaxation vs. DC
approximation,” arXiv preprint arXiv:1703.05893, 2017.
[24] H. Nagarajan, E. Yamangil, R. Bent, P. Van Hentenryck, and S. Backhaus, “Optimal resilient transmission grid design,” in 2016 Power
Systems Computation Conference (PSCC). IEEE, 2016, pp. 1–7.
[25] K. Zheng, D. Boteler, R. J. Pirjola, L.-g. Liu, R. Becker, L. Marti,
S. Boutilier, and S. Guillon, “Effects of system characteristics on
geomagnetically induced currents,” Power Delivery, IEEE Trans. on,
vol. 29, no. 2, pp. 890–898, 2014.
[26] “C57.163-2015 - IEEE guide for establishing power transformer capability while under geomagnetic disturbances.”
[27] P. Wong, P. Albrecht, R. Allan, R. Billinton, Q. Chen, C. Fong,
S. Haddad, W. Li, R. Mukerji, D. Patton et al., “The ieee reliability
test system-1996. a report prepared by the reliability test system task
force of the application of probability methods subcommittee,” Power
Systems, IEEE Trans. on, vol. 14, no. 3, pp. 1010–1020, 1999.
[28] I. Dunning, J. Huchette, and M. Lubin, “JuMP: A modeling language
for mathematical optimization,” arXiv preprint arXiv:1508.01982, 2015.
[29] R. Walling and A. Khan, “Characteristics of transformer exciting-current
during geomagnetic disturbances,” Power Delivery, IEEE Trans. on,
vol. 6, no. 4, pp. 1707–1714, 1991.
[30] D. Boteler, “Geomagnetically induced currents: present knowledge and
future research,” Power Delivery, IEEE Trans. on, vol. 9, no. 1, pp.
50–58, 1994.
[31] A. Pulkkinen, E. Bernabeu, J. Eichner, C. Beggan, and A. Thomson,
“Generation of 100-year geomagnetically induced current scenarios,”
Space Weather, vol. 10, no. 4, 2012.
[32] Review of the GMD benchmark event in tpl0071. http://www.
energy.gov/oe/downloads/review-gmd-benchmark-event-tpl-007-1. Accessed: 2016-11-26.
[33] A. Morstad, “Grounding of outdoor high voltage substation: Samnanger
substation,” 2012.
[34] D. K. Molzahn, C. Josz, I. A. Hiskens, and P. Panciatici, “A laplacianbased approach for finding near globally optimal solutions to opf
problems,” IEEE Trans. on Power Systems, vol. 32, no. 1, pp. 305–315,
2017.
[35] H. Nagarajan, M. Lu, E. Yamangil, and R. Bent, “Tightening McCormick relaxations for nonlinear programs via dynamic multivariate
partitioning,” in International Conference on Principles and Practice of
Constraint Programming. Springer, 2016, pp. 369–387.
[36] H. Nagarajan, M. Lu, S. Wang, R. Bent, and K. Sundar, “An adaptive,
multivariate partitioning algorithm for global optimization of nonconvex
programs,” arXiv preprint arXiv:1707.02514, 2017.

