# Geomagnetic Disturbance Uncertainty Quantification Modeling: An Electromagnetic Transient and Steady-State Simulation based Approach

Arturo Bretas, Kiana Brown, Varoozhan Hartoonian, Jonathan Snodgrass

► To cite this version:

Arturo Bretas, Kiana Brown, Varoozhan Hartoonian, Jonathan Snodgrass. Geomagnetic Disturbance Uncertainty Quantification Modeling: An Electromagnetic Transient and Steady-State Simulation based Approach. 2023 IEEE PES Innovative Smart Grid Technologies Europe (ISGT EUROPE), Oct 2023, Grenoble, France. pp.1-5, 10.1109/ISGTEUROPE56780.2023.10407839. [hal-04429149](https://hal.science/hal-04429149v1)

HAL Id: hal-04429149

https://hal.science/hal-04429149v1

Submitted on 31 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Geomagnetic Disturbance Uncertainty Quantification Modeling: An Electromagnetic Transient and Steady-State Simulation based Approach

Arturo Bretas*‡, Kiana Brown*, Varoozhan Hartoonian*, Jonathan Snodgrass†

* Pacific Northwest National Laboratory, 99354 Richland, USA

† Texas A&M University, 77843 College Station, USA

‡ G2Elab, Grenoble INP, CNRS, Université Grenoble Alpes, 38000 Grenoble, France

arturo.bretas@pnnl.gov, kiana.pitman@pnnl.gov, varoozhan.hartoonian@pnnl.gov, snodgrass@tamu.edu

**Abstract**—Geomagnetic disturbances have been shown to disrupt the operation of the bulk electrical system through low-frequency effects in the earth's magnetic field that in turn induce changing electrical fields on the earth's surface. As a result, geomagnetically induced currents flow in transmission lines, introducing the risk for widespread damage to high-voltage transformers and voltage collapse due to induced reactive power loss. In a previous work, the authors showed how certain modeling assumptions can lead to results that are unacceptable in certain edge case scenarios. In this paper, the authors build on their previous results using a probabilistic approach for uncertainty quantification in order to compare the results from electromagnetic transient program modeling software, ATP, to the positive-sequence calculation software, PowerWorld.

**Index Terms**—electromagnetic transient program, geomagnetic disturbances, geomagnetically induced currents, positive-sequence calculation, uncertainty quantification modeling.

## I. INTRODUCTION

As observed through geomagnetic disturbance (GMD) events, GMD can severely affect critical infrastructure and its operation through induced quasi-direct currents [1]. The power industry has been working towards understanding the impact of geomagnetically-induced currents (GICs) on the grid and in the development of tools to help with mitigation over the last decade. This work has led to the creation of a North American Electric Reliability Corporation (NERC) Standard TPL-007-4 [2], which is being implemented and has led to the developed GMD assessment processes that consider physics-based models.

When considering a physics-based model, one must also consider uncertainties within the model [3]. This uncertainty can be quantified [4], and analyzed statistically [5] in order to guide the decisions regarding whether to accept certain modeling assumptions. Several GIC-oriented modeling tools have already been implemented, including harmonic generation and propagation, transformer thermal models, commercial power flows [6], and dynamic simulations [7]. This wide array of tools is crucial for modeling the power system to ensure safe and reliable operation, but it is also crucial to know which tools are appropriate to use in particular circumstances and to understand when a tool may not provide the necessary resolution. The latter includes situations due to individual

modeling uncertainties such as operating condition, transformer saturation, grounding resistance, and other system characteristics which should all be appropriately considered [8]–[10].

In this work, implicit modeling uncertainties are quantified through a probabilistic approach for the purpose of understanding the advantages, limitations, and impacts of different types of models used for GMD events. Further, statistical analysis is performed on GMD synthetic data, generated using the Alternative Transient Program (ATP) [11]. Data from ATP is then compared with results generated by the modeling software PowerWorld [12] for similar scenarios, and modeling differences are statistically analyzed.

The specific contributions of this work towards the state-of-the-art are:

* GMD implicit modeling uncertainty quantification.

* Statistical analysis of GMD synthetic data.

* Error quantification between different modeling tools.

The remainder of the paper is organized as follows. Section II provides theoretical background on uncertainty quantification and model calculations. Section III presents the conceptual framework of this research. A case study based on the IEEE 9-bus is shown in Section IV. Finally, Section V presents concluding remarks from this work.

## II. THEORETICAL BACKGROUND

### A. Uncertainty Quantification

Uncertainty propagation can be understood as the quantification of uncertainties in the system outputs from uncertain models or inputs. Modeling uncertainty can be either structural/topological or parametric [13]. If one assumes that there is no numerical uncertainty, or uncertainty in system inputs, then one can assume that the source of uncertainty in outputs are from the model.

Considering uncertainty quantification, there are two types of approaches [14]: forward propagation and inverse assessment. The following will summarize the first approach, which was leveraged for this work.

The forward propagation approach for uncertainty quantification focuses on the effect of the outputs from the variability listed in the sources of uncertainty. The aim of the analysis can be:

1) Evaluate mean and variance of outputs.

<footnote>This work was funded by the U.S. Department of Energy under Contract DE-AC05-76RL01830.</footnote>
2) Evaluate reliability of the outputs.

3) Assess the probability distribution of the outputs.

The arithmetic mean, is equal to the sum of all of the samples divided by the number of samples. This can be mathematically expressed by:

$$ \bar{x} = \frac{1}{n} \left( \sum_{i=1}^{n} x_i \right) = \frac{x_1 + x_2 + \cdots + x_n}{n}. \quad (1) $$

where, $\bar{x}$ is the arithmetic mean of samples $x_1, x_2, \dots, x_n$ and $n$ is the number of samples. The variance of a collection of $n$ equally likely values can be otherwise mathematically expressed by:

$$ \sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2. \quad (2) $$

were $\sigma^2$ is the variance and $\sigma$ is the standard deviation. With these two metrics available, reliability of the outputs can be provided through statistical hypothesis testing [15].

Considering reliability evaluation of outputs, the classical formulation towards such is the quasi-static state estimation problem, modeled through the Weighted Least Squares (WLS) method [14]. In this classical formulation, the system is modeled as a set of non-linear algebraic equations, known as the measurement model:

$$ \mathbf{z} = h(\mathbf{x}) + \mathbf{e}. \quad (3) $$

where $\mathbf{z} \in \mathbb{R}^{1 \times d}$ is the measurement vector, $\mathbf{x} \in \mathbb{R}^{1 \times N}$ is the vector of state variables, $h: \mathbb{R}^{1 \times N} \to \mathbb{R}^{1 \times d}$ is a continuously non-linear differentiable function, and $\mathbf{e} \in \mathbb{R}^{1 \times d}$ is the measurement error vector. Each measurement error, $e_i$ is assumed to have zero mean, standard deviation $\sigma_i$ and Gaussian probability distribution. $d$ is the number of measurements and $N$ is the number of states.

In the classical WLS approach, the best estimate of the state vector in (3) is found by minimizing the weighted norm of the residual, or the cost function $J(\mathbf{x})$:

$$ J(\mathbf{x}) = \|\mathbf{z} - h(\mathbf{x})\|_{R^{-1}}^2 = [\mathbf{z} - h(\mathbf{x})]^T \mathbf{R}^{-1} [\mathbf{z} - h(\mathbf{x})] \quad (4) $$

where $\mathbf{R}$ is the covariance matrix of the measurements. In [5], it is shown that in the gross error analytic process, a two-step approach should be adopted. In the first step, all measurements should be weighted equally proportional to the measurement magnitude and gross error analysis is performed. After gross error processing, in the second step, meter precision can be restored and state estimation can be executed.

Regarding the solution of (4), it is obtained through the Newton-Raphson method. The linearization of (3) equals: In this paper, we consider the standard deviation of each measurement to be 1% of the measurement magnitude, which has been shown to improve the detection of bad data [5]. In order to solve this problem, (3) is linearized at a certain point $\mathbf{x}^*$ in (5) and the optimal states are found through an iterative process.

$$ \Delta \mathbf{z} = H \Delta \mathbf{x} + \mathbf{e} \quad (5) $$

where $H = \frac{\partial h}{\partial \mathbf{x}}$ is the Jacobian matrix of $h$ at the current state estimate $\mathbf{x}^*$, $\Delta \mathbf{z} = \mathbf{z} - h(\mathbf{x}^*) = \mathbf{z} - \mathbf{z}^*$ is the correction of the measurement vector and $\Delta \mathbf{x} = \mathbf{x} - \mathbf{x}^*$ is the correction of the state vector. The WLS solution can be seen,

geometrically [16], as the projection of $\Delta \mathbf{z}$ onto the Jacobian space by a linear projection matrix $P$, i.e. $\Delta \hat{\mathbf{z}} = P \Delta \mathbf{z}$. Letting $\mathbf{r} = \Delta \mathbf{z} - \Delta \hat{\mathbf{z}}$ be the residual vector, the $P$ matrix that minimizes $J(\mathbf{x})$ will be orthogonal to the Jacobian range space and to $\mathbf{r}$; $\Delta \hat{\mathbf{z}} = H \Delta \hat{\mathbf{x}}$. This is in the form:

$$ \langle \Delta \hat{\mathbf{z}}, \mathbf{r} \rangle = (H \Delta \hat{\mathbf{x}})^T R^{-1} (\Delta \mathbf{z} - H \Delta \hat{\mathbf{x}}) = 0. \quad (6) $$

Solving (6) for $\Delta \hat{\mathbf{x}}$:

$$ \Delta \hat{\mathbf{x}} = (H^T R^{-1} H)^{-1} H^T R^{-1} \Delta \mathbf{z}. \quad (7) $$

At each iteration, a new incumbent solution $\mathbf{x}_{new}^*$ is found and updated following $\mathbf{x}_{new}^* = \mathbf{x}^* + \Delta \hat{\mathbf{x}}$. (7) is solved at each iteration until $\Delta \hat{\mathbf{x}}$ is sufficiently small to claim convergence of the solution.

The projection matrix $P$ is the idempotent matrix that has the following expression:

$$ P = H(H^T R^{-1} H)^{-1} H^T R^{-1}. \quad (8) $$

As shown in [15], the geometrical position of the measurement error in relation to the range space of $H$ provides another way of interpreting the state estimation. Hence, as the measurement vector can be decomposed into two subspaces, and thus it is possible to decompose the measurement error vector into two components, as follows:

$$ \mathbf{e} = \underbrace{P \mathbf{e}}_{\mathbf{e}_U} + \underbrace{(I - P) \mathbf{e}}_{\mathbf{e}_D}. \quad (9) $$

The component $\mathbf{e}_D$ is the detectable error, which is the residual in the classical WLS model, while the component $\mathbf{e}_U$ is the undetectable error, also known as the masked component of the error. $\mathbf{e}_D$ is in the orthogonal space to the range space of Jacobian whereas $\mathbf{e}_U$ is hidden in the Jacobian space.

$$ \mathbf{e}^2 = \mathbf{e}_D^2 + \mathbf{e}_U^2. \quad (10) $$

The error vector $\mathbf{e}$ in (10) is called the Composed Measurement Error (CME). In order to quantify the undetectable error, the Measurement Innovation Index (II) is introduced [17]:

$$ II_i = \frac{e_D^i}{e_U^i} = \frac{\sqrt{1 - P_{ii}}}{\sqrt{P_{ii}}}. \quad (11) $$

Low Innovation index translates into a large component of error that is not reflected in the residual. Therefore, the residual will be very small even if there is a gross error. By using (10) and (11), the CME can be expressed in terms of the residual and the innovation index:

$$ CME_i = r_i \left( \sqrt{1 + \frac{1}{II_i^2}} \right). \quad (12) $$

The CME values are thus estimated from real measurements taken from the SG, and can then be used for bad data analysis [16]. Real measurements in the context of this work are defined as measurements obtained from power systems simulators, such as ATP [11].

Regarding gross error analysis, a $\chi^2$ hypothesis testing is used for the detection of bad data in the measurement set. On such, the $J(\hat{\mathbf{x}})$ based objective function value (13) is compared to a $\chi^2$ threshold, and is the estimated state vector. The latter is based on a chosen probability $p$ (typically
$p = 0.99$) and the degrees of freedom $d_f$ of the measurement model:

$$J(\hat{x}) = \sum_{i=1}^{d_f} \left[ \frac{e_i}{\sigma_i} \right]^2 > \chi_{d_f, p}^2. \quad (13)$$

If the value of $J(\hat{x})$ is greater than the $\chi^2$ threshold, then a gross error is detected, and the signal does not follow the hypothesized distribution, i.e., modeling assumptions are not supported. Still, following [3], bad data is identified through analysis of the normalized composed measurement error, and corrected through a relaxed model [13]. This methodology is used in comparing results from this work.

## B. Modeling Methodology

PowerWorld is a power system simulation package with many analytical features. The two features of focus in this research are power flow calculations and the GIC add-on.

Power flow calculations are performed by solving the power flow equations:

$$P_k = V_k \sum_{j=1}^{N-1} [V_j[G_{kj}cos(\theta_k - \theta_j) + B_{kj}sin(\theta_k - \theta_j)]] \quad (14)$$

$$Q_k = V_k \sum_{j=1}^{N-1} [V_j[G_{kj}sin(\theta_k - \theta_j) - B_{kj}cos(\theta_k - \theta_j)]] \quad (15)$$

Because of the non-linearity in the power flow equations, they must be solved using an iterative method such as Newton-Raphson [14]. PowerWorld performs this calculation to solve power system models. Conversely, ATP is an electromagnetic transient program (EMT) which solves power system models through numerical integration of ordinary differential equations. While the input to both software is nearly identical, the internal modeling of each is unique, and therefore can be analyzed using the forward propagation approach presented previously.

PowerWorld also has a modeling feature that allows for simulation of geomagnetically induced currents. Induced currents that flow through transmission lines cause a DC bias, which can saturate transformers. This saturation, specifically half-cycle saturation, causes an increase in transformer reactive power loading that has been shown to increase linearly with GIC through the transformer [18]. Because of this approximate linear relationship, reactive power losses of the transformer can be scaled linearly with relation to the GIC flow. This is implemented in PowerWorld through an additional reactance in parallel with shunt magnetizing reactances [12]. GIC effects are further modeled in PowerWorld through the use of DC voltage sources in series with transmission lines, wherein the magnitude of voltage varies with the expected electric field based on geographical location. Such modeling features require the use of standard modeling inputs such as line resistance, transformer configuration, and parameters, as well as more obscure parameters such as substation grounding resistance, transformer coil resistance, winding characteristics, and geographical coordinates [12].

GMD effects are similarly modeled in ATP, through the inclusion of DC voltage sources in series with transmission lines near affected transformers. ATP requires similar parameters to PowerWorld such as line resistance, transformer configuration and parameters. However, inputs such as geographical coordinates are not required in ATP. Thus, GIC

characteristics are matched as closely as possible through the measured transformer neutral current in each model.

In the following, the conceptual framework towards modeling uncertainty quantification is presented.

## III. CONCEPTUAL FRAMEWORK

The goal of this work is GMD implicit modeling uncertainty quantification. For such, we demonstrate and quantify differences between synthetic data generated in ATP and PowerWorld software, considering GMD events. Utilities and the industry at large may tend towards using only one software for their modeling needs. However, in some cases, the use of only one type of software that incorporates only one kind of model, as positive sequence power flow of PowerWorld, may mean that important information about the system is omitted, especially in edge cases. In turn, this may result in misoperations or failures to operate due to unforeseen system conditions or parameters that cannot be modeled in a specific software.

The question then is, when are these omitted system information acceptable? To answer this question in this work we demonstrate and quantify the differences between synthetic data generated by ATP and PowerWorld, through a gross error analysis $\chi^2$ hypothesis testing. On such, it is assumed that $e_i$, i.e., the residual between both software generated data, considering the same modeling assumptions, has a Gaussian normal distribution with zero mean and known standard deviation equal to 1% of the measurement magnitude [4]. Considering such, a null hypothesis $H_0$ and an alternative hypothesis $H_1$ are made:

$$H_0 \rightarrow E[J(\hat{x})] = m - n \quad (16)$$

$$H_1 \rightarrow E[J(\hat{x})] > m - n \quad (17)$$

Thus, the hypothesis test is done as follows:

* If $J(\hat{x}) > C$, then reject the hypothesis $H_0$.

* If $J(\hat{x}) < C$, then reject the hypothesis $H_0$.

For the hypothesis test, $C$ is a constant to be determined and given by:

$$C = \chi_{d_f, p}^2 \quad (18)$$

The $C$ value depends then on the degrees of freedom of the measurement model ($d_f$) and the confidence level chosen ($p$). This means that if $H_0$ is true, the probability of $J(\hat{x}) > C$ is $1-p$. Thus, if the $H_0$ is rejected, the synthetic data generated by ATP and PowerWorld will have unacceptable differences, thus the modeling assumptions made are unacceptable. Modeling cases are observed at varying severity and analyzed to determine if such modeling features are able to be ignored.

In both ATP and PowerWorld models, system parameters for the IEEE 9 bus system [14], are matched as closely as possible in order to isolate modeling differences due to calculation methods rather than parameter discrepancies. Transformer grounding resistances and GIC magnitudes are varied to determine severity of error between each model for a variety of scenarios. System imbalance is further introduced in the ATP model and compared with positive sequence PowerWorld models as an additional test of modeling limitations.
<table><caption>TABLE I
CASES CONSIDERED</caption><thead><tr><th rowspan="2"></th><th colspan="4">Use-Cases</th></tr><tr><th>2Ω</th><th>10Ω</th><th>25Ω</th><th>40Ω</th></tr></thead><tbody><tr><th>0.1A</th><td>5a.1.,2.,3</td><td>6a.1.,2.,3</td><td>7a.1.,2.,3</td><td>8a.1.,2.,3</td></tr><tr><th>0.25A</th><td>5b.1.,2.,3</td><td>6b.1.,2.,3</td><td>7b.1.,2.,3</td><td>8b.1.,2.,3</td></tr><tr><th>0.5A</th><td>5c.1.,2.,3</td><td>6c.1.,2.,3</td><td>7c.1.,2.,3</td><td>8c.1.,2.,3</td></tr><tr><th>1A</th><td>5d.1.,2.,3</td><td>6d.1.,2.,3</td><td>7d.1.,2.,3</td><td>8d.1.,2.,3</td></tr><tr><th>10A</th><td>5e.1.,2.,3</td><td>6e.1.,2.,3</td><td>7e.1.,2.,3</td><td>8e.1.,2.,3</td></tr><tr><th>25A</th><td>5f.1.,2.,3</td><td>6f.1.,2.,3</td><td>7f.1.,2.,3</td><td>8f.1.,2.,3</td></tr><tr><th>50A</th><td>5g.1.,2.,3</td><td>6g.1.,2.,3</td><td>7g.1.,2.,3</td><td>8g.1.,2.,3</td></tr><tr><th>75A</th><td>5h.1.,2.,3</td><td>6h.1.,2.,3</td><td>7h.1.,2.,3</td><td>8h.1.,2.,3</td></tr><tr><th>100A</th><td>5i.1.,2.,3</td><td>6i.1.,2.,3</td><td>7i.1.,2.,3</td><td>8i.1.,2.,3</td></tr></tbody></table>

## IV. CASE STUDY

A Fast-Fourier Transform (FFT) was applied to time-series data from ATP in order to calculate phasor values of voltage for each phase at each bus within the IEEE 9-bus system. These phasor values were then converted to per unit using the system voltage and compared to the per-unit PowerWorld phasor calculations of voltage at each bus. Test cases included different phase unbalances conditions, GIC values and grounding resistances values. In the following, test results are presented.

<figure><img src="image_1.png" /><figcaption>Fig. 1. $J(\hat{x})$ Case 6i</figcaption></figure>

Table I presents the cases analyzed. Abnormally high model values corresponding to predicted unacceptable modeling assumptions are in red, and predicted acceptable assumptions are in black. Regarding Table I, dotted notation represents cases with different phase unbalance conditions of 0, 1, 2 and 3%. Second to last rows represents different GIC currents injected on transformers secondary. First row represents different transformer grounding resistances.

As an example, case 5a.1 represents injected GIC of 0.1 A, 1% load imbalance and 2 Ω of transformers grounding resistance. Considering the $\chi^2$ test, $d_f$ was equal to 100, and $p$ equal to 0.99. This translated to a $C$ value equal to 61.754.

Figure 1 illustrates the $J(\hat{x})$ for case 6i. As noted, two values are above the threshold and two are below. The two cases above the threshold correspond to predicted unacceptable modeling assumptions being upheld. For the other cases, one can see that the $J(\hat{x})$ is below the threshold value, and therefore predicted unacceptable modeling assumptions

<figure><img src="image_2.png" /><figcaption>Fig. 2. Load Imbalance Case 6i</figcaption></figure>

<table><caption>TABLE II
CONFUSION MATRIX</caption><thead><tr><th>n = 144</th><th>Actual Acceptable</th><th>Actual Unacceptable</th></tr></thead><tbody><tr><th>Predicted Acceptable</th><td>TP = 40</td><td>FP = 32</td></tr><tr><th>Predicted Unacceptable</th><td>FN = 47</td><td>TN = 25</td></tr></tbody></table>

are not being upheld. The latter represents false negatives, which are not as dangerous to GMD system planning as false positives, which represent predicted acceptable modeling assumptions not upheld.

Figure 2 illustrates the normalized residual mean for comparison 6i. As one can see from Figure 2, the values of residual increase with the load imbalance. This is of course expected, as PowerWorld does not model imbalance system operation while ATP does.

Table II presents a confusion matrix for all cases considered. As one can see, the number of true positives and true negatives, i.e., modeling assumptions that directly correspond with their predicted value, make up 45.1% of cases considered. That said, the percentage of modeling assumptions that are incorrectly supported or not, 54.9%, are considerable. Specifically, the number of false positive (FP) cases make up almost a quarter of total cases at 22.2%. Such cases are those in which modeling assumptions are expected to be acceptable, but are above the acceptable $\chi^2$ threshold, and are therefore not supported.

## V. CONCLUSIONS

This paper presents GMD modeling uncertainty quantification as well as statistical analysis of synthetic generated data. A probabilistic approach is used towards implicit modeling uncertainties quantification. Considering such, several modeling assumptions are made with different conditions simulated, including grounding resistance, voltage imbalance, GIC values and transformer saturation. Statistical analysis of synthetic data highlights that the actual modeling assumptions that correspond to their predicted conditions are less than half, 45.1%, while a significant portion of cases, 22.2%, correspond to predicted acceptable modeling assumptions that are actually not acceptable. It should be noted that the latter cases could potentially have a high impact on system operation. Further analysis is currently being made
on the specific false positive cases and uncertainty propaga-
tion potential consequences to system operation, as well as
mitigation strategies.

REFERENCES

[1] A. G. McNish, "The Magnetic Storm of March 24, 1940," *Terrestrial Magnetism and Atmospheric Electricity*, vol. 45, no. 3, pp. 359-364.

[2] NERC, "TPL-007-4 Transmission System Planned Performance for Geomagnetic Disturbance Events," https://www.nerc.com/pa/Stand/Reliability%20Standards/tpl-007-4.PDF, 2020.

[3] A. S. Bretas, N. G. Bretas, and B. E. Carvalho, "Further Contributions to Smart Grids Cyber-Physical Security as a Malicious Data Attack: Proof and Properties of the Parameter Error Spreading Out to the Measurements and a Relaxed Correction Model," *International Journal of Electrical Power & Energy Systems*, vol. 104, pp. 43-51, 2019.

[4] A. Bretas, N. Bretas, S. Braunstein, A. Rossoni, and R. Trevizan, "Multiple Gross Errors Detection, Identification and Correction in Three-Phase Distribution Systems WLS State Estimation: A Per-Phase Measurement Error Approach," *Electric Power Systems Research*, vol. 151, pp. 174-185, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0378779617302109

[5] N. G. Bretas and A. S. Bretas, "A Two Steps Procedure in State Estimation Gross Error Detection, Identification, and Correction," *International Journal of Electrical Power Energy Systems*, vol. 73, pp. 484-490, 2015. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0142061515002495

[6] 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 *2012 North American Power Symposium (NAPS)*, 2012, pp. 1-7.

[7] T. R. Hutchins and T. J. Overbye, "Power System Dynamic Performance during the Late-Time (E3) High-Altitude Electromagnetic Pulse," in *2016 Power Systems Computation Conference (PSCC)*, 2016, pp. 1-6.

[8] U. Bui, T. J. Overbye, K. Shetye, H. Zhu, and J. Weber, "Geomagnetically Induced Current Sensitivity to Assumed Substation Grounding Resistance," in *2013 North American Power Symposium (NAPS)*. IEEE, 2013, Conference Proceedings, pp. 1-6.

[9] M. Kazerooni, H. Zhu, and T. J. Overbye, "Improved Modeling of Geomagnetically Induced Currents Utilizing Derivation Techniques for Substation Grounding Resistance," *IEEE Transactions on Power Delivery*, vol. 32, no. 5, pp. 2320-2328, 2016.

[10] A. S. Bretas, K. Brown, V. Hartoonian, T. McDermott, and J. Dagle, "Towards Improved Decision Making for the Smarter Grid: Geomagnetic Disturbance Implicit Modeling Uncertainty Quantification," in *IEEE Power and Energy Society General Meeting*, 2023, pp. 1-5.

[11] H. K. Høidalen, "ATPDraw," https://www.atpdraw.net/, 2020.

[12] "PowerWorld," https://www.powerworld.com/, accessed: 2023-01-27.

[13] T. Zou, A. S. Bretas, C. Ruben, S. C. Dhulipala, and N. Bretas, "Smart Grids Cyber-Physical Security: Parameter Correction Model Against Unbalanced False Data Injection Attacks," *Electric Power Systems Research*, vol. 187, p. 106490, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0378779620302935

[14] A. S. Bretas, N. G. Bretas, J. B. London, and B. E. Carvalho, *Cyber-Physical Power Systems State Estimation*. Elsevier, 2021. [Online]. Available: https://www.sciencedirect.com/book/9780323900331/cyber-physical-power-systems-state-estimation?via=ihub=

[15] N. G. Bretas and A. S. Bretas, "The Extension of the Gauss Approach for the Solution of an Overdetermined Set of Algebraic Non Linear Equations," *IEEE Transactions on Circuits and Systems II: Express Briefs*, vol. 65, no. 9, pp. 1269-1273, 2018.

[16] N. G. Bretas, S. A. Piereti, A. S. Bretas, and A. C. P. Martins, "A geometrical view for multiple gross errors detection, identification, and correction in power system state estimation," *IEEE Transactions on Power Systems*, vol. 28, no. 3, pp. 2128-2135, 2013.

[17] N. Bretas, A. Bretas, and S. Piereti, "Innovation concept for measurement gross error detection and identification in power system state estimation," *IET Generation, Transmission Distribution*, vol. 5, pp. 603-608(5), June 2011.

[18] X. Dong, Y. Liu, and J. Kappenman, "Comparative Analysis of Exciting Current Harmonics and Reactive Power Consumption from GIC Saturated Transformers," in *IEEE PES 2001 Winter Meeting*, 2001.