A Hands-On Python Example
Geomagnetically Induced Currents (GIC) are quasi-DC currents driven into power transmission networks during geomagnetic storms. When the Sun ejects a coronal mass ejection toward Earth, the resulting disturbance in the geomagnetic field induces an electric field at ground level. Long transmission lines act like giant antennas for this field, and the induced currents flow into the earth through transformer neutral grounding points. Because power transformers are designed for 50/60 Hz AC, this quasi-DC current can saturate the transformer core half-cycle by half-cycle, causing overheating, harmonic distortion, and — in the worst historical cases (Quebec 1989, Halloween Storm 2003) — cascading blackouts.
Utilities mitigate this risk by installing series blocking devices (capacitors or resistors) in the transformer neutral-to-ground path at selected substations. Since blocking devices are expensive, the practical engineering question is:
Given a limited budget of k blocking devices, which substations should receive them so that the worst-case GIC across every possible storm direction is minimized?
This article walks through a complete, runnable example: modeling a small transmission network with the classic Lehtinen–Pirjola (LP) method, then solving the placement problem as a robust (minimax) combinatorial optimization.
1. The Physics: Lehtinen–Pirjola Network Model
A transmission network is modeled as a resistive DC circuit. Each substation i has a transformer neutral grounding resistance $R_{g,i}$, and each transmission line $(i,j)$ has a series resistance $R_{ij}$ carrying an induced EMF $E_{ij}$ caused by the geoelectric field.
The nodal voltage equation (Kirchhoff’s Current Law at every substation) is:
$$
(\mathbf{Y} + \mathbf{Y}_e),\mathbf{V} = \mathbf{J}
$$
where $\mathbf{Y}$ is the network admittance matrix built from line conductances $g_{ij} = 1/R_{ij}$, $\mathbf{Y}e = \mathrm{diag}(1/R{g,i})$ is the grounding admittance, and $\mathbf{J}$ is the current injection vector caused by the storm’s electric field:
$$
J_i = \sum_{j \in \mathcal{N}(i)} g_{ij}, \big(\mathbf{E}\cdot \mathbf{L}_{ij}\big)
$$
Here $\mathbf{E}$ is the uniform horizontal geoelectric field vector (V/km) and $\mathbf{L}_{ij}$ is the displacement vector (in km) from node i to node j.
Once $\mathbf{V}$ is solved, the GIC flowing through each transformer’s neutral (the quantity that damages equipment) is simply:
$$
I_{g,i} = \frac{V_i}{R_{g,i}}
$$
A blocking device is modeled by setting $R_{g,i}$ to an effectively infinite value, which drives $I_{g,i}$ toward zero at that node.
The Key Trick: Linearity in Storm Direction
If the storm direction is $\theta$, then $\mathbf{E} = E_0(\cos\theta, \sin\theta)$, and since $\mathbf{J}$ is linear in $\mathbf{E}$:
$$
\mathbf{J}(\theta) = \cos\theta \cdot \mathbf{A} + \sin\theta \cdot \mathbf{B}
$$
where $\mathbf{A}$ and $\mathbf{B}$ are constant vectors (independent of $\theta$). Because the system matrix $(\mathbf{Y}+\mathbf{Y}_e)$ does not depend on $\theta$ either, the solution is also linear in $\theta$:
$$
\mathbf{V}(\theta) = \cos\theta \cdot \mathbf{V}_A + \sin\theta \cdot \mathbf{V}_B
$$
This means we never need to re-solve the linear system for every storm angle — we solve it twice (for $\mathbf{A}$ and $\mathbf{B}$) and reconstruct every angle by simple trigonometric scaling. This is the speed optimization used below.
2. The Concrete Example
We model an 8-substation regional grid (a ring topology with a central hub, which is a common real-world structure). We test every 1-degree storm direction (360 angles) and evaluate every possible placement of 2 blocking devices ($\binom{8}{2} = 28$ combinations), selecting the placement that minimizes the worst-case GIC across all storm directions and all unprotected substations.
3. Python Implementation (Google Colab Ready)
1 | # ============================================================ |
4. Code Walkthrough
Section 1 — Topology: Eight substations are placed on a 2D plane (in kilometers) and connected as a ring with a central hub, which mirrors how many regional transmission networks are actually built for redundancy. Line resistance is computed from geographic distance using a realistic resistance-per-kilometer figure for extra-high-voltage lines.
Section 2 — Grounding parameters: Every transformer neutral starts with a small, realistic grounding resistance (0.6 Ω). A “blocked” node is simulated by replacing that resistance with an enormous value (1 MΩ), which forces near-zero current through that neutral — exactly what a physical blocking capacitor does for DC/quasi-DC current.
Section 3 — Admittance matrix construction: This builds the graph Laplacian-style matrix $\mathbf{Y}$ from line conductances, plus two auxiliary vectors $\mathbf{A}$ and $\mathbf{B}$ that capture how much current each node “collects” per unit of eastward or northward geoelectric field, respectively. These vectors are computed once and reused for every storm direction.
Section 4 — The fast solver: Instead of solving the linear system once per storm angle, we solve it exactly twice (np.linalg.solve with a 2-column right-hand side) to get $\mathbf{V}_A$ and $\mathbf{V}_B$. Every angle’s voltage vector is then reconstructed instantly via cos(θ)·V_A + sin(θ)·V_B, vectorized across all 360 angles simultaneously with np.outer.
Section 5 — Robust placement search: For each of the 28 possible 2-node placements, the code computes the worst-case (maximum absolute) GIC across all substations and all 360 storm directions, then keeps the placement with the smallest worst case — a classic minimax robust design.
Section 6 — Visualization: Builds one combined figure with four panels described below.
5. Why This Runs Fast (Speed Optimization Explained)
A naive implementation would loop over all 360 storm angles inside the loop over 28 placements, rebuilding and re-solving an 8×8 linear system 10,080 times. Because the system matrix $(\mathbf{Y}+\mathbf{Y}_e)$ does not change with storm angle — only the right-hand side does — we exploit the linearity shown in Section 1 of the physics discussion to solve the system only twice per placement (56 solves total instead of 10,080), then reconstruct all 360 angles algebraically. For this small 8-node network the difference is invisible in wall-clock time, but the technique matters enormously at realistic scale: a continent-scale grid model with hundreds of substations and thousands of candidate placements would make the naive approach computationally prohibitive, while this vectorized approach scales linearly with the number of placements rather than the product of placements and angular resolution.
6. Execution Results
Number of candidate placements evaluated : 28 Angular resolution : 360 steps Total computation time : 0.0177 s Baseline worst-case GIC (no mitigation) : 346.79 A Best blocking-device placement (nodes) : (2, 7) Worst-case GIC after mitigation : 323.44 A Risk reduction : 6.7 %
7. Interpreting the Graphs
Top-left (3D surface, baseline): This shows the magnitude of GIC at every substation (y-axis) as the storm direction sweeps through 360 degrees (x-axis), with no mitigation installed. The tall ridges reveal which substations are most exposed, and at which storm orientation — typically substations aligned with the dominant line direction see the highest peaks, since the induced EMF is strongest when the geoelectric field is parallel to the transmission line.
Top-right (3D surface, mitigated): The same surface after installing blocking devices at the optimal two substations. The peaks at the protected nodes are flattened to near zero, and — importantly — the peaks at the remaining nodes are also reduced, because blocking a node changes the current distribution throughout the whole network, not just at that node.
Bottom-left (worst-case curve): This directly compares, for every storm direction, the single highest GIC value anywhere in the network, before (red) and after (green) mitigation. The dashed horizontal line marks the guaranteed worst-case ceiling achieved by the optimal placement — this is the number a grid operator would actually report as “maximum expected GIC under any storm orientation.”
Bottom-right (network map): A geographic view of the grid, with node size proportional to its worst-case GIC after mitigation and the two selected blocking-device substations highlighted in red. This panel is the most operationally useful output: it tells an engineer exactly where to schedule equipment installation.
8. Takeaways for Real-World Grid Operators
This example uses a small, illustrative network, but the underlying method — the Lehtinen–Pirjola model combined with a linearity-exploiting solver and a robust minimax search over storm direction — is exactly the kind of approach used in academic and utility-industry studies of GIC mitigation. For larger real-world grids (hundreds of substations), exhaustive combination search becomes impractical once the number of devices $k$ grows, and the same linear-algebra speedup would typically be paired with a greedy or integer-programming heuristic rather than brute-force enumeration. Nonetheless, the physics, the math, and the vectorization trick shown here remain the foundation of that larger-scale analysis.
