The Minimum Subset Vertex Cover Problem for the Optimal Placement of E₃ HEMP Arrestors

1. Introduction

In the recent years, there has been a growing interest in the mitigation of disturbances due to high-altitude electromagnetic pulses (HEMP) on the power grid [1,2]. The HEMP threat is traditionally split into three components: the $E_1$, $E_2$, and $E_3$. Each component poses their own unique threat to the grid. The $E_1$ and $E_2$ portions are high power, short-time span insults, whereas the $E_3$ portion is a lower frequency, lower power, longer time span threat. All portions of the HEMP disturbance are capable of saturating transformer cores by inducing current through transmission lines and out through ground on the neutral paths of the transformers. A similar phenomenon was seen in the 1989 Hydro-Quebec incident wherein a coronal mass ejection led to geomagnetically induced currents. Transformers were saturated and the grid eventually collapsed [3]. Core saturation can also lead to overheating of transformers and the breakdown of components within transformers. This paper demonstrates how a power grid may be fully protected from an $E_3$ HEMP disturbance using $E_3$ arrestor devices placed at substations that make up any vertex cover of the grid.

Mitigation devices can be categorized under four main archetypes, seen in Figure 1.

The high-side controller sits on the transmission line side of a transformer and pushes geomagnetically induced currents away from the associated substations. The neutral blocking controller is placed on the neutral path to ground of a transformer. It works by rejecting current through the transformer core, preventing magnetic core saturation of the individual transformer. The load-side controller works by drawing power from the $E_3$ HEMP disturbance away from the neutral path to ground and into a storage device on the load-side of the transformer. The final controller archetype is the one being studied in the paper. This is the $E_3$ arrestor and it is a current source that is attached to the bus of a substation. It works by redirecting current from the $E_3$ HEMP insult through safe paths to ground and away from transformer neutral paths to ground. The device is akin to a lightning arrestor, but instead of remaining high impedance at relatively low (60 Hz) frequencies and low at high frequencies, the device is low impedance at low frequency (<1 Hz) and high at 60 Hz. The remainder of this introduction will serve as a literature review on the currently proposed solutions. As it stands, there are no widespread deployed mitigation solutions for the $E_3$ HEMP disturbance problem on North American grids.

Literature on the topic often advocates for blocking capacitors (passive, integral controllers) on transformer paths to ground. In fact, it is typically stated that these neutral blocking capacitors are the optimal device for $E_3$ HEMP mitigation [4,5,6,7,8,9]. The neutral blocking capacitors are a suitable option for geomagnetic disturbance (GMD) mitigation due to the low-frequency nature of the waveform. For a pseudo-DC signal, the blocking capacitor makes an effective open circuit for the disturbance, but allows the [round-precision=0]60 power signals to pass through with minimal impedance. However, it has been shown in [10,11] that the $E_{3A}$ portion of the HEMP waveform contains sufficiently high-frequency components such that the capacitor does not appear to be an open circuit to the disturbance. Additionally, the blocking capacitor adds an integration term to the system. By including this integral term in series with the inductance of the transformer core, the dynamics of the system increase from first order to second order; as such, impulses to the system will generate ringing. This is seen in [10,11], where the $E_{3A}$ portion of the HEMP waveform has sufficiently high frequency to trigger ringing in the system. This phenomenon is called ferroresonance and has been documented in the literature [12,13].

An additional issue with this solution has to do with the phenomenon of GIC-shifting. Capacitors are not dissipative and will not remove HEMP energy from the system [14]. In general, power will only be dissipated by resistive elements and devices performing work. If the energy is not dissipated from the system, it must be moved elsewhere. Often, this is onto other transformers either in the same substations or in neighboring substations. This can exacerbate already present issues on the power grid: protect one transformer at the cost of another [10].

The final issue associated with neutral blocking capacitors is the cost. Beyond a handful of test cases in the US, neutral blocking capacitors are not installed on transformers. As such, thousands of transformers would require upgrades to be compliant with this paradigm. A single ABB SolidGround device, which would protect only one transformer, costs nearly $750k to US taxpayers, adjusted to inflation [15]. Furthermore, the price of one SolidGround unit has not tracked with inflation, since it was initially proposed as a solution. Since it was proposed in 2013, the cost of a single SolidGround unit has increased by over 650% [16], so it is difficult to project the future cost of implementation. Even if only the 6000 most vulnerable transformers are protected by neutral blocking capacitors, it is projected that these upgrades will cost the US taxpayers $4.36B, adjusting for inflation [15].

Another potential solution is the high-side blocking solution [17]. Devices on the high-side include unified power flow controllers (UPFC), in-line transmission level capacitors, and static synchronous compensators (STATCOM). It has been shown that the high-side blocking device can be modified to cooperate at a global level, reducing the required device placements on the grid [18]. Additionally, there is a decentralized approach that has been proven effective in New Zealand against geomagnetic disturbances on the grid [19]. By switching relays, grid operators have been able to redirect unwanted geomagnetically induced current (GIC) around the grid to either be dissipated off through transmission line or absorbed at transformer loads.

The third archetype discussed in this paper is the load-side mitigation device. These devices include energy storage at electrical vehicle charging stations, energy storage at substation upgrade deferral stations, and devices similar to train wayside energy collectors [20]. These devices have the potential to regulate other grid disturbances or are currently being used in the regulation of grid disturbances by storing and supplying power as grid demands change.

The final archetype is the $E_3$ HEMP arrestor. The $E_3$ arrestor is inspired by conventional lightning arrestors already on the grid. A lightning arrestor is placed on the grid in locations vulnerable to lightning strikes. They operate as open circuits until a power surge, such as lightning, is present. The resultant overvoltage creates breakdown in the metal oxide varistor in the arrestor and the resistance drops to near zero. Current from the lightning travels through the arrestor and out to ground, protecting sensitive grid assets. The $E_3$ arrestor works inversely and instead remains high impedance at 60 Hz to prevent the main power signal from pass straight through the bus to ground, but is low impedance at low frequency (<1 Hz).

$E_3$ arrestor placement must also be handled differently than lightning, since a lightning strike is a localized event, whereas HEMP insults are a widespread field over vast swaths of the grid. Lightning strikes at a single point in the grid topology, whereas the $E_3$ field affects large regions of the grid. The combination of complex intertanglings of transmission lines and substations generates large, fractal antennas, capable of receiving power from an $E_3$ HEMP disturbance [21,22]. Since the overvoltaging will appear across the grid, as opposed to in one place, it is necessary to place multiple $E_3$ arrestors throughout the grid. When the HEMP disturbance is detected, actively controlled current sources attached to substation buses will route the $E_3$ portion of the insult out to ground, similar to how a lightning arrestor behaves. The analysis in this paper is done technologically agnostically to allow for grid designers to use the method presented here to assess their own grids. That is, any device of this archetype that meets the specifications laid out in this paper will suffice for mitigation protection. In this paper, the technology agnostic device is a controlled current source.

This paper looks at the optimal placement of $E_3$ arrestor controllers for $E_3$ HEMP mitigation. It has been shown on the neutral and high-side that different placements of archetypal controllers influences the efficacy of the mitigation strategy [18]. For instance, it has been shown that placing neutral blocking devices at critical nodes is an optimal configuration for grid protection against $E_3$ HEMP disturbances [23]. For the purpose of this paper, it is assumed that devices will be added to the system and that the grid will have no pre-existing assets. This serves two purposes. The first is that it allows for better insight into how optimal placement works. The second is that it allows grid designers to verify if sufficient devices are already placed to meet mitigation requirements.

This paper is split into 6 main sections. The first section is the introduction and it provides information on why $E_3$ HEMP disturbances should be mitigated on the grid. This section also discusses current proposed methods for mitigation. The next section presents the systems models. The models in this section are the per-phase transformers model, the 18 and 150-bus synthetic grids, and the $E_3$ HEMP disturbance model. Following that is a short section of the linear quadratic regulator (LQR) and how it is implemented as a global controller for $E_3$ HEMP disturbance mitigation. The fourth section describes the Minimum Subset Vertex Cover algorithm. This is a modified version of the more conventional minimum vertex cover problem. This modified version may only build a minimum vertex cover from a subset of the nodes in the graph. This is because it may not be possible to place an $E_3$ arrestor on all substations for one reason or another. The next section is the results section which shows the outcome of the minimum vertex cover placement on the 18 and 150-bus models. It also provides a comparison against a genetic algorithm for optimal placement. The last section is the conclusion which summarizes the problem and findings. The proposed algorithm in this paper is successful in finding the bounds for placements of energy storage devices to mitigate the $E_3$ HEMP disturbance from the grid. Additionally, it is shown that any vertex cover is sufficient for $E_3$ HEMP mitigation. This provides the ability to analyze the efficacy of pre-existing assets for the mitigation of $E_3$ HEMP disturbances.

2. System Models

The following section describes the system models used within the paper. This include a model of the per-phase transformer, an 18-bus synthetic grid model, a 150-bus synthetic grid model, and the $E_3$ HEMP waveform model [24]. The results from a collection of computational experiments using these models can be found in Section 5.

3. System Models

The following section describes the system models used within the paper. This include a model of the per-phase transformer, an 18-bus synthetic grid model, a 150-bus synthetic grid model, and the $E_3$ HEMP waveform model [24]. The results from a collection of computational experiments using these models can be found in Section 5.

3.1. Per-Phase Transformer Model

The per-phase transformer model used in this paper has been modified from the model presented in [10]. Figure 2 shows the per-phase transformer model used in this paper. The corresponding dynamic system of equations is given below in (1).

$$\left[\begin{array}{c}{L}_1{\dot{i}}_1\\ L_2{\dot{i}}_2\\ {\dot{\lambda }}_m\\ C_B{\dot{v}}_B\end{array}\right]=\left[\begin{array}{cccc}-(R_1+R_c)& -R_c& 0& 1\\ -R_c& -(R_2+R_c+R_L)& 0& 0\\ R_c& R_c& 0& 0\\ -1& 0& 0& -1/R_B\end{array}\right]\left[\begin{array}{c}{i}_1\\ i_2\\ {\lambda }_m\\ v_B\end{array}\right]+\left[\begin{array}{c}-R_c{i}_m\\ -R_c{i}_m\\ -R_c{i}_m\\ i_{tx}+u_B\end{array}\right]$$

(1)

The state vector for the per-phase model is taken as

$$x={\left[\begin{array}{cccc}{i}_1& i_2& {\lambda }_m& v_B\end{array}\right]}^T,$$

where $v_B$ is the high-side bus voltage. Here, $L_1$ and $L_2$ represent the primary and secondary side energy storage due to winding inductance; $R_1$ and $R_2$ are the primary and secondary side winding losses; $R_c$ denotes the core losses; $R_L$ is the resistive load; $v_B$ is the high-side bus voltage; and $C_B$ and $R_B$ are the bus capacitance and resistance, respectively. The $v_B$ term is the bus voltage on the primary side of the transformer and $u_B$ is the controlled current source attached to the bus, acting as the $E_3$ arrestor. The term $i_m$ represents a nonlinear current through the core, approximated as

$$i_m=k_1\tan \left(k_2{\lambda }_m\right)$$

(2)

where $k_1=0.0016$ and $k_2=1.2879$ [10]. The corresponding plot of the current is shown in Figure 3.

Figure 3 demonstrates that as long as a transformer’s magnetic core flux is prevented from saturating, it operates linearly. Within the linearized (unsaturated) region, the ratio between $i_m$ and ${\lambda }_m$ is approximately $1/500$, implying that $L_m\approx 500H$[10], where $L_m$ is the linearized inductance of the transformer core. This linear approximation is used solely for the development of the LQR feedback gains. With a properly chosen Q matrix (i.e. one where the terms in Q outweigh those of R), the transformer magnetic core flux tends to remain low enough such that the transformer remains unsaturated and linear. Thus, with the properly selected LQR gains, the LQR remains an optimal controller and prevents the system from becoming nonlinear and saturated.

3.2. The 18-Bus Synthetic Grid Model

This test case is comprised of 15 transformers across 8 substations and 18 buses. Controllers may be added to any of the buses in the system. This is a synthetic grid that was designed to emulate a real grid [25]. The transformer reduced order model has been presented. In this model, all substations (18 buses) have an $E_3$ arrestor controller. These controllers are controllable current sources that connect directly to a substation’s bus.

3.3. The 150-Bus Synthetic Grid Model

This test case is comprised of 60 transformers across 98 substations and 150 buses. Controllers are limited to the buses which are stated to have a load. This is a synthetic grid that was designed to emulate a real grid [26,27,28]. The transformer reduced order model has been presented. In this model, 90 out of the 98 substations (90 buses) have an $E_3$ arrestor controller.

3.4. The $E_3$ HEMP Waveform Model

The $E_3$ HEMP waveform used in this set of simulations was obtained from the 2021 DOE memo [24]. This waveform is approximately twice the size in magnitude from the 1996 IEC 61000-2-9 standard [29], but matching the waveform magnitude of the 2025 iteration of IEC 61000-2-9. This doubling is used to provide an additional margin of safety when designing mitigation techniques. Traditionally, the HEMP insult has been defined as having three distinct time domains: $E_1$, $E_2$, and $E_3$ [24]. This paper is focused on protecting transformers from the effects of $E_3$. This portion of the HEMP insult is the portion that will send a transformer into saturation for the longest time period [30]. As a result, the risks associated with overheating are more commonly associated with the $E_3$ portion. The $E_1$ portion of the HEMP insult exhibits the highest frequencies and field strengths out of the three time-domain components. The second portion, the $E_2$, is lower in frequency and field strength than the $E_1$ portion of the insult. The final, late-time component, is the $E_3$ portion of the insult. The $E_3$ is comprised of the lowest frequencies and weakest field strengths out of the three components. The late-time portion of the HEMP insult is further divided into two sub-components. The first is the $E_{3A}$, or blast component. The $E_{3A}$ field is described by the equation [24]:

$$E_{3A}\left(t\right)=\alpha (\gamma t-\delta t^2)exp(-t/\beta )$$

(3)

where $\alpha =9.5$, $\beta =1.4$, $\gamma =26$, and $\delta =8.9$. The second component, known as the heave or $E_{3B}$ portion, is defined by the equation:

$$E_{3B}\left(t\right)=-\alpha t^2{\displaystyle \frac{{\beta }^{3}t-3{\beta }^3\gamma +t^4}{\gamma {({\beta }^3+t^3)}^2}}exp(-t/\gamma )$$

(4)

where $\alpha =1.3·{10}^6$, $\beta =200$, and $\gamma =20$. The total field strength is expressed as:

$$E_3\left(t\right)=E_{3A}\left(t\right)+E_{3B}\left(t\right),$$

(5)

where $E_{3A}$ and $E_{3B}$ are defined in (3) and (4), respectively. These parameters are derived from curve fitting based on simulation work conducted by Oak Ridge National Laboratories [24]. The waveform of the insult is depicted in Figure 6.

The apparent voltage along any given transmission line is then given by

$$V\left(t\right)=E_3\left(t\right)·l·cos\left(\theta \right),$$

(6)

where l is the transmission line’s length and $\theta$ is the angle between the transmission line and $E_3$ HEMP field orientation. As a HEMP field becomes more parallel to a transmission line, it couples more. Additionally, voltage scales linearly with the length of the transmission line, with the longest of transmission coupling more. The studies presented in this paper perform a sweep of HEMP fields from 0^o to 180^o [18].

4. The Linear Quadratic Regulator

This paper uses an LQR as the controller for HEMP mitigation. It balances system error against control effort. It should be noted that this is not necessarily an energy storage minimized controller [31]. Further, an LQR only produces an optimal control for linear systems. In this scenario, since the transformers remain unsaturated, the system behaves linearly, and the LQR produces an optimal control [18].

The global LQR (G-LQR) is the full-state feedback controller designed using fullstate knowledge of the grid system. The G-LQR is the absolute most information that the system can have and so it represents the best case scenario for the $E_3$ arrestor. The $E_3$ arrestor has access to all states in the system and the controls can benefit from destructive interference [17]. Let the system dynamics be described by the linear dynamics given in (7) as

$$\dot{\overrightarrow{x}}=A\overrightarrow{x}+B\overrightarrow{u}+D\overrightarrow{v},$$

(7)

where $\overrightarrow{x}$ is the state vector of the local bus, A is the linear dynamics matrix of the bus, B is the control matrix, D is the exogenous input matrix, and $\overrightarrow{v}$ is the time-varying disturbance vector from the HEMP insult. The control vector $\overrightarrow{u}=-K\overrightarrow{x}$, where K is the optimal gains matrix obtained by solving the algebraic Riccati equation for specified Q and R weighting matrices. These weighing matrices are defined based on the grid. These weightings are discussed in Section 5.

5. Graph Theory Approach to Placement

This following section presents a physics informed graph theory method for the placement of $E_3$ arrestor devices. Graph theory interpretations of system dynamics allows for more readily scalable solutions to the $E_3$ HEMP mitigation device placement problem [23]. This section is split into three subsections. The first section discusses the bounding case for placing $E_3$ arrestor mitigation devices. The second provides the physics-based justification for the use of graph theory and the third lays out how the algorithm operates.

5.1. The Bounding Case Solution

Before determining placement algorithms, the bounding case solution must be formed. Bounding case solutions have so far been found for neutral and high side blocking devices [18]. The bounding case solution for HEMP mitigation using the $E_3$ arrestor is an $E_3$ arrestor device on all grid buses. This would permit the entire grid to be controlled by the mitigation devices. Placing $E_3$ arrestors on all buses in the system form a vertex cover of the grid. That is, every transmission line of the system is connected to a substation with at least one $E_3$ arrestor. These bounds are the requirements for the controllers to prevent a HEMP disturbance from having controllability over the system.

However, non-controllability of the HEMP insult is not necessary for HEMP mitigation to occur. Specifically, [32] defines controllability as the ability to transfer some system’s state from an initial state $x_0$ to $x_f$ in a finite period of time, using an unconstrained control vector u. By definition, the HEMP disturbance is a constrained input/control vector on the system. Therefore, the HEMP may not truly have controllability over the system. This naturally raises the question of: what is the minimum required number of $E_3$ arrestor devices to mitigate a HEMP $E_3$ insult, given that the fully protected grid exceeds typical system requirements? It has been found that any vertex cover provides sufficient protection. This is discussed further in Section 5. Therefore, if a vertex cover is necessary for protection, then the vertex cover can be minimized and the number of required $E_3$ arrestor devices can be reduced. Additionally, knowing that a vertex cover is more than sufficient for $E_3$ disturbance mitigation provides an easy check if a grid is already sufficiently protected; if pre-existing $E_3$ arrestor devices form a vertex cover of any sorts on a grid, then the placement is sufficient for mitigation.

It is stated in [18] that the bounding case requirement for transmission line blocking devices is an actuator on each transmission line. The bounding case entirely removes all coupling from the $E_3$ disturbance on the line. A vertex cover comprised of $E_3$ arrestors guarantees that all transmission lines in the system have at least one controller associated with them. Thus, the minimum vertex cover is the minimum number of controllers to the meet the bounding case requirements.

The process for determining the minimum bounding case solution follows. A possible vertex cover solution is the placement of an $E_3$ arrestor at every possible bus location. Doing so forms a vertex cover on the power grid. An $E_3$ arrestor controls bus voltages by routing exogenous power to ground. As such, it is able to directly control the states on adjacent transmission lines, effectively blocking $E_3$ HEMP signals from entering into the transformers. In fact, any vertex cover formed by $E_3$ arrestors can entirely eliminate $E_3$ HEMP coupling from the grid because it controls all of the transmission lines. The next step is to determine where assets already exist on the grid. If pre-existing $E_3$ arrestor assets already form a vertex cover, then the grid has sufficient placements to entirely decouple the exogenous power. If a vertex cover is not formed by pre-existing assets, then the design continues. The list of pre-existing assets are used to seed the generation of a minimum vertex cover (i.e. the smallest number of devices required to form a vertex cover). Additionally, nodes where devices cannot be placed are removed from the set that can form a vertex cover. In the event that no pre-existing assets are on the grid, then the minimum vertex cover is found from the set of nodes where devices may be placed. This minimum vertex cover is the smallest set of nodes where $E_3$ arrestors may be placed to ensure complete decoupling of exogenous power from the disturbance by controlling all transmission lines. This is the minimum bounding case solution.

5.2. Physics-Informed Justification

The power grid may be described using graph theory with transmission lines acting as edges and substations acting as nodes in the network [23,33]. The benefit of this approach is the increased simplicity of the network, while still maintaining knowledge of the core grid topology. This provides the framework for more readily computationally scalable approaches to grid configuration and optimization across a range of applications. This paper uses a physics informed graph theoric algorithm to place technology agnostic energy storage devices on the buses of substations to mitigate the effects on $E_3$ HEMP on the grid.

The approach in this paper is developed around the core concept that individual transmission lines form a complex set of closed loops called cycles that start and end at the same node, without crossing over themselves [34]. These cycles form a complex network of fractal antennas [35], capable of receiving power from HEMP electromagnetic disturbances. Each transmission line acts as an edge within a series of closed loop cycles. The complex interaction of transmission lines forms patches, akin to fractal antennas with edges and holes [21]. This complex interaction allows the grid to absorb power from HEMP insults, akin to other fractal antennas used today to receive RF signals [22]. In the field of RF engineering, fractal antenna impedances, resonances, and gains can be configured through the manipulation of self-similar loops [22,36].

To protect transformers from the $E_3$ HEMP disturbance, there are two commonly studied methods. The first is the breaking of cycles through power rejection. Fractal antennas may be broken with neutral blocking devices [37] on transformer neutral paths within substations or with high-side blocking devices, like relays [19] or unified power flow controllers (UPFCs) [17], on transmission lines. This is one way to tune the power grid as a fractal antenna to reject power from the $E_3$ HEMP disturbance. The alternative $E_3$ arrestor, which is presented in this paper, draws power away from transformers down paths to ground that are safe to receive the impulse [38]. This power can be stored in energy storage devices, like batteries and flywheels, or shunted off as heat [39]. This is equivalent to retuning the grid level fractal antenna to send exogenous power down safe paths to ground. Instead of removing all fractal antennas from the system, the $E_3$ arrestor retunes the existing system to make a new, high-gain antenna to reroute power away from transformers and sensitive power grid assets.

$E_3$ arrestor devices on substation buses can reroute current to ground from every transmission line that they are connected to. For instance, in Figure 4, an $E_3$ arrestor device on substation 1 may draw current from 2 different transmission lines. In effect, exogenous current from those transmission lines is rerouted through as safe path to ground, in lieu of passing through a transformer neutral on substation 1, 2, or 4. This leads to the idea of using vertex covers to prevent exogenous power from flowing through transformer neutral paths to ground. By placing $E_3$ arrestor devices at substations that form a vertex cover C, every transmission line (graph edge) is incident to at least one controlled substation (vertex). In graph-theoretic terms, the complement $V\setminus C$ is an independent set, meaning there are no edges between any two vertices in $V\setminus C$. Therefore, no transmission line exists whose endpoints are both uncontrolled; every edge has at least one endpoint with an available arrestor actuator. In the physical interpretation adopted here, this guarantees that each line has an adjacent controlled bus capable of providing a low-frequency shunt path, enabling rerouting of quasi-DC disturbance-driven current away from transformer neutral paths and toward designated safe paths to ground (or energy storage), subject to actuator limits and controller tuning.

5.3. Minimum Subset Vertex Cover Algorithm

The minimum vertex cover of a network determines the smallest set of nodes adjacent to every connection [40,41]. More specifically, given an undirected graph $G=(V,E)$ where V denotes the set of vertices (substations) and E denotes the set of edges (transmission lines), a vertex cover is a set $C\subseteq V$ such that for every edge $e=(u,v)\in E$, at least one endpoint of the edge is in C, i.e.,

$$\forall (u,v)\in E,u\in C\mathrm{or}v\in C.$$

Equivalently, in indicator-variable form,

$$\forall (u,v)\in E,x_u+x_v\ge 1.$$

While computing an arbitrary vertex cover is trivial, computing the minimum vertex cover is expensive and NP-hard.

The minimum vertex cover can be found by solving a computationally expensive integer linear program (ILP)

$$\begin{array}{ccc}\hfill \min & & \sum _{x\in V}{x}_v\hfill \\ \hfill \mathrm{st}& & x_u+x_v\ge 1,\forall (u,v)\in E\hfill \\ \hfill & & x_v\in \{0,1\},\forall v\in V\hfill \end{array}$$

(8)

In practical settings, arrestors may be installable only on a subset of substations (or buses). Let $V_a\subseteq V$ denote the set of allowable vertices for placement. The minimum subset vertex cover problem is then

$$\begin{array}{cc}\hfill \min & \sum _{v\in V_a}{x}_v\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& x_u+x_v\ge 1,\forall (u,v)\in E\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & x_v\in \{0,1\},\forall v\in V_a\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & x_v=0,\forall v\in V\setminus V_a.\hfill \end{array}$$

(12)

A feasible solution exists if and only if $V_a$ itself is a vertex cover of G, i.e., every edge has at least one endpoint in $V_a$. Here, each vertex v corresponds to a binary variable $x_v$. If $x_v$ is one, then it is included in the vertex cover and when it is zero, it is not.

For a large graph, solving this formulation can be computationally prohibitive, so an approximation must be found. A factor-n approximation denotes that the size of the approximate vertex cover is no more than n-times the optimal size. It has been shown that the best polynomial time approximation algorithm to the minimum vertex cover has factor-1.3606 [42].

A much more efficient factor-2 approximation can be found by relaxing the above ILP into a linear program (LP)

$$\begin{array}{ccc}\hfill \min & & \sum _{x\in V}{x}_v\hfill \\ \hfill \mathrm{st}& & x_u+x_v\ge 1,\forall (u,v)\in E\hfill \\ \hfill & & 0\le x_v,\forall v\in V\hfill \end{array}$$

(13)

A vertex v is included in the vertex cover when $x_v\ge {\displaystyle \frac{1}{2}}$ [43]. In fact, an extreme point solution to this formulation, which is a solution that cannot be expressed as a convex combination of two other feasible solutions, possesses a half-integral property. This means that each $x_v$ must take the value of $\{0,{\displaystyle \frac{1}{2}},1\}$. Such extreme point solutions can be found by solving the LP formulation using the simplex method.

The factor can be further reduced using additional approximation methods. One method to reduce the factor to 1.5 uses 4-Color Theorem (4CT) to partition vertices. The LP formulation above is a 2-factor approximation because, as a relaxation to the ILP, its objective value must be less than or equal to the ILP. If we assume that we find an extreme point solution, rounding the variables with value ${\displaystyle \frac{1}{2}}$ to 1 at most doubles the objective value. Call this new objective value FIP. Since $LP\le ILP$ then $FIP\le 2LP\le 2IP$, which gives the 2-factor. This bound can be improved by reducing the number of variables that need to be rounded.

A graph is called planar if it can be drawn on a plane where none of the edges intersect outside its vertices. This situation occurs commonly with networks like power grids. Next, a coloring of a graph labels each vertex with a color, so that no two adjacent vertices share the same color. In a planar graph, the Four Color Theorem (4CT) [44] guarantees that a coloring of the graph exists with no more than four colors. In the context of a set cover, all vertices of a single color can be removed from consideration and a set cover will still exist since no edge is adjacent to two vertices of the same color. This information can be combined with the LP formulation in the following manner to form the 4CTLP approximation.

1.

Color the planar graph with at most four colors.

2.

Solve the LP formulation above using the simplex method in order to achieve an extreme point solution.

3.

Determine the color whose vertices contain the most number of fractional, ${\displaystyle \frac{1}{2}}$, values. For each vertex, v, in this color, set each of the fractional $x_v$ variables to zero.

4.

Set the remaining fractional $x_v$ variables from the other colors to 1

5.

All vertices with value $x_v=1$ are included in the vertex cover. The remaining are not.

Since the above process selects the color with the most fractional values and there are only four colors, there can not be any more than ${\displaystyle \frac{3}{4}}$ of the remaining values that need to be rounded to 1. This increases the objective value by at most $2{\displaystyle \frac{3}{4}}={\displaystyle \frac{3}{2}}$, which gives a factor-3/2 approximation.

From practical or heuristic reasons, it may be desirable to further restrict the number of vertices allowed the vertex cover. For example, this paper uses the vertex cover to determine which buses should be protected from an HEMP insult and not all buses can install protection. The set of buses that can install protection can be represented as a subset of all vertices, $V^{\prime }\subseteq V$. As long as the remaining allowable vertices are incident to all edges in the graph, a vertex cover still exists. This paper uses the term Minimum Subset Vertex Cover when a minimum vertex cover from a subset $V^{\prime }\subseteq V$ is required. The above algorithms can be easily modified to accommodate this change by restricting the vertices to $V^{\prime }$ rather than V.

It is possible to use this minimum subset vertex cover to verify if pre-existing assets meet the protection criterion. In fact, it is a linear time check to verify if pre-existing assets meet the protection criterion. If the pre-existing assets form a vertex cover, then sufficient energy storage sizing and individual controller specifications will guarantee protection against $E_3$ HEMP.

5.3.0.1. Linear-time verification of a vertex cover.

Given a candidate protected set $C\subseteq V$, verification that C is a vertex cover is $O\left(\right|E\left|\right)$:

1.

For each edge $(u,v)\in E$:

2.

If $(u\notin C)$ and $(v\notin C)$, then C is not a vertex cover; terminate.

3.

If all edges pass the test, then C is a vertex cover.

To solve this, the following steps are taken:

1.

Color the graph using 4 colors such that no neighboring nodes share a common color.

2.

Partition the vertices into four subsets, $V_1,V_2,V_3,V_4$, determined by color.

3.

Perform the ILP algorithm, minimizing ${\sum }_{v\in V}c\left(v\right)x_v$, subject to the constraints $x_u+x_v\ge 1,\forall \{u,v\}\in E$ and $x_v\in \{0,1\},\forall v\in \eta$. If a vertex is not in the permissible set ($v_i\notin \eta$), then its upper and lower bounds are set to 0.

4.

Identify the fractional vertices from the optimization problem

5.

Find the color, $C_i$, with the most fractional values.

6.

Form the vertex cover S, including all vertices where $x_v=1$ and $x_v=0.5$, except the vertices in the color class $C_i$

The resulting vertex cover will be an approximation of the true minimum vertex cover, comprised only of members of the permissible set $\eta$. Figure 7 shows the Minimum Subset Vertex Cover algorithm as a flowchart.

6. Results

The following case studies mitigate a HEMP insult using the global LQR scheme described in section 3 and presented more completely in [18]. This is an optimal controller in the linear regime of a system. The global controls allow for cooperation between controllers to bring error states down lower. The A and B matrices are estimated using MATLAB’s built-in $\mathtt{jacobianest}$ function. LQR gains are found using MATLAB’s built-in $\mathtt{lqry}$ function. The objective of all controllers is to minimize the maximum per-unit flux of all transformers. However, it is a permissible solution for placements and control laws to produce an absolute maximum transformer flux of $|{\lambda }_m|<1.2pu$ across all transformers. This constraint prevents saturation of the transformer magnetic core.

6.1. 18-Bus Model

For the Horton et al. model, it is assumed that energy storage controls may be placed on any bus within the system. These controllers are technology agnostic and provide a control effort in the form of current. The Q and R matrices for this were set to multiples of the identity, $Q={10}^3{\mathcal{I}}_{n\times n}$, where n is the number of transformers in the system and $R={\mathcal{I}}_{m\times m}$, where m is the number of controllers. Simulations were truncated to running only for 4 seconds. This was decided since all transformers would reach peak $|{\lambda }_m|$ values within the first 4 seconds of the simulation. This drastically reduced run times for the all-permutations, graph theory, and genetic algorithm (GA) test cases. All test cases are exposed to a sweep of 16 different HEMP field orientations, ranging from 0^o to 180^o [18].

6.1.1. The All-Permutations Comparison

This section provides the baseline test case to compare the optimization techniques against. These results were found by running all possible permutations of controller configurations and recording the best performing configurations for each controller count. Results are shown in Table 1. Results in the all-permutations case are not subject to any actuator constraints and so future current and power limits may alter the optimization schemes.

6.1.2. The Minimum Subset Vertex Cover Results

Since all buses in the system may contain a controller, the 4CTLP algorithm is used to find a vertex cover of all vertices in the connectivity graph. This algorithm found that 5 out of 18 substation buses must be protected. Through experimentation, it was found that, for the sake of the vertex cover problem, as long as one bus was protected at a substation, the substation was considered protected. The best performing case yielded a maximum absolute flux of $|{\lambda }_m|=0.3472pu$. This solution was found in approximately 8 ms and verified in 706 ms with a dynamic simulation, putting the total runtime of the graph theory algorithm as 714 ms. The graph theory approach was performed on a single thread of a 12th Gen Intel(R) Core(TM) i7-12700 @ 2.10 GHz. For the all-permutations case, the script was provided with 256 threads on the Sandia National Labs high-performance computing cluster. The cores used in the simulation are the 32 socket Intel Xeon Platinum 8176 @ 2.10GHz base (3.8GHz turbo) chip. The cores were also hyperthreaded, with 2 threads per core.

For comparison, a GA optimizer was used to place the mitigation devices, as well. The GA was run on the Sandia National Labs high-performance computing cluster. The cluster is comprised of 32 socket Intel Xeon Platinum 8176 @ 2.10GHz base (3.8GHz turbo) chips. The GA used 256 threads running in parallel. The runtime of this GA was 73.26 minutes, or approximately 6100× the runtime of the graph theory based approach. Figure 8 presents the Pareto front for the placements generated by the GA. The GA found a solution that required 4 mitigation devices and the minimum vertex cover required 5 mitigation devices. This controller count from the minimum vertex cover approximation agrees with the findings of the full NP-hard computation, which found that minimum vertex cover comprised of 5 nodes. The approximation agrees with the NP-hard problem, providing evidence to the claim that the approximation will provide at most $3/2\times$ the number of nodes in the minimum vertex cover as the NP-hard problem. While the GA found a more optimal solution to the problem, the Minimum Subset Vertex Cover results provide important insight into the problem, as well as a fast check to determine protection of a grid using pre-existing assets (i.e. checking if pre-existing assets is a linear scaling problem). This is discussed further in the 150-Bus Model subsection of the Results.

6.2. 150-Bus Model

The 150-bus model differs from the 18-bus Horton et al. model in that not all buses are stated to have a controller. For the sake of this problem, it is assumed that each substation on the grid may be equipped with a controllable energy storage device on at least one of its buses. Due to the number of possible placements of the $E_3$ arrestor, the all-permutations verification case is omitted ($2^{90}$ possible combinations). For the same reasons, the GA is omitted from this section. Once again, the nature of energy storage device actuator is technology agnostic. The Q and R matrices for this were set to $Q={10}^{15}{\mathcal{I}}_{n\times n}$, where n is the number of transformers in the system and $R={\mathcal{I}}_{m\times m}$, where, m is the number of controllers. While the weighting of the Q matrix may seem high, it has been shown in other EMP mitigation work that a large ratio of values in Q to R in LQR controls tends to yield an approximate energy storage optimized controller [11,37,45]. For the 150-bus model, it is stated for the sake of this example that 90 out of the total 150 buses may be controlled. Therefore, the minimum subset vertex cover is found on 90 of the 98 member long V set. Using the 4CTLP algorithm, it was found that the 150-bus model may be protected with no more than 57 devices. The Minimum Subset Vertex Cover was found in 486 ms of runtime and the dynamic simulation was verified in 40.56 s.

Remark 1. 

Initial studies show that any vertex cover is capable of sufficiently protecting the grid from $E_3$ HEMP. A selection of 250 vertex covers from the permissible set η were generated at random. These vertex covers ranged from 58 to 72 members in size. It should be noted that the minimum vertex cover found previously contained 57 members. In every test case, the maximum $\parallel {\lambda }_m\parallel <<1.2$ pu. In fact, the worst $\parallel {\lambda }_m\parallel$ from this set was $6.244·{10}^{-5}$ pu. For sufficiently tuned controllers, any vertex cover is sufficient for protection of the grid from $E_3$ HEMP disturbances since all transmission lines are controlled. This is due to the fact that placing controllers at each node in a vertex cover of a grid ensures that power absorbed by the fractal antenna structure of the grid is rerouted through the $E_3$ arrestor devices, effectively retuning the fractal antenna structure of the grid, as discussed in Section 4.1.

While finding a minimum vertex cover is an NP-hard problem, verifying if a given set of vertices is a vertex cover is a linear problem that scales with $\mathcal{O}\left(\right|V|+|E\left|\right)$. If pre-existing assets on the grid are considered in this problem, then it becomes a simple linearly scaled check to verify that the $E_3$ arrestor device placements are sufficient for $E_3$ HEMP mitigation. That is, if the pre-existing assets form a vertex cover on the grid, it can be easily verified with minimal computational effort if new devices must be added to the grid to mitigation $E_3$ HEMP disturbances. Even if a vertex cover is not achieved by pre-existing devices, the grid may still be protected against $E_3$ HEMP disturbances. More work needs to be carried out on the topic of optimizing $E_3$ arrestor device placement, as well as on the topic of incorporating $E_3$ arrestor devices with additional controller types, such as the high-side and neutral blocking devices.

7. Conclusions

This paper presents a novel concept of $E_3$ arrestors analogous to lightning arrestors. It was found that any vertex cover placement of $E_3$ arrestors can almost completely remove all disturbance from an $E_3$ HEMP waveform applied to the power grid since all transmission lines are controlled. Therefore, this paper contends that the minimum number of controllers to completely mitigate the $E_3$ HEMP waveform (i.e. meet the high-side bounding requirements in [18]) is the minimum vertex cover. Since the minimum vertex cover is prohibitively expensive to solve and most power grids can be represented by planar graphs, this paper uses an approximation algorithm based on the four-color theorem combined with an LP approximation to find a 3/2-factor approximation to the minimum vertex cover. If pre-existing $E_3$ arrestor assets can be utilized for $E_3$ HEMP mitigation, then no additional devices may be required. It has also been shown that this placement algorithm of mitigation devices is many orders of magnitude faster than using a GA to optimize the placements, while still yielding similar results.

The main drawback to this approach is that it does not tell the true minimum required number of controllers to maintain the flux of all transformers such that $\left|\right|{\lambda }_m\left|\right|<1.2pu$. Since there is tolerance in some level of DC flux on the transformers, there is room to further optimize the controller placements and use these results as a homotopy for a scalable optimizer. Additionally, future work should address the decentralized control problem. This paper assumes that the controllers have fullstate knowledge of the system. In practice, the mitigation strategy would become significantly more robust with a decentralized control scheme. However, decentralized controls schemes typically under-perform compared to centralized control schemes. This future work would look at optimizing between sensor placement, controller placement, computational power associated with state estimation, control efforts, and energy storage requirements.