Optimal Blocking Devices Placement for Geomagnetic Disturbance Mitigation Based on Sensitivity of Induced Geoelectric Fields

When geomagnetic storms occur, geomagnetically induced currents (GICs), as a kind of quasi-DC, flow through the transformer and cause the increase of reactive power loss (<inline-formula> <tex-math notation="LaTeX">$Q_{\mathrm {GIC}}$ </tex-math></inline-formula>). Several mitigation actions (e.g., changing grid topology or installing blocking devices) exist that can reduce the harmful GIC effects on the grids. Making such decisions can be challenging, because the magnitude and direction of the induced geoelectric fields (IGFs) are uncertain and non-stationary. In this paper, the sensitivity of voltage and GIC to IGF is first calculated, and the joint voltage/GIC sensitivity feature space is constructed based on principal component analysis. On this basis, scenario clustering is conducted to obtain voltage/GIC sensitivity scenarios. This scenario not only reflects the uncertainty of IGF itself, but also reflects the impact of geomagnetic disturbance (GMD) on the system operation state. Combined with the structural parameters of Jiangsu power grid, the placement of blocking devices (BDs) is optimized. Compared with the constant IGF scenario and the IGF intensity scenario, the expected value of <inline-formula> <tex-math notation="LaTeX">$Q_{\mathrm {GIC}}$ </tex-math></inline-formula> based on the sensitivity scenario decreases more.


I. INTRODUCTION
When geomagnetic disturbances (GMDs) occur, the neutral point potentials of transformers in different geographical locations are different, and the geomagnetically induced currents (GICs) flow in the circuit composed of transmission lines, grounding transformers and the earth, resulting in the increase of reactive power loss (Q GIC ) of transformers [1]. Although the induced geoelectric field (IGF) intensity in the middle and low latitudes is weaker than that of in high latitude, with the application of large cross-section and multi-split conductors and the development of ultra-high power grid [2], the reduction of line resistance will lead to the increase of GICs, which will lead to the increase of Q GIC .
The associate editor coordinating the review of this manuscript and approving it for publication was Gab-Su Seo .
It is of practical significance to study the impact of Q GIC with large amount and strong volatility on voltage stability [3], [4].
Barnes P.R. pointed out that the reason for the Quebec blackout was that the reactive power loss of the transformer caused by GIC increases, and the overload of the capacitor bank is forced to withdraw, resulting in unbalanced reactive power and voltage collapse [5]. During the geomagnetic storm on May 10, 1992, the monitoring data of the 345kV/115kV transformer of a substation by the Electric Power Research Institute (EPRI) company showed that the Q GIC was very large [6]. Even though GMD events are rare in comparison with other more common uncertain events, it has been identified as a ''high-impact, low-frequency'' event causing risk to the power systems as stated in the report from the U.S. Department of Energy [7]. Geomagnetic storms occur simultaneously on a global scale, so the Q GIC has the characteristics of mass generation and suddenness in the whole network. There are many transformers in the power grid, and the total Q GIC is very large, which can easily cause the voltage drop of the power grid, so that the problem of voltage stability becomes prominent [8]. As a natural phenomenon, GMDs will continue to occur. With the continuous expansion of the scale of the power grid, especially the large-scale construction of ultra-high voltage (UHV) power grids and the widespread application of singlephase autotransformers, these factors will increase the possibility of the power grid being infringed by GIC, resulting in a wider range of Q GIC fluctuations and a deeper impact.
Although the geomagnetic field is usually recorded in eastward B x and northward B y directions as a function of time, the actual magnitude and direction of geomagnetic field are uncertain. Moreover, the correlation between B x and B y is unknown [9]. Several mitigation actions can reduce the harmful GIC/Q GIC effects on the electric transmission grids. First, installing direct current (DC) blocking devices at substations can prevent GICs. Second, changing the network topology by deactivating a few transmission components (e.g., transmission lines, transformer, and generators) can prevent damage due to GICs. Both methods achieve the governance effect by changing the flow path of the GIC. This study is aimed at the installation of capacitive BDs devices. In [10], the objective function is to minimize the installation cost of BDs, and the node voltage and generator reactive power limit are used as constraints to optimize BDs placement. In [11], the amplitude of the IGF is set as 1 V/km, and the maximum value of the system Q GIC under different IGF directions is calculated in combination with the GIC benchmark model, which is used to optimize BDs placement. The IGF introduced by GMDs has been assumed to be a constant to simplify GIC ang BDs analysis in previous studies. However, both the amplitude and direction of GMDs usually keep changing throughout time. As a result, it is more realistic and accurate to evaluate GMDs damage and solve optimal BDs placement problems considering the uncertainty of IGFs. In [12], the GMD event from March 13 to 14, 1989 was used as the ''benchmark'' event to establish a time-varying IGF model, and the duration of a certain IGF value was used as a weight to optimize BDs placement. In [13], the distribution law of IGFs and GIC of multiple GMDs events has been studied, but the optimal treatment research based on multiple GMDs has not been carried out. In [14], the uncertain IGF is modelled using the distributionally robust optimization approach that worst-case expectation of the system cost is minimized.
In this paper, the IGFs is calculated by using the geomagnetic field data recorded by the geomagnetic observatory in multiple geomagnetic storm events, and the frequency statistics is carried out to establish the probability distribution model of the IGF. The voltage/GIC sensitivity scenario model is introduced, and the scenario clustering results in the sensitivity space are used to guide the selection of the placement of BDs, which not only reflects the randomness of the IGF itself, but also reflects its influence on the system operation. The optimal BDs placement is carried out based on constant IGF scenario, time-varying IGF scenario and sensitivity scenario, and then regenerate IGFs samples to verify the optimization effect. It is found that the Q GIC based on sensitivity scenario is smaller than other scenarios.

II. ANALYSIS OF IGF BASED ON MULTIPLE GEOMAGNETIC STORM EVENTS A. PROBABILITY DISTRIBUTION CHARACTERISTICS OF IGF
The layered model of earth resistivity in Jiangsu area and plane wave method [16] as shown in Figure 1 are used to calculate IGF.
The uncertainty of IGF during geomagnetic storms makes Q GIC fluctuate. Determining the probability distribution function (PDF) of IGF is the basis for studying the impact of multiple geomagnetic storms on system stability. For power grids in middle and low latitudes, it is generally believed that only strong GMDs will affect their stable operation. As a standard for evaluating GMD intensity, the geomagnetic storms with D st index less than −100nT are strong geomagnetic storms. In this paper, the geomagnetic storm event [13] with D st index less than or equal to 100 nT in the 23rd solar activity cycle is taken as typical cases to evaluate the impact of GMDs on power system. The IGFs are calculated by using the layered earth conductivity model and the plane wave method [15], and its frequency statistics are carried out according to the interval, as shown in Figure 2.
To quantitatively compare the fitting effects of several common distribution functions, the fitting index is defined as shown in (1).
where M is the number of groups of frequency distribution histogram; y i = f C i , N i and C i are the height and center position of the i-th straight square column respectively; f (·) VOLUME 10, 2022  is the fitted probability density function; y i is the value corresponding to the fitting probability density function at the center position C i . The smaller the fitting index I , the higher the fitting accuracy. Using different distribution functions to fit the IGF component of 29 geomagnetic storm events, the fitting index values are obtained. The average, maximum and minimum values are shown in Table 1. Compared with logistic distribution, normal distribution and Cauchy distribution, no matter E x or E y , the fitting index value of t location-scale distribution is the smallest, so the fitting effect is the best. The probability density distribution function of the t location-scale distribution is shown in (2).
where (·) is the gamma function. µ is the position parameter. σ is the scale parameter. υ is the shape parameter.

B. IGF SAMPLING
Since the PDF of the IGF is known, and to avoid repeated sampling, Latin hypercube sampling (LHS) is used. LHS consists of two steps: sampling and sorting. Sampling can make the values obtained from the model cover all the variables evenly, and sorting can solve the correlation between variables.
For E x and E y , the probability distribution is divided into N intervals, and the corresponding variable values are obtained by inverse transformation of these N intervals. The random sampling process is shown in Figure 3.
This LHS sample covers all the information. To break the correlation of random variables, sorting is required. The generated 1∼N randomly arranged matrix L, calculate the correlation coefficient of each row of elements in L, and then form the correlation coefficient matrix ρ. The unit matrix G = D −1 L is obtained by Cholesky decomposition ρ = DD T , and the elements of each row in are rearranged to obtain independent samples.

III. OPTIMAL BDS PLACEMENT BASED ON SENSITIVITY SCENARIO
A. VOLTAGE SENSITIVITY For the system with N B nodes, the power flow equation expressed in polar coordinates can be linearized, as shown in (3).
where elements J Pθ , J PV , J Qθ and J QV of augmented Jacobian where S VQ is voltage/reactive power sensitivity matrix, S VQ = . S VQ reflects the linearized incremental relationship between system node voltage and reactive power, and each column of the matrix reflects the influence of each node on system voltage. Therefore, the elements in each column are added and processed to obtain a set of 1×(N B -1) dimension vector, called sensitivity vector, which is used to measure the impact of reactive power injection at each node on system voltage, as shown in (5).
According to the GIC equivalent model of line and transformer, the N -node network as shown in Figure 4 is established. The equivalent resistance between any two nodes j and k is R jk , and the equivalent voltage source is shown in (6).
where E is the IGF vector. dl is the micro element of the line length. E x and E y are values of northward and westward IGF respectively. l x and l y are the northward and westward length of the line respectively. Define grounding resistance matrix Z of N ×N , where N is the number of nodes, then the relationship between GIC I = [I 1 , I 2 , · · · , I N ] T and the voltage U = [U 1 , U 2 , · · · , U N ] T is shown in (7).
where Z is the diagonal matrix, and its element is the equivalent resistance of each grounding branch. According to the Kirchhoff theorem and the relationship between circuit variables in Fig.1, equation (8) and (9) can be obtained.
Based on the method proposed by Lehtinen and Pirjola [16], the N × N network admittance matrix is introduced, and its element is shown in (11) and (12).
Define N × M matrix α and β, and its element is shown in (16) and (17).
Combine two N × M matrix α and β into a N × 2M matrix A p = [α, β]. Combine two M × 1 matrix E x and E y into a 2M × 1 matrix E = E x , E y T . The equation (18) can be described as (19), and (20) is obtained by substituting (19) into (14).
The sensitivity of the GIC flowing through the node to the IGFs is shown in (21).

C. OPTIMAL BDs PLACEMENT BASED ON SENSITIVITY SCENARIO
The IGF samples are substituted into the system, and the voltage sensitivity S 1 and GIC sensitivity S 2 are calculated. The first s principal components whose contribution rate is greater than a certain set value are obtained by principal component analysis (PCA) [17], which constitutes the joint sensitivity feature space of voltage and GIC. Sensitivity scenarios related to system voltage stability are then generated through cluster analysis. The obtained scenarios are mapped to IGFs clustering to obtain IGFs sensitivity scenarios. Combined with K-means clustering method and joint sensitivity information, a sensitivity scenario model is established. The process of optimizing BDs placement based on sensitivity scenarios is shown in Figure 5. The specific process is as follows: 1) A probability distribution model of the induced geoelectric field is established. Taking the GMDs in the 23rd solar cycle as typical events, the IGFs are calculated using the geomagnetic field data recorded by the Beijing Geomagnetic Observatory, combined with the layered ground resistance model and the plane wave method.
2) Get an IGF samples. The IGF samples are obtained by Latin hypercube sampling, the GIC is calculated in combination with the network structure parameters, and the Q GIC is calculated by the K-value method [18].
3) Sensitivity calculation. The joint sensitivity is calculated using the methods in Sections III-A and B.
4) Scenario clustering. The sensitivity feature space is formed by PCA method, the optimal number of scenarios is determined by the clustering index K DBI [19], and m joint sensitivity scenarios are obtained by K-mean clustering method. 5) Optimize BDs placement. The BDs optimization problem based on the joint sensitivity scenario can be described as: under the constraint that the GIC of all substations is less than the limit value (18A is the requirement of DL/T 621-1997 AC Electrical Installations Grounding), determine the BDs placement and minimize the number of BDs. A binary-coded genetic algorithm is used to optimize the placement of BDs. The 0 element on a chromosome indicates that BDs are not installed, and 1 indicates that BDs are installed. The objective function is to minimize Q GIC , as shown in (22). Constraints include power balance, active/reactive power output, voltage limit, and node flow through GIC maximum limit.
where T is the number of scenarios. p j is the proportion of

IV. EXAMPLE ANALYSIS A. INFLUENCE OF GMDs ON VOLTAGE STABILITY
Taking the Jiangsu power grid as shown in Figure 6 as an example, the influence of GMDs on voltage stability and the optimization BDs placement based on sensitivity scenarios   are studied. Compared with other regions, the transmission lines of Jiangsu power grid are complex, and the coupling degree between power plants and stations of different voltage levels is relatively high, including 6 1000 kV substations, 110 500 kV substations, and 635 220 kV substations. The transformer parameters of each voltage level are shown in Table 2. GIC and Q GIC are calculated based on the structural parameters of Jiangsu power grid. Some substations with larger Q GIC are shown in Table 3.  70.52% and 63.55%. It can be seen from the calculation data that the Q GIC of 500kV substation is much larger than that of 1000kV and 220kV substation, and the proportion is more than 50%. Although the resistance of 500kV transmission line and the equivalent resistance of transformer are less than that of 220kV network, the number of lines and substations are far more than that of 220kV network. Although the equivalent resistance of 1000kV network is less than 500kV network, the number of 1000kV transmission lines and transformers is far less than 500kV, so 500kV transformers produce the most Q GIC . The voltage of each node when GMDs occur or does not occur is shown in Table 4.
The drop of node voltage of 500 kV substation is significantly greater than that of 1000 kV and 220 kV substation, which is mainly due to the large Q GIC of 500 kV transformers, resulting in serious reactive power shortage.

B. OPTIMAL BDs PLACEMENT BASED ON SENSITIVITY SCENARIO
The sensitivity of GIC flowing through the neutral point of transformer to IGFs is shown in Table 5.
It can be seen from Table 4 that the sensitivity of the GIC flowing through the 500 kV transformer to IGF is significantly greater than that of the GIC flowing through the 1000 kV and 220 kV transformer. Therefore, GMDs have a greater impact on 500 kV substations. In addition, for the 220 kV and 500 kV substations, the sensitivity of the Yunhe and Xinfeng substations at the ''corner'' position is significantly greater than that of other substations. Therefore, special attention needs to be paid to the influence of GMD on the substation at the end of the line.
Using the PCA method, the results of the principal component analysis can be obtained. As shown in Table 6, it can be seen that the contribution rate of the first principal component reaches 35.87%, and the cumulative contribution rate of the first five principal components reaches 98.75%. Therefore, it can be considered that the sensitivity scenario clustering space composed of the first seven principal components can be used for cluster analysis.
IGF is calculated by multiple GMD events combined with earth resistivity model, and then IGF intensity scenario is obtained by cluster analysis. Calculate the voltage sensitivity and GIC sensitivity to IGF. Then the joint sensitivity space is constructed from voltage sensitivity and GIC sensitivity, and the sensitivity scenario is finally obtained through principal component analysis and cluster analysis. The number of scenarios is determined according to David-Bouldin index K DBI . The relationship between K DBI and cluster number is shown VOLUME 10, 2022  in Figure 7. For the IGF intensity scenario and the sensitivity scenario, when the number of clusters is 13 and 9 respectively, K DBI gets the minimum value. The cluster centers and their proportions of different scenarios are shown in Table 7.
The IGF intensity scenarios are classified only according to the amplitude of E x and E y . The sensitivity scenario disturbs the distribution law of the IGFs according to the intensity. Compared with the IGF intensity scenario, the sensitivity scenario highlights the impact on voltage stability, which not only reflects the uncertainty of IGFs, but also reflects the characteristics of system operation.
In the genetic algorithm, the initial population size is 100, the crossover probability is 0.9, the mutation probability is 0.2, and the maximum evolution algebra is 100. The minimum fitness is defined as the minimum number of BDs installed in all populations in each evolution when the maximum GIC of all substations is less than 18 A, and the average fitness is defined as the average number of BDs installed in all populations in each evolution. The placement of BDs is optimized based on constant IGF scenario (1V/km), IGF intensity scenario and sensitivity scenario. Constant IGF scenario is to assume that IGF is a constant value. IGF intensity scenarios only consider the randomness of IGF, and conduct scenario analysis for IGF. Sensitivity scenarios cluster analysis is conducted for voltage sensitivity and IGF sensitivity, and IGF scenario division is based on the influence of GMD on system operation characteristics. For the constant IGF scenario (E x = 1 V/km), according to the results of optimized installation placement, BDS is installed in Anji substation, Dongming substation, Maoshan substation, Jinling substation and Mengxi substation, respectively. For the IGF intensity scenario, BDs are installed at Anji substation, Tianmuhu substation, Fangxian substation, and Xijindu substation, respectively. For sensitivity scenarios, BDs are installed at nodes Yangxian substation, Fangxian substation, Tianmuhu substation and Minzhu substation, respectively. To verify the governance effect of BDs on uncertain GMDs, 100 sets of IGF samples are regenerated, and the changes of total Q GIC before and after the installation of BDs are calculated as shown in Table 8.
Optimizing the BDs placement based on the three scenarios reduces the Q GIC and improves the voltage stability. However, for IGFs with different amplitudes and directions, the average Q GIC after optimizing BDs placement based on sensitivity scenario is the smallest after installation. Therefore, the sensitivity scenario is more suitable for solving the problem of optimizing the BDs placement.

V. CONCLUSION
In view of the uncertainty of IGF when geomagnetic storms occur, taking multiple geomagnetic storms in the 23rd solar activity cycle as a typical case, three parameter t distribution is used to fit the IGF distribution function, and IGF samples are obtained by Latin hypercube sampling method. A sensitivity scenario probability model is proposed, and the placement of BDs is optimized based on this model. Through simulation analysis, the following conclusions are drawn: From the model itself, the IGF sensitivity scenario can not only reflect the uncertainty of the IGF, but also reflect the influence of the IGF on the system voltage stability. From the actual effect, compared with the constant IGF scenario and the varying IGF intensity scenario, after optimizing the placement of BDs based on the sensitivity scenario, the Q GIC is reduced more and the risk of voltage instability is reduced more. Therefore, the optimized BDs placement based on sensitivity scenarios is more advantageous.
Compared with the constant IGF scenario and IGF intensity scenario, because the sensitivity scenario comprehensively considers the uncertainty of IGF and the impact of GMD on system stability, the optimal configuration of BDs based on the sensitivity scenario not only avoids the problem of too large IGF and insufficient governance, but also avoids the problem of small IGF and waste of governance devices. For IGF with different amplitude and direction, good governance effect is achieved.