The Lehtinen–Pirjola method modified for efficient modelling of geomagnetically induced currents in multiple voltage levels of a power network

The need for accurate assessment of the geomagnetic hazard to power systems is driving a requirement to model geomagnetically induced currents (GIC) in multiple voltage levels of a power network. The Lehtinen–Pirjola method for modelling GIC is widely used but was developed when the main aim was to model GIC in only the highest voltage level of a power network. Here we present a modification to the Lehtinen–Pirjola (LP) method designed to provide an efficient method for modelling GIC in multiple voltage levels. The LP method calculates the GIC flow to ground from each node. However, with a network involving multiple voltage levels, many of the nodes are ungrounded, i.e. have infinite resistance to ground, which is numerically inconvenient. The new modified Lehtinen–Pirjola (LPm) method replaces the earthing impedance matrix [Z^e] with the corresponding earthing admittance matrix [Y^e] in which the ungrounded nodes have zero admittance to ground. This is combined with the network admittance matrix [Yⁿ] to give a combined matrix ([Yⁿ] + [Y^e]), which is a sparse symmetric positive definite matrix allowing efficient techniques, such as Cholesky decomposition, to be used to provide the nodal voltages. The nodal voltages are then used to calculate the GIC in the transformer windings and the transmission lines of the power network. The LPm method with Cholesky decomposition also provides an efficient method for calculating GIC at multiple time steps. Finally, the paper shows how software for the LP method can be easily converted to the LPm method and provides examples of calculations using the LPm method.

1 Introduction

Geomagnetic disturbances produce geoelectric fields that drive geomagnetically induced currents (GIC) in power networks. These GIC flow along transmission lines and through transformer windings, where they can cause half-cycle saturation leading to harmonic generation, increased consumption of reactive power and transformer heating. These, in turn, can cause misoperation of protective relays and voltage sag and, in extreme cases, damage to transformers and system collapse (Kappenman and Albertson, 1990; Bolduc, 2002; Molinski, 2002; Kappenman, 2012; Guillon et al., 2016). A key requirement for understanding the impact of geomagnetic disturbances on power networks is the ability to model the GIC produced in a network by specified geoelectric fields. In 1985, Lehtinen and Pirjola published a landmark paper that provides the first description of a stand-alone method for modelling GIC. The Lehtinen–Pirjola (1985) method (hereafter referred to as the “LP method”) has been widely used in the geophysics and space weather community and provided the basis for GIC studies in many countries (e.g. Pirjola and Lehtinen, 1985; Mäkinen, 1993; Mäkinen et al., 1993; Thomson et al., 2005; Wik et al., 2008; Viljanen et al., 2012; Torta et al., 2014; 2017; Divett et al., 2018).

The LP method was designed at a time when mostly only the highest voltage levels of a power network were considered in GIC calculations. This was because the transmission lines at the lower voltage levels have higher resistance and so will experience smaller GIC values. However, in a desire to provide more comprehensive modelling of GIC in a power network, many modern studies are now looking to model GIC in multiple voltage levels in a power network. The LP method has been effectively used for such studies (e.g. Mäkinen, 1993; Mäkinen et al., 1993; Viljanen et al., 2012; Divett et al., 2018); however, using the LP method for multiple voltage levels involves many ungrounded nodes, thus having infinite resistance to ground, which is numerically inconvenient. Also, the main focus of the LP method was the GIC flow to ground through the transformer primary windings, which was the desired output when modelling a single voltage level of a power network. However, models for multiple voltage levels require calculation of the nodal voltages which are then used to calculate the GIC in the transformer windings (Boteler and Pirjola, 2014, 2017).

In this paper we show how the LP method can easily be modified to efficiently model GIC in multiple voltage levels of a power network by converting the LP method to calculate the nodal voltages directly. First we summarize the steps in the LP method and then show how these can be modified to give the modified Lehtinen–Pirjola method (hereafter referred to as the “LPm method”). We also show that the LPm method involves inversion of a matrix that is symmetric positive definite, allowing the use of efficient methods including sparse matrix techniques. Then we show how software for GIC calculations using the LP method can easily be converted to the LPm method and provide example calculations for the benchmark model introduced by Horton et al. (2012), including tables of values at intermediate steps, to help people transitioning their modelling from the LP method to the LPm method.

2 Lehtinen–Pirjola method

The GIC modelling method derived by Lehtinen and Pirjola (1985), the “LP” method, is produced by starting with Kirchhoff's current law that the net current flowing into a node, k, on branches from other nodes, n, is equal to the current flowing to ground from node k (LP Eq. 8),

$$\begin{array}{}\text{(1)}& i_k=\sum _{n=\mathrm{1}}^N{i}_{nk}=-\sum _{n=\mathrm{1}}^N{i}_{kn},\end{array}$$

and relates the current in a branch i_kn to the driving electromotive force (emf) e_kn (if there is one), the voltage difference between the nodes at the ends of the line v_k and v_n, and the admittance of the branch (LP Eq. 7) y_kn:

$$\begin{array}{}\text{(2)}& i_{kn}=y_{kn}[e_{kn}+(v_k-v_n\left)\right].\end{array}$$

Substituting Eq. (2) into Eq. (1) gives (LP Eq. 9)

$$\begin{array}{}\text{(3)}& i_k=-\sum _{n=\mathrm{1}}^N{y}_{kn}[e_{kn}+(v_k-v_n\left)\right].\end{array}$$

Note that, when considering multiple voltage levels, branches in the network consist of not just transmission lines, but also transformer windings. The transmission lines experience the driving emf produced by the magnetic field variations, whereas the transformer windings do not.

The driving emf in each transmission line is represented by an equivalent current source:

$$\begin{array}{}\text{(4)}& j_{kn}=e_{kn}{y}_{kn}.\end{array}$$

The equivalent current sources are then summed to give the current source directed into each node (LP Eq. 13).

$$\begin{array}{}\text{(5)}& J_k^{\mathrm{e}}=-\sum _{\begin{array}{l}n=\mathrm{1}\\ n\ne k\end{array}}^N{j}_{kn}\end{array}$$

Making this substitution in Eq. (3) gives

$$\begin{array}{}\text{(6)}& i_k=J_k^{\mathrm{e}}-\sum _{\begin{array}{l}n=\mathrm{1}\\ n\ne k\end{array}}^N(v_k-v_n)y_{kn}.\end{array}$$

Thus,

$$\begin{array}{}\text{(7)}& i_k=J_k^{\mathrm{e}}-v_k\sum _{\begin{array}{l}n=\mathrm{1}\\ n\ne k\end{array}}^N{y}_{kn}+\sum _{\begin{array}{l}n=\mathrm{1}\\ n\ne k\end{array}}^N{v}_n{y}_{kn}.\end{array}$$

The first summation represents the dependence of current i_k on voltage v_k and so gives diagonal elements of a network admittance matrix

$$\begin{array}{}\text{(8)}& Y_{kk}^{\mathrm{n}}=\sum _{\begin{array}{l}n=\mathrm{1}\\ n\ne k\end{array}}^N{y}_{nk}.\end{array}$$

The second summation represents the dependence of current i_k on all the other nodal voltages v_n and so gives the off-diagonal elements of the (symmetric) network admittance matrix

$$\begin{array}{}\text{(9)}& Y_{kn}^{\mathrm{n}}=-y_{kn},n\ne k\end{array}$$

(note that the superscript n on the left-hand side is not an index). Combining the above equations gives (LP Eq. 11)

$$\begin{array}{}\text{(10)}& i_k=J_k^{\mathrm{e}}-\sum _{n=\mathrm{1}}^N{v}_n{Y}_{kn}^{\mathrm{n}}.\end{array}$$

This can be written in matrix form:

$$\begin{array}{}\text{(11)}& \left[{\mathbf{I}}^{\mathrm{e}}\right]=\left[{\mathbf{J}}^{\mathbf{e}}\right]-\left[{\mathbf{Y}}^{\mathrm{n}}\right]\left[{\mathbf{V}}^{\mathrm{n}}\right],\end{array}$$

where the elements of column matrix [I^e] are the currents i_n, and the elements of column matrix [Vⁿ] are the voltages v_n.

LP makes the substitution

$$\begin{array}{}\text{(12)}& \left[{\mathbf{V}}^{\mathrm{n}}\right]=\left[{\mathbf{Z}}^{\mathrm{e}}\right]\left[{\mathbf{I}}^{\mathrm{e}}\right],\end{array}$$

where [Z^e] is the earthing impedance matrix. Thus,

$$\begin{array}{}\text{(13)}& v_k=\sum _{n=\mathrm{1}}^N{Z}_{kn}^{\mathrm{e}}{i}_n.\end{array}$$

Substituting Eq. (12) into Eq. (11) gives a matrix equation involving only the node to ground currents [I^e] as the unknowns:

$$\begin{array}{}\text{(14)}& \left[{\mathbf{I}}^{\mathrm{e}}\right]=\left[{\mathbf{J}}^{\mathbf{e}}\right]-\left[{\mathbf{Y}}^{\mathrm{n}}\right]\left[{\mathbf{Z}}^{\mathrm{e}}\right]\left[{\mathbf{I}}^{\mathrm{e}}\right].\end{array}$$

Gathering terms in [I^e] gives

$$\begin{array}{}\text{(15)}& \left(\right[\mathbf{1}]+[{\mathbf{Y}}^{\mathrm{n}}\left]\right[{\mathbf{Z}}^{\mathrm{e}}\left]\right)\left[{\mathbf{I}}^{\mathrm{e}}\right]=\left[{\mathbf{J}}^{\mathbf{e}}\right],\end{array}$$

where [1] is the unit matrix with size equal to the number of nodes in the model network. Equation (15) can be solved by matrix inversion to give the currents flowing to ground (LP Eq. 12):

$$\begin{array}{}\text{(16)}& \left[{\mathbf{I}}^{\mathrm{e}}\right]=\left(\right[\mathbf{1}]+[{\mathbf{Y}}^{\mathrm{n}}\left]\right[{\mathbf{Z}}^{\mathrm{e}}\left]{)}^{-\mathrm{1}}\right[{\mathbf{J}}^{\mathbf{e}}].\end{array}$$

The values of [I^e] were the desired output when modelling a single voltage level of a power network. However, if there was more than one transformer at a substation (as usually occurs), it was necessary to split the current in proportion to the admittances of the transformer windings to determine the fraction of the current that flowed in each transformer winding.

Now, when modelling the GIC in multiple voltage levels of a power network, many of the nodes are ungrounded. However, the LP method needs to specify an earthing impedance for each node. This is done by adding “virtual” connections to ground from each ungrounded node (Mäkinen, 1993; Pirjola, 2005). These virtual earthing connections have infinite resistance, but this cannot be represented in the earthing impedance matrix [Z^e], so a high value is used instead. The LP method then gives the current flow to ground from each node, including small current values through the virtual earthing connections. It is necessary to use the [I^e] values and the earthing impedance [Z^e] in Eq. (12) to calculate the nodal voltages [Vⁿ]. The nodal voltages are then used to calculate the GIC flow in the branches using Eq. (2). This is the equation to use for the GIC in the transmission lines. For branches of the network that are transformer windings, there is no driving emf, so Eq. (2) reduces,

$$\begin{array}{}\text{(17)}& i_{kn}=y_{kn}(v_k-v_n),\end{array}$$

to give the GIC flow in the transformer windings.

3 Lehtinen–Pirjola modified method

When modelling GIC in multiple voltage levels of a power network, it is necessary to calculate the nodal voltages before calculating the GIC in the transmission lines and transformer windings. In the LPm method the matrix equations are modified to provide a solution in terms of the nodal voltages. This also has the advantage that there is no need to add virtual earthing connections to ground from the ungrounded nodes.

To convert the currents flowing to ground [I^e] provided by the LP method to nodal voltages [Vⁿ], start with Eq. (15) and make the substitution from Eq. (12) (Pirjola, 2007),

$$\begin{array}{}\text{(18)}& \left[{\mathbf{I}}^{\mathrm{e}}\right]=\left[{\mathbf{Z}}^{\mathrm{e}}{]}^{-\mathrm{1}}\right[{\mathbf{V}}^{\mathrm{n}}],\end{array}$$

which gives

$$\begin{array}{}\text{(19)}& \left[{\mathbf{J}}^{\mathbf{e}}\right]=\left(\right[{\mathbf{Z}}^{\mathrm{e}}{]}^{-\mathrm{1}}+\left[{\mathbf{Y}}^{\mathrm{n}}\right]\left)\right[{\mathbf{V}}^{\mathrm{n}}].\end{array}$$

The LP method allows for the [Z^e] matrix to have off-diagonal elements representing the voltage produced at node i by currents flowing to ground from other nodes (Pirjola, 2008). However, if the circuit is set up with a node at the neutral point of each substation, this does not happen (see Boteler and Pirjola, 2014). In this case, [Z^e] becomes diagonal with elements equal to the earthing resistances r_i of the nodes, and the inverse of [Z^e] is simply the earthing admittance matrix [Y^e] given by

$$\begin{array}{}\text{(20)}& Y_{ij}^{\mathrm{e}}=\left\{\begin{array}{ll}{y}_i=\frac{\mathrm{1}}{r_i},& i=j,\\ \mathrm{0},& i\ne j.\end{array}\right.\end{array}$$

Then Eq. (19) can be rewritten as

$$\begin{array}{}\text{(21)}& \left[{\mathbf{J}}^{\mathbf{e}}\right]=\left(\right[{\mathbf{Y}}^{\mathrm{e}}]+[{\mathbf{Y}}^{\mathrm{n}}\left]\right)\left[{\mathbf{V}}^{\mathrm{n}}\right].\end{array}$$

The voltages of the nodes are then found by taking the inverse of the sum of the admittance matrices and multiplying by the nodal current sources:

$$\begin{array}{}\text{(22)}& \left[{\mathbf{V}}^{\mathrm{n}}\right]=\left(\right[{\mathbf{Y}}^{\mathrm{e}}]+[{\mathbf{Y}}^{\mathrm{n}}\left]{)}^{-\mathrm{1}}\right[{\mathbf{J}}^{\mathbf{e}}].\end{array}$$

These node voltages can then be substituted into Eqs. (2) and (17) to give the GIC in the transmission lines and the transformer windings. The GIC flow to ground is simply given from Ohm's law using the neutral point node voltage and the admittance to ground (Eq. 18 with [Z^e]${}^{-\mathrm{1}}=\left[{\mathbf{Y}}^{\mathrm{e}}\right]$).

The LPm method involves inversion of a matrix ([Y^e] + [Yⁿ]) which is symmetric and positive definite and can thus be solved using a particularly efficient case of lower–upper (LU) decomposition, the Cholesky decomposition (Press et al., 2007). Note that most of the nodes within the network are unconnected, meaning that [Yⁿ] has many zeros. This is also the case with [Y^e], so the Cholesky decomposition enables the use of sparse matrix methods (Stott and Alsaç, 1987; Press et al., 2007), thus providing an efficient way to model GIC in multiple voltage levels of a power network.

4 Calculation of GIC time series

GIC modelling is now being used, not just to assess the GIC for specified electric field values, but also to determine the variation of GIC throughout a geomagnetic disturbance. If the network configuration does not change during that time (not always the case), then the matrix inversion does not need to be recalculated at every time step.

If the electric field is assumed to be uniform across the network, then linear superposition can be used to calculate the GIC (Boteler, 2013). (A uniform electric field would be produced, e.g. if calculations are made using data from a single magnetic observatory and a one-dimensional (1-D) earth conductivity model.) The GIC modelling can be made for two cases: (i) a northward electric field of 1 V km⁻¹ and (ii) an eastward electric field of 1 V km⁻¹. For each location k in the network, this modelling gives reference GIC values α_k and β_k for the northward and eastward electric fields that can be scaled by the actual electric field values at each time step and then combined to give the time series of GIC values at that location.

$$\begin{array}{}\text{(23)}& i_k\left(t\right)={\mathit{\alpha }}_k{E}_N\left(t\right)+{\mathit{\beta }}_k{E}_E\left(t\right)\end{array}$$

This concept can be extended for using two magnetic observatories (Boteler et al., 2014), but this still requires use of a single 1-D earth conductivity model for the whole network.

In practice there is considerable variability in the earth conductivity structure across a power network. There are many modelling techniques for calculating the electric fields in such cases, ranging from use of multiple 1-D earth models (Marti et al., 2014) to use of magnetotelluric transfer functions and 3-D earth conductivity models (Weigel, 2017). In these cases, the electric fields across the network can change from place to place and from one time step to the next. This will result in a different set of nodal current sources [J^e] for each time step. However, for much of the time the network configuration may be unchanged; thus, once the inverted matrix ([Yⁿ] + [Y^e])⁻¹ has been calculated, it does not need to be recalculated at each time step and can be used with the nodal current source [J^e] for each time step to calculate the nodal voltages [Vⁿ] and hence the time series of GIC values.

However, for GIC calculations using the LPm method, even more efficient time series calculations are possible. The solution of a matrix equation such as Eq. (21) can be accomplished using LU decomposition, as explained in Press et al. (2007). This involves writing the matrix ([Yⁿ] + [Y^e]) as a product of two matrices:

$$\begin{array}{}\text{(24)}& \left[\mathbf{L}\right]\left[\mathbf{U}\right]=\left(\right[{\mathbf{Y}}^{\mathrm{n}}]+[{\mathbf{Y}}^{\mathrm{e}}\left]\right),\end{array}$$

where [L] is a lower triangular matrix and [U] is an upper triangular matrix.

For a positive-definite symmetric matrix, as is obtained with the LPm method, the [L] matrix can be chosen such that the [U] matrix is the transpose of [L]. In this case we can write (Eq. 24) using the Cholesky decomposition

$$\begin{array}{}\text{(25)}& \left[\mathbf{L}\right]\left[{\mathbf{L}}^{\mathrm{T}}\right]=\left(\right[{\mathbf{Y}}^{\mathrm{n}}]+[{\mathbf{Y}}^{\mathrm{e}}\left]\right).\end{array}$$

This Cholesky decomposition to solve the linear set is

$$\begin{array}{}\text{(26)}& \begin{array}{rl}\left(\right[{\mathbf{Y}}^{\mathrm{n}}]+[{\mathbf{Y}}^{\mathrm{e}}\left]\right)\left[{\mathbf{V}}^{\mathrm{n}}\right]& =\left(\right[\mathbf{L}\left]\right[{\mathbf{L}}^{\mathrm{T}}\left]\right)\left[{\mathbf{V}}^{\mathrm{n}}\right]\\ & =\left[\mathbf{L}\right]\left(\right[{\mathbf{L}}^{\mathrm{T}}\left]\right[{\mathbf{V}}^{\mathrm{n}}\left]\right)=\left[{\mathbf{J}}^{\mathbf{e}}\right],\end{array}\end{array}$$

by first solving for the column matrix [P] such that

$$\begin{array}{}\text{(27)}& \left[\mathbf{L}\right]\left[\mathbf{P}\right]=\left[{\mathbf{J}}^{\mathbf{e}}\right]\end{array}$$

and then solving

$$\begin{array}{}\text{(28)}& \left[{\mathbf{L}}^{\mathrm{T}}\right]\left[{\mathbf{V}}^{\mathrm{n}}\right]=\left[\mathbf{P}\right].\end{array}$$

The advantage is that the solution of a triangular set of equations is quite trivial, as Eq. (27) can be solved by forward substitution and Eq. (28) can be solved by back substitution. Also, once the Cholesky decomposition has been done, it is possible to solve with as many right-hand sides as required, one at a time. Thus the LPm method provides a much faster and versatile way of calculating a time series of GIC values.

5 Comparison between the LP and LPm methods

To illustrate the differences between the LP method and the LPm method, consider the circuit for a substation with a two-winding transformer and three autotransformers as shown in Fig. 1. The LP method requires the addition of virtual connections to ground from nodes 1 and 2, as explained above. However, in the LPm method the connection to ground is expressed as an admittance value. For the ungrounded nodes the admittance to ground is zero, which can easily be included in the earthing admittance matrix without having to add virtual connections to the circuit.

Figure 1 Substation with a two-winding transformer and three autotransformers and the equivalent circuits for the LP and LPm methods.

The steps involved in calculating GIC in multiple voltage levels of a power network using the LP method and the LPm method are summarized in Fig. 2. In the LPm method, because it involves only admittances and calculates the nodal voltages directly, there is no need to add virtual connections to ungrounded nodes, and then there is no need to convert the currents through the virtual connections to nodal voltages.

Figure 2 Comparison of the steps involved in the LP and LPm methods.

Figure 2 also shows how easy it is to convert from the LP method to the LPm method. Many steps in the process are the same. The only changes are to set up the earthing admittance matrix [Y^e] instead of the earthing impedance matrix [Z^e]. This removes the need to add virtual connections to ungrounded nodes. Then the combined admittance matrix ([Yⁿ] + [Y^e]) is formed instead of the matrix ([1] + [Yⁿ][Z^e]). After that, the matrix inversion is done and multiplied by the current source vector [J^e], the same as in the LP method (but note the comments about faster inversion methods in the previous section). An advantage of the LPm method is that the matrix calculation yields the nodal voltages directly without the need to obtain them from the earthing currents as required in the LP method. Finally, the same step is performed in the LP and LPm methods to use the nodal voltages to calculate the GIC in the transmission lines and transformer windings.

Figure 3 Single-line diagram of the benchmark test case of Horton et al. (2012).

6 Example calculation using the LPm method

To illustrate the use of the LPm technique, we present the calculation of GIC in the benchmark model of Horton et al. (2012) shown in Fig. 3. The following tables will also provide values for testing when converting software from the LP method to the LPm method.

To construct the network admittance matrix [Yⁿ] and the earthing admittance matrix [Y^e], it is first necessary to assign node numbers to the buses and neutral points at each substation. The modelling does not depend on any particular choice of node numbers. Here we use the node number assignment shown in Table 1. Note that some buses are not included in the model for calculating the GIC because they are connected to transformer windings in a delta configuration so there is no path for the GIC to flow. Also, there is no neutral point node at substation 7 because there are no transformers at this site as it is just a switching station connecting transmission lines.

Table 1 Assignment of node numbers.

Table 2 Network admittance matrix [Yⁿ] values in Siemens for the benchmark model shown in Fig. 3.

The network admittance matrix [Yⁿ] is constructed using the transmission line information and transformer information presented in Tables II and III of Horton et al. (2012). In the network admittance matrix [Yⁿ] the values are given by Eqs. (8) (for the diagonal elements) and (9) (for the off-diagonal elements). Note that the values are presented for a single phase of the power network, and it is assumed that the other two phases are identical. The network admittance matrix values for the benchmark model are given in Table 2.

The earthing admittance matrix [Y^e] is constructed using the substation grounding resistance, r_g, presented in Table I of Horton et al. (2012). Note that the GIC from all three phases flow through the substation grounding resistance, so the voltage drop here is 3 times that produced by a current from a single phase. To account for this in a single-phase model, the earthing admittance is given by

$$\begin{array}{}\text{(29)}& y^{\mathrm{e}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}{r}_{\mathrm{g}}}}.\end{array}$$

Table 3 Earthing admittance matrix [Y^e] values in Siemens for the benchmark model shown in Fig. 3.

In the earthing admittance matrix, most nodes are not connected to ground, so their earthing admittance values are zero, and there are only non-zero admittance values for the six neutral point nodes. Earthing admittance matrix values for the benchmark model are shown in Table 3.

Note that the general theory is expressed in terms of impedance and admittance, which can have reactive components, but, in practice, at the frequencies applicable to GIC the reactive components are negligible, and the network characteristics can be described as purely resistive or conductive.

Table 4 Inverted matrix ([Yⁿ] + [Y^e])⁻¹ values in Ohms for the benchmark model shown in Fig. 3.

The resulting inversion of the matrix gives ([Yⁿ] + [Y^e])⁻¹ shown in Table 4.

Horton et al. (2012) consider two cases: a northward electric field of 1 V km⁻¹ and an eastward electric field of 1 V km⁻¹. They give values for the input induced emf in each line and show the calculated output values in terms of both the nodal (bus) voltages and the GIC (A/phase) in the transmission lines and transformer windings.

Table 5 Voltages in the transmission lines and equivalent current sources for northward and eastward electric fields of 1 V km⁻¹.

The voltage source in each transmission line and the equivalent current source, calculated from Eq. (4), are shown in Table 5. These current sources are then summed (Eq. 5) to give the nodal current sources [J^e] shown in Table 6.

Table 6 Nodal current sources for northward and eastward electric fields of 1 V km⁻¹.

The nodal current sources are then combined with the inverted matrix (Eq. 22) to give the nodal voltages shown in Table 7. For nodes 1–11 these give the bus voltages shown in Table V of Horton et al. (2012). For nodes 12–18, combining the nodal voltage and the substation grounding resistance gives the GIC flow to ground for each substation shown in Table VII of Horton et al. (2012).

Table 7 Nodal voltages produced by northward and eastward electric fields of 1 V km⁻¹.

The nodal voltages substituted into Eqs. (2) and (17) give the GIC values for the transmission lines and transformer windings shown in Tables VI and VIII of Horton et al. (2012).

7 Discussion

The above example calculation shows that the LPm method provides GIC values that exactly match those for the benchmark model provided by Horton et al. (2012). The values in Tables 2 to 7 can also be used to check intermediate steps in the software used for the LPm calculations.

Any software developed to model GIC should be able to exactly match the values provided by the Horton et al. (2012) paper. The results presented in that paper are not an average of modelling results nor an approximation to the correct values but are the identical values obtained using four different software implementations. However, initial calculations involving the four different software implementations provided similar but slightly different results. Further investigation showed that the origin of the differences was in the way that distances between substations were being calculated in the different implementations. Some versions used formulas based on a spherical earth and some used formulas taking account of its non-spherical shape. It was then decided to standardize on substation latitudes and longitudes based on the WGS84 ellipsoid model of the Earth which is used by the global navigation satellite system (GNSS) for geolocation. After this, all the calculations gave exactly the same results. To get the source voltage values presented in Table 5 (and hence match the GIC results for the benchmark model) thus requires using the formulas presented in the Appendix of Horton et al. (2012) for calculating distances between substations.

Many people have used the LP method for calculating GIC in the highest voltage level of their power networks. With the increasing requirement to calculate GIC in multiple voltage levels of a power network, it is hoped that the new LPm method described above will provide an easy way to convert existing LP software. The conversion is a simple process. Just replace the earthing impedance matrix [Z^e] with the corresponding earthing admittance matrix [Y^e], form the new matrix ([Y^e] + [Yⁿ]) and do the matrix inversion. This directly gives the nodal voltages which are required to calculate the GIC in the transmission lines and transformer windings. There is no need for any “virtual” nodes or connections. Also, it is more efficient as ([Y^e] + [Yⁿ]) is symmetric positive definite and so can be solved using Cholesky decomposition, which is a special case of LU decomposition (Press et al., 2007) with U = L^T (i.e. the upper triangular factor is the transpose of the lower triangular factor).

Power networks, on average, have three transmission lines and one or two transformer windings connected to a bus, so a typical row in the admittance matrix has only five or six non-zero elements, independent of the overall network size. Thus, for larger networks, where node numbers can be in the thousands, the admittance matrix will have over 99 % of its values equal to zero. Cholesky factorization takes advantage of this fact by making use of sparse matrix methods (Stott and Alsaç, 1987; Press et al., 2007), thus additionally reducing memory usage and computation time. To examine how this affects the GIC modelling, we performed calculations for two networks using both the LP and LPm methods. The networks modelled were (1) the benchmark network of Horton et al. (2012), which has 18 nodes, and (2) the nation-wide Spanish Power Grid operated by Red Eléctrica de España (REE), which has 1388 nodes. GIC in the 400 kV part of the REE system were considered by Torta et al. (2014); for this study we include both the 400 and 220 kV levels of the REE network (see Torta et al., 2021, for reference).

Table 8 Properties of the matrices to be inverted using the LP and LPm methods for different power networks, namely the Horton et al. (2012) benchmark and REE.

Tests we did showed that the LP and LPm methods both produce matrices that are sparse, so there is potential for sparse matrix techniques to be applicable. Table 8 shows the calculation times and memory usage for GIC calculations using the LP and LPm methods. These show that memory usage was drastically reduced when using sparse matrix techniques, with the reduction being more significant with the larger REE network. The time for the matrix inversion is significantly affected, as expected, by the size of the matrix involved. For the Horton network (18 × 18 matrix) the change to sparse techniques actually increased the inversion time. However, the sparse techniques applied to the REE network (1388 × 1388 matrix) produced an approximately order of magnitude reduction in inversion time. This reduction in inversion time was greatest for the LPm method, which was nearly an order of magnitude faster than the LP method.

The “Inversion time” column in Table 8 reflects the time required to compute [Vⁿ] from Eqs. (25), (27) and (28), thus including the Cholesky decomposition and the forward and back substitutions in the LPm method; note, in consequence, that it is not strictly an inversion, though we will refer to it as such. Also note that, when referring to LP, this column reflects the time to compute [I^e], including the decomposition of $\mathbf{M}=\left[\mathbf{1}\right]+\left[{\mathbf{Y}}^{\mathrm{n}}\right]\left[{\mathbf{Z}}^{\mathrm{e}}\right]$. However, M is not even symmetric in the LP method, so the decomposition of M is indeed an upper–lower (UL) factorization. The difference in speed of the inversion between the LP and LPm methods is that LPm involves inversion of a symmetric positive-definite matrix which allows the use of the technique of Cholesky decomposition, that significantly reduces the time of the inversion process (note that LU requires the determination of more unknowns). The parameters of the calculations presented in Table 8 are obtained from GIC modelling using programs in MATLAB. Specific inversion times and memory usage will vary with the programming language used, but it is expected that the general results presented here will apply regardless of the programming language used.

8 Conclusions

We have presented a new version of the LP method, modified for efficient modelling of GIC in multiple voltage levels of a power system. In the LPm method the earthing impedance matrix, [Z^e], is replaced with an earthing admittance matrix [Y^e] that is added to the LP network admittance matrix [Yⁿ] to give a new combined matrix ([Yⁿ] + [Y^e]) to be inverted.

Multiplication of the inverted matrix ([Yⁿ] + [Y^e])⁻¹ (or equivalently recycling its Cholesky decomposition) by the nodal current sources [J^e] provides a direct calculation of the nodal voltages [Vⁿ]. These nodal voltages are then used to calculate the GIC in transmission lines and transformer windings.

Guidance is provided for converting software from the LP method to the LPm method and an example calculation using the benchmark model of Horton et al. (2012) is presented to provide a set of values for testing GIC calculation software.

Calculations of GIC using the LPm method involve a matrix that is symmetric positive definite. This enables a solution to be obtained by Cholesky decomposition, (a specific case of LU decomposition), which is numerically more accurate than computing the matrix inversion itself. The factorization of Cholesky decomposition can always be implemented using sparse matrix techniques, speeding up the calculations for large networks.

Thus the LPm method provides an efficient method for calculating GIC in multiple voltage levels in a power network that provides a valuable tool for assessing the geomagnetic hazard to power systems.