A multi-objective reactive power optimization strategy for mitigating voltage fluctuation in power network caused by geomagnetic storm

Abstract

Geomagnetically induced currents through the transformer windings produce biases of transformers, increase transformer reactive power losses, and lead to a significant lack of reactive power support in the power system. As a result, the geomagnetically induced currents cause voltage fluctuations and affect the safe and stable operation of the power network. This study focuses on the problem of multi-objective reactive power optimization strategy to mitigate the negative effects of geomagnetically induced current on the power network. Using the stochastic response surface method, this article establishes a model of the probabilistic power flow calculation, combines a particle swarm optimization algorithm based on the niche chaotic mutation, and sets up a multi-objective reactive power optimization model. The introducing of chaos mutation avoids the result falling into a local optimal solution and improves the global search ability. The geomagnetically induced current–benchmark simulation verified the accuracy and effectiveness of the model and algorithm, and proved that the strategy can effectively reduce the voltage deviation caused by the geomagnetically induced current.

Keywords

Geomagnetic storm geomagnetically induced currents voltage fluctuation reactive power optimization stochastic response surface particle swarm optimization

Introduction

Solar storms caused by violent outburst of explosive activities on the sun may lead to geomagnetic disturbances (GMDs). GMDs cause geomagnetically induced currents (GICs) flowing through long-distance transmission lines, transformer neutral groundings, and the earth.¹ GICs cause transformer magnetic flux saturation, lead to transformer heating, increase in transformer reactive power demand, and harmonics. A severe GIC may lead to transformer damage or a power grid blackout.² After the Quebec blackout caused by a GMD in 1989, PR Barnes pointed out that the blackout was caused by two factors: (1) GIC caused the increase in transformer reactive power demand and (2) the existence of the harmonics forced the capacitor banks to drop out, which caused a voltage collapse triggered by the reactive power imbalance.³ With the monitoring data at a 345/115-kV transformer during a GMD on 10 May 1992, Electric Power Research Institute (EPRI) in the United States showed that the transformer reactive power demand caused by the GIC (GIC-Q for short) is significant.⁴ GMDs occur approximately at the same time in the global context; therefore, the GIC-Q has characteristics such as cut-in clusters and abruptness. Because the number of power transformers is large, the total amount of the GIC-Q will be large enough to cause power system voltage instability issues or large-scale voltage collapse.⁵

After the Quebec blackout, many countries conducted research to determine the numerical relationship of the GIC-Q. The mechanism of transformer GIC-Q was revealed comparatively early by R Pirjola.⁶ J Berge⁷ described the linear relationship of the GIC-Q by means of experimental simulation and field data. Marti et al.⁸ deduced the linear relationship between the GIC-Q through effective susceptance. Liu et al.⁹ calculated the numerical relations of the GIC-Q using the V – I curves provided by the transformer manufacturers and the $ψ - I$ curve of the iron core material, and derived the linear relations of the GIC-Q. Those research results revealed the demand mechanism and the number relations of GIC-Q, and explained the distribution of power grid reactive power demand under a GMD, which provides the basis for power grid reactive power optimization control, voltage stability analysis, and risk assessment of voltage instability.

In recent years, with the development of power system, many GMD issues were identified in the mid-low latitude area including South Africa, Brazil, Japan, China, and some other places.^10–14 GMDs are natural phenomenon and happen from time to time. The continuous expansion of the power grid in China, especially the large-scale construction of ultra-high-voltage (UHV) power grid and the use of single-phase autotransformers, will increase the probability of GIC events, and thereby cause the grid GIC-Q having a wider effective range and a stronger influence.^15,16 Therefore, the power company must not only take reasonable defensive measures against GMD but also study the mitigating schemes for GIC-Q. Neutral point blocking devices can effectively reduce the reactive power demand caused by the GIC,^17,18 but their impacts on power system stability still need further research. Making full use of the existing reactive power compensation equipment, reasonable allocation of the blocking devices and optimization of the scheduling are all feasible defensive safeguards against the GMD.

However, a precondition is to forecast accurately the solar storms for defending GIC-Q through the effective scheduling. For this reason, since 1989, when the solar storm caused the Quebec power grid blackout, the solar storm forecasting has gained great importance.¹⁹ The Advanced Composition Explorer (ACE) satellite, launched in August 1997, is designed to monitor the solar storms at the sun–earth gravitational equilibrium point. It can predict the occurrence of the magnetic storms 15–60 mins in advance, and promptly remind the power company to attract attention, so that it provides a protection for the power grid from the space weather disaster.²⁰

For the target of the minimum expectations of the offset voltage and reactive power compensation capacity, by considering the uncertainty of the GIC-Q, this article proposes a reactive power optimization strategy, which is the combination of the stochastic response surface method (SRSM) and multi-objective particle swarm optimization (MOPSO) algorithm based on niche chaotic mutation. Simulation results with GIC standard examples show the practicability and validity of the method.

The rest of this article is organized as follows. After a general description of the calculation method of the GIC and GIC-Q in section “GIC-Q modeling,” and the establishment of a probabilistic power flow calculation model for the case of considering the effect of GIC-Q in section “Probability trend model based on SRSM,” the multi-objective reactive power optimization algorithms are presented in section “Multi-objective optimization algorithm,” and a standard numerical simulation analysis of GIC is presented in section “Example analysis.” Finally, section “Conclusion” summarizes the article.

GIC-Q modeling

Primarily, there are two risks associated with GIC in the power system. The first is the transformer’s permanent damage, primarily the high-voltage (HV) transformers. The second is the demand of reactive power support leading to the potential frisk of a voltage collapse. This article focuses on the second risk. The main problem of reactive power demand caused by the GIC represents a significant threat to the power grid. Considering the GIC-Q is of great significance for the study of reactive power optimization.

GIC modeling

Since the current system of magnetosphere, ionosphere, is complex, using deep magnetotelluric survey data, establishing equivalence earth resistivity model, and calculating inductively field are common methods. Simultaneously, the feature of the ionospheric current direction is not obvious in the mid-low latitude. The electric field may occur at any direction. For this reason, without accurate earth electrical and structural data, we can select the electric field of 1 V/km to compute GIC. So, the transmission line equivalent to DC voltage U_ij is only related to the location of line ends

U_{i j} = E_{N} L_{N} + E_{E} L_{E}

(1)

where E_N/E_E is the North/East component of the electric field, respectively (V/km), and L_N/L_E is the North/East distance of the line between the two end points, respectively (km).

The equivalent node is obtained by the line-induced DC voltage U_ij, and the node is injected with a DC current source. Furthermore, with the DC network parameters of the components including transmission lines, transformer windings, and the ground resistance of the substations, the GIC in each branch of the network can be calculated from the network node voltage. The detailed calculation can be seen in Zheng et al.²¹

Equivalent GIC calculation of a transformer

The path of GIC is represented as a DC network which consisted of the transmission lines, transformer windings, and the ground grid of the substations. Thus, a DC model of the transformer is needed to study for GIC calculation. The corresponding transformers can be generally categorized as autotransformers and step up/down transformer which are divided into two types, YN0-YN0 and YN0-D; their DC models are illustrated in Figure 1.

Figure 1.

DC model of an (a) autotransformer, (b) step up/down (YN0-YN0), and (c) step up/down (YN0-D) transformers.

In the case of autotransformers shown in Figure 1(a), the HV and low-voltage (LV) terminals of the transformer are electrically connected through the series winding, featuring the DC resistance R_S. The path to the ground is provided through the common winding of the autotransformer, with the DC resistance denoted by R_C. The two windings are the flow paths of the GIC.

During the period of the GMD, the flow path of the GIC is different between the series winding and common winding of an autotransformer. The equivalent GIC is

I_{S} N_{S} + I_{C} N_{C} = I_{e f f} (N_{S} + N_{C})

(2)

Namely

I_{e f f} = \frac{[(N - 1) I_{S +} I_{C}]}{N}

(3)

where I_eff is the equivalent GIC, I_S is the GIC flowing into the series winding, I_C is the GIC flowing into the common winding, N = V_H/V_L, and V_H and V_L are the terminal voltages of HV side and LV side, respectively.

In the case of a YN0-YN0 transformer shown in Figure 1(b), the GICs will flow in two windings. The same equivalent GIC flowing into the transformer is

I_{e f f} = \frac{I_{H} + I_{L}}{N}

(4)

where I_H and I_L are the GIC flowing into HV winding and LV winding through the DC resistance R_HV and R_LV, respectively, and N is the same as in equation (3).

In the case of a YN0-D transformer shown in Figure 1(c), the HV and LV windings are separate, and only the HV winding is modeled. The reason is that the LV winding can be ungrounded. In this case, the effective GIC of the transformer winding is one-third of the GIC through the transformer’s neutral point.

GIC-Q modeling

The reactive power demand of a transformer caused by GIC is calculated based on the method in Dong et al.⁷ and Rezaei-Zare²² which is

Q_{G I C} = K U_{P U} I_{G I C}

(5)

where K is the GIC’s reactive loss coefficient obtained from experimental testing with the unit of MVar/A; the coefficients K = [1.8, 1.18, 0.33, 0.29, 0.66] correspond to a single-phase autotransformer, single-phase-transformer group, three-phase shell-type transformer, three-phase three-column core transformer, and three-phase five-column transformer, respectively; and U_PU is the per-unit voltage at the transformer’s terminal.

Voltage fluctuation theory caused by geomagnetic storm

Figure 2 shows the two-node system which forms a DC network model to calculate GIC, where V_i and V_j are, respectively, both ends of the transmission line voltage; U_dc is the DC line-induced voltage caused by the GMD; R_dc is the DC resistance of the line; R_S/R_C are, respectively, autotransformer series winding and common winding DC resistance; and R_g is the transformer’s neutral grounding DC resistance.

Figure 2.

Diagram of the two-node network.

Taken ${\overset{\cdot}{V}}_{j}$ as a reference axis, the line voltage drop is expressed as follows

Δ V = \frac{(P_{j} + Δ P) R_{dc} + (Q_{j} + Q_{G I C} - Q_{C j}) X}{V_{j}}

(6)

where node j are, respectively, injected active and reactive power by P_j and Q_j; Q_Cj is the charging reactive power of the line; X is the line reactance; and $μ_{GICi}$ P is the power loss caused by transformer’s reactive power loss due to GIC, so

Δ P = \frac{(2 Q_{j} Q_{G I C} + Q_{G I C}^{2}) R_{dc}}{V_{j}^{2}}

(7)

Therefore, under the condition that the transmission power and parameters of the line are unchanged, the voltage drop during storms is only related to Q_GIC.

Probability trend model based on SRSM

GIC is quasi-DC with characteristics of deep penetration into the earth, however, to obtain the earth’s accurate conductivity over the large scope and depth of 400–500 km is quite difficult. Different land areas have different values of electrical conductivity. Currently, there is still a lack of accurate ground conductivity model, so that it is difficult to accurately calculate the induced electric field caused by the geomagnetic storm and GIC in the power grid. In this situation, the calculation results of the reactive power losses of transformers caused by the GIC must have some uncertainty. There is a deviation in the voltage distribution calculation which makes it not a good guide to the actual power network operation. Therefore, it is necessary to use the probability method to calculate the voltage distribution caused by the magnetic storm, so that the calculation results can be better to guide the practice. In this article, the SRSM is used to model the probability flow.

The SRSM is a typical method to analyze the uncertainty, whose basic principle is that the probability distribution of an n-dimensional random variable is mapped to a Hermite chaos polynomial of known parameters. By selecting a few typical samples to determine the undetermined coefficients of the polynomial, the probability distribution of the output response is achieved.

In terms of the uncertainty of reactive power demand of the transformer caused by a GMD, there is still no detailed model of the probability distribution. Consequently, this article uses a normal distribution probability model to describe the GIC-Q uncertainty existing in the calculation process. For node i, the GIC-Q is called Q_GICi, and the probability density function is

f (Q_{G I C i}) = \frac{1}{\sqrt{2 π} σ_{G I C i}} \exp (- \frac{{(Q_{G I C i} - μ_{G I C i})}^{2}}{2 σ_{G I C i}^{2}})

(8)

where $μ_{GICi}$ and $σ_{GICi}$ are the mean and standard deviation of Q_GICi, respectively, taking the calculation value as the mean and 10% as standard deviation.

The GIC-Q is noted as the uncertain parameter X, and the solutions of the power flow calculation are recorded as output response Y of the system. N random input variables $X = [x_{1}, x_{2}, \dots, x_{n}]$ and output responses Y have a mapping relationship as follows

Y = F (X) = F (x_{1}, x_{2}, … x_{n})

(9)

where F(X) is the distribution function of x.

The requirements of the SRSM are not only that the random input variables are independent of each other but also that they must obey the same standard distribution. Hence, this study switches the input random variables to random variables obeying the standard normal distribution.

The different distributions of the random variables are represented as

x_{i} = F_{i}^{- 1} (Φ (ξ_{i}))

(10)

where $Φ (ξ_{i})$ is the distribution function of $ξ_{i}$ ; $ξ_{i}$ is a random variable obeying the standard normal distribution; and $F^{- 1} (\cdot)$ is the inverse function of $F (\cdot)$ . The output response Y is expressed as Hermite chaos polynomial

\begin{array}{l} Y = a_{0} + \sum_{i_{1} = 1}^{n} a_{i_{1}} H_{1} (ξ_{i_{1}}) + \sum_{i_{1} = 1}^{n} \sum_{i_{2} = 1}^{i_{1}} a_{i_{1} i_{2}} H_{2} (ξ_{i_{1}}, ξ_{i_{2}}) \\ + \sum_{i_{1} = 1}^{n} \sum_{i_{2} = 2}^{i_{1}} \sum_{i_{3} = 3}^{i_{2}} a_{i_{1} i_{2} i_{3}} H_{3} (ξ_{i_{1}}, ξ_{i_{2}}, ξ_{i_{3}}) + … \end{array}

(11)

where n is the number of variables, and $a_{0}, a_{i 1}, \dots$ is an undetermined constant term; the number of undetermined coefficients is

N = \frac{(m + n)!}{m! n!}

(12)

The calculation formula of an m-order, n-dimensional Hermite polynomial is

H_{m} (ξ_{i_{1}}, \dots, ξ_{i_{n}}) = {(- 1)}^{m} e^{\frac{1}{2} ξ^{T} ξ} \frac{\partial^{n}}{\partial ξ_{i_{1}} \partial ξ_{i_{2}} \dots \partial ξ_{i_{n}}} e^{- \frac{1}{2} ξ^{T} ξ}

(13)

When the order number of the Hermite polynomial is higher, the simulation accuracy of the output response Y is better. At the same time, the number of undetermined coefficients needed to confirm is greater. The measured results show that when m = 2, the fitting precision of a second-order Hermite polynomial can meet the requirements. Therefore, this study uses second-order polynomial chaos

Y = a_{0} + \sum_{i_{i} = 1}^{n} a_{i, 2} ξ_{i} + \sum_{i = 1}^{n} a_{ii, 2} {(ξ_{i}^{2} - 1)}^{2} + \sum_{i = 1}^{n - 1} \sum_{j > i}^{n} a_{ij, 2} ξ_{i} ξ_{j}

(14)

The process is as follows. Select the configuration points of $ξ$ , determine the sampling points of X, calculate the model of the sample, and in turn, determine the undetermined coefficients in equation (14). To improve the calculation speed for the undetermined coefficients, this study adopts configuration point selection based on the principle of linear independence²³ and combines the characteristics of the reactive power optimization under the GMD. The calculation process is illustrated as a flowchart in Figure 3.

Figure 3.

Flowchart of the collocation points selection method.

This process is as follows: calculate the response values of the configuration points based on the configuration points obtained through the power flow calculation model; in turn, determine the undetermined coefficients of the Hermite polynomial by solving the linear system of equations; and finally, the probability distributions of all variables are obtained according to the combination of the Hermite polynomial and the power flow equation. The equations of n nodes are as follows

{\begin{matrix} P_{i} = P_{Gi} - P_{Li} \\ = U_{i} \sum_{j \in i} U_{j} (G_{ij} \cos δ_{ij} + B_{ij} \sin δ_{ij}) \\ Q_{i} = Q_{Gi} - Q_{Li} - Q_{GICi} \\ = U_{i} \sum_{j \in i} U_{j} (G_{ij} \sin δ_{ij} - B_{ij} \cos δ_{ij}) \end{matrix}

(15)

Because of assuming a uniform geoelectric field, we can determine the reactive power demand value of each transformer and set it to be a PQ node in the power flow calculation.

Multi-objective optimization algorithm

Because of the uncertainty of GIC-Q, this study uses the SRSM to describe the output response as a probability distribution function. Then, we combine the MOPSO algorithm to make the reactive power optimization process closer to actual operation and improve the adaptability and reliability of the optimized scheme under the GMD.

Multi-objective optimization problem

The multi-objective optimization problem can be defined as

{\begin{matrix} min F_{i} (X) = [f_{1} (X), f_{2} (X), \dots f_{M} (X)] \\ s . t . g_{j} (X) \leq 0 (or \geq 0) \\ h_{k} (X) = 0 \end{matrix}

(16)

where $X = [x_{1}, x_{2}, \dots, x_{n}] \in S$ ; n is the number for decision-making; S means a D-dimensional decision space; and M is the target number.

Several definitions associated with this discussion are given below.

Definition 1

Pareto dominance: there are two feasible solutions, which are x₁ and x₂, considering x₁ dominates x₂ to meet all $f (x_{1}) \leq f (x_{2})$ , and there are at least one $f (x_{1}) < f (x_{2})$ , and $i, j = 1, 2, \dots, n$ , write for $x_{1} ≻ x_{2}$ .

Definition 2

The optimal solution of Pareto, a solution that exists, $x * \in Ω$ ( $Ω$ is an M-dimensional object space), the optimal solution or non-dominated solution of Pareto under the circumstance is given by

\neg \exists x \in Ω : x ≻ x^{*}

(17)

Definition 3

The optimal solution set and cutting edge of Pareto. All sets of Pareto solutions are referred to as the Pareto optimal solution set. Namely

P = {x^{*} | \neg \exists x \in Ω : x ≻ x^{*}}

(18)

The area formed by all the Pareto optimal solutions corresponding to the objective function value is referred to as the Pareto frontier or equilibrium surface.

In this study, the GIC-Q is used as the research object. Under the condition where the structure and capacity of the grid have been defined, this study makes full use of the existing adjustable reactive power output of generator and on-load voltage regulating transformer, and at the same time, collocates reactive power compensation equipment at appropriate nodes to deal with voltage fluctuations caused by the GIC-Q to ensure the safe and stable operation of the power system under the GMD. The objective functions and constraint functions of the multi-objective are as follows.

1. Objective function

For the target of the minimum cost of the node voltage deviation and reactive power compensation equipment, this article sets up the multi-objective reactive power optimization function considering the influence of GIC

{\begin{matrix} min f = E (f_{VD}, f_{Q}) \\ f_{VD} = \sum_{i \in N} | \frac{U_{i} - U_{i 0}}{Δ U_{max}} | \\ f_{Q} = \sum_{S \in N_{C}} Q_{S} \end{matrix}

(19)

where $U_{i}$ and $U_{i 0}$ are, respectively, the actual voltage and the nominal voltage of the PQ node; $Δ U_{max}$ is the maximum deviation (namely, the difference value of the upper and lower limit values of the node voltage); N is the number of the PQ node; N_C is the set of adjustable points of reactive power compensation capacity; Q_S is the reactive power compensation capacity; and E(·) is the expected value.

2. Constraint condition

The constraint condition is namely the power flow equation in the process of the probabilistic power flow calculation

f (P_{G}, Q_{G}, P_{L}, Q_{L}, Q_{C}, Q_{GIC}) = 0

(20)

The control variable range constraints are

{\begin{matrix} U_{Gi min} \leq U_{Gi} \leq U_{Gi max}, i = 1, 2, \dots, N_{G} \\ T_{j, min} \leq T_{j} \leq T_{j, max}, j = 1, 2, \dots, N_{T} \\ Q_{CK min} \leq Q_{C} \leq {Q_{C}}_{k max}, k = 1, 2, \dots, N_{C} \end{matrix}

(21)

The state variable range constraints are

{\begin{matrix} P_{r} (U_{i, min} \leq U_{i} \leq U_{i, \max}) > P_{r, \min} \\ Q_{Gi, min} \leq Q_{Gi} \leq {Q_{Gi}}_{, max} \end{matrix}

(22)

where $P_{G} and Q_{G}$ are, respectively, the active power and reactive power of the generator bus; $P_{L} and Q_{L}$ are, respectively, the active power and reactive power of the load bus; Q_GIC is the GIC-Q loss; Q_C is the reactive power offered by a parallel capacitor; $U_{G}$ is the voltage of the generator; $T$ is the transformer tap; $U_{i}$ is the voltage of the PQ bus; $N_{G}$ is the number of the generator bus; $N_{T}$ is the number of the transformer branch; $N_{C}$ is the number of the reactive power compensation bus; $P_{r} (\cdot)$ is the probability; $P_{r, min}$ is the threshold voltage needed to be reached by the voltage-qualified probability; and the subscripts “max” and “min” denote the upper limit and lower limit, respectively.

Particle swarm optimization algorithm

The particle swarm optimization (PSO) algorithm was first proposed by Kennedy. It is a calculation method based on swarm intelligence evolution. The basic principle is searching the optimal solution through iteration. At each iteration, particles tracking individual optimal positions (P_best) and the global optimal position (G_best) update themselves. The speed and position formulas of the particles named i in the kth iteration and the first d components are

v_{i_{d}}^{k + 1} = {wv}_{i_{d}}^{k} + c_{1} r_{1} (p_{i_{d}}^{k} - x_{i_{d}}^{k}) + c_{2} r_{2} (p_{g_{d}}^{k} - x_{i_{d}}^{k})

(23)

x_{i_{d}}^{k + 1} = x_{i_{d}}^{k} + v_{i_{d}}^{k + 1}

(24)

where $w$ is the inertia weight; $c_{1}$ and $c_{2}$ are the accelerated constants that, respectively, adjust the step length of the global optimal particle and individual optimal particle flight directions; $r_{1}$ and $r_{2}$ are uniform random numbers in [0,1]; $p_{id}$ is the optimal location of the individual; $p_{gd}$ is the global optimal position; K is the number of iterations.

Niche chaotic mutation MOPSO algorithm

Because the PSO algorithm has the disadvantages of precocious phenomena and late convergence speed, this article proposes MOPSO²⁴ by considering the niche chaotic mutation.

The niche algorithm forms a niche for each particle based on the distance between the particles. They become a particle swarm, and then a particle swarm algorithm is used for each particle group to update the speed and position values. The updated particle swarm uses the sharing mechanism to improve fitness of the particles, and ensure optimization of the group.

To save the Pareto optimal solution which is found in the process of evolution, an external pre-strike technology is established in the process of iteration to adjust the external archive and improve the computational efficiency. This study uses the crowded distance sorting archived external set of particles.²⁵

The distance between the particles and the fitness is determined according to the following equations. The distance between the particles is

d_{ij} = ‖ X_{i} - X_{j} ‖ = \sqrt{\sum_{k = 1}^{n} {(x_{ik} - x_{jk})}^{2}}

(25)

where $X_{i}$ and $X_{j}$ are the two particles of the D-dimensional decision space. For a given parameter (niche radius, $R_{ch}$ ), if $R_{ch} > d_{ij}$ , the individual will join the niche groups called X_p.

Fitness is

f'_{i} = \frac{f_{i}}{\sum_{j = 1}^{n} s (d_{ij})}

(26)

where $f_{i}$ is the original fitness of an individual of the first i; n is the individual number of the niche; and $s (d_{ij})$ is a sharing function

s (d_{ij}) = {\begin{matrix} 1 - {(d_{ij} / R_{h})}^{λ} & d_{ij} < R_{h} \\ 0 & else \end{matrix}

(27)

where $λ$ is the parameter controlling the sharing function shape.

A too fast convergence speed will make the MOPSO algorithm to search for a local optimum, which is the Pareto optimal frontier instead of the global optimal frontier. For this reason, the chaotic mutation is brought in MOPSO in this study, which generates a chaotic disturbance for the current particle to return the best individual and to bring the solution out of the local extremum interval. Here, the chaotic mutation mapping logistic iterative formula is given by

β_{j}^{k + 1} = μ β_{j}^{k} (1 - β_{j}^{k})

(28)

where k is the number of iterations; $β_{j}$ is the first j chaotic variables of the corresponding particle, named $X_{i}$ ; $β_{j} \in (0, 1)$ ; and $μ \in (2, 4]$ .

Flowchart of the niche chaotic mutation MOPSO is shown in Figure 4.

Figure 4.

Flowchart of the niche chaotic mutation MOPSO.

Example analysis

Example overview

The simulation calculation of this study uses a GIC standard test system, as shown in Figure 5.

Figure 5.

Diagram of the GIC-benchmark case.

The system includes buses, generators, switching stations, transmission lines, and transformers in the quantities of 18, 7, 1, 15, and 15, respectively. Bus 1 is selected as a balancing bus, the load of which is 4700 MW+1800 MVar; other parameters can be found in Pirjola²⁶ and Liu and Wu.²⁷

GIC-Q calculation

In the simulation, GMD caused by the solar activity impacts the current systems in the mid-low-latitude ionosphere. At present, no research has yet proven the direction characteristics of the current systems of the mid-low-latitude ionosphere. It is thought that the geoelectric fields at the mid-low-latitude levels are 1 V/km and can occur in any direction.²⁸

The GICs through the transformer windings and transmission lines are calculated for directions $θ \in [0 \circ, 360 \circ]$ , with an increment of 5°, and under the 1 V/km uniform electric field. Figure 6 shows the calculation results of the HV winding single-phase GIC in the T6 transformer for 24 different directions, with the increment of 15°. With East as the reference direction, the direction of the electric field is shown using the direction of the line; the magnitude of GIC is expressed by the length of the line.

Figure 6.

GIC at high-voltage winding for different directions of geoelectric fields.

According to the basic principle of GIC calculation and the calculation results of Figure 6, the direction of the electric field between [0°, 360°] of GIC is a 180° symmetry, which means the magnitude of GIC is equal and the polarity of GIC is opposite. The loss value of GIC-Q is related only to the magnitude of the GIC and is not related to polarity. For this reason, while calculating the value of GIC-Q, we considered the direction of the electric field between [0°, 180°]. The calculation results are shown in Figure 7.

Figure 7.

GIC-Q losses for different directions of geoelectric fields.

It is known by the calculation that the GIC-Q achieves a maximum value of 914 MVar with the direction of geoelectric fields at 135°. The reactive power loss of each transformer is illustrated in Table 1.

Table 1.

Reactive power loss of transformers.

T_i	2	3	4	5	6	7	8
Q_T (MVAR)	26.9	53.9	53.9	172.2	72.8	72.8	58.5
T_i	9	10	11	12	13	14	15
Q_T (MVAR)	58.5	13	13	59.7	26.9	59.7	172.2

With the calculated value of GIC-Q at each transformer in Table 1 as the mean value, 10% of the mean value as the standard deviation, and the uncertainty of GIC-Q described by the normal distribution function, the specific expression is shown in equation (6).

Figure 8 shows the voltage cumulative distribution function of node-20 through the SRSM and Monte Carlo method (MCM).

Figure 8.

Cumulative distribution function of the node-20 voltage amplitude.

Figure 8 shows that the cumulative distribution function obtained through the SRSM is the same as that obtained with the MCM. However, the computing speed and precision of the SRSM are superior, thus verifying the effectiveness of the SRSM.

Consequently, the limit probability value of each node voltage below 0.95 pu can be obtained using the SRSM probability flow model proposed in this article, which is shown in Table 2.

Table 2.

Low-voltage limit probability (%) of the node.

Bus	1	2	3	4	5	6	7	8	11
Probability value	0	12.1	45.2	100	87.6	64.3	53.9	50.2	46.2
Bus	12	13	14	15	16	17	18	19	20
Probability value	32.1	13.8	10.5	100	78.4	52.9	10.8	12.1	100

It is seen from Table 2 that under the effect of a uniform geoelectric field of 1 V/km, the limit risk of nodes 4, 15, and 20 is up to 100%, which seriously threatens the safe and stable operation of the power grid, and therefore, the reactive power must be optimized.

Multi-objective optimization calculation

The process of this calculation is as follows: ensure the largest reactive power output of the existing generator and the largest adjustment range of the on-load regulating transformer; select nodes 4, 15, and 20 of the risk of the low-pressure limit to the parallel capacitor group; and proceed with reasonable reactive power optimization to make each node voltage within a reasonable range and ensure the safe and stable operation of power grid.

In this study, we use the SRSM method to establish the probability calculation model of reactive power optimization. The node voltage offset and reactive power compensation equipment’s minimum cost are the optimization goals. Pareto optimal solutions of these two targets are obtained through niche chaotic MOPSO.

Parameter settings of the niche chaotic MOPSO algorithm are as follows: inertia weight, $w = 0.84$ ; accelerating factor, $c_{1} = c_{2} = 1.545$ ; the scale of the niche archive is 70; the number of population particles is 50; and the maximum number of iterations is 100.

To illustrate the effectiveness of the optimization algorithm in this study to increase the solution of niche MOPSO and MOPSO, with the largest number of iterations as calculation termination conditions, the convergence process and convergence speed of the three types of algorithms are compared, as shown in Figure 9.

Figure 9.

Diagram of the iteration convergence.

It is shown in Figure 9 that the rate of convergence of the early stage of the iterations and the number of optimal solution searches of the niche chaotic mutation MOPSO are significantly faster and more than niche MOPSO and MOPSO, respectively. The individual number of Pareto optimal solution sets has reached the upper limit when the iteration is 27 and realizes the algorithm convergence. Niche MOPSO and MOPSO, respectively, realized the algorithm convergence when the iterations were 38 and 44, respectively. By contrast, the search ability and search speed of the algorithm proposed in this article are higher than niche MOPSO and MOPSO.

The target is the minimum expectations of offset voltage and reactive power compensation capacity, using niche chaotic MOPSO to calculate the Pareto frontier optimal solution set, which is shown in Figure 10.

Figure 10.

Diagram of the distribution of Pareto optimal solution set.

It is shown in Figure 10 that the optimal solution set is searched by the niche chaotic MOPSO, which forms the optimal solution of multiple optimization objectives in different situations. It is also known by the calculation that the optimal compromise solution of niche chaotic MOPSO is 3.67 pu and 960 MVar. The optimal solution of the optimization of the proceeds is that nodes 4, 15, and 20, respectively, input capacitor banks whose capacities are 210, 350, and 400 MVar. After the reactive power compensation optimization, node 20, whose voltage expectation is 0.964 pu, and its limited risk increase by 44.15%.

Conclusion

Using an uniform electric field of 1 V/km, this study calculated the GIC-Q of GIC-benchmark case in the direction of the electric field between [0°, 180°] with a 5° increment. The calculation results show that the GIC-Q reaches a maximum value of 914 MVar when the direction of the electric field is 135°. Such reactive power loss will inevitably cause a power grid reactive power shortage. For this reason, the reactive power compensation equipment should be reasonably collocated in the power grid and optimized.

Considering the uncertainty of GIC-Q of the power grid, this study set up a power system probabilistic power flow calculation model including the GIC-Q using the SRSM. The calculated results were compared with the MCM, showing that the SRSM used to address the uncertain random issue remains the same level as the MCM, but has the advantages of high precision and fast computing speed.

With the minimum expectations of the offset voltage and reactive power compensation capacity as the targets, this study established a multi-objective reactive power optimization calculation model and then solved the reactive power optimization model through the niche chaotic MOPSO algorithm. The example calculation results show that the convergence rate of the proposed method is faster and the optimal solution set searched is closer to the optimal frontier of ideal Pareto; thus, it has better optimization performance. The optimal compromise solution given by the niche chaotic MOPSO algorithm provided a reactive power optimization configuration that can improve the grid voltage dramatically and ensure the safe and stable operation of power grid under the GMD.

Footnotes

Academic Editor: ZW Zhong

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Takasu

Oshi

Miyawaki

et al . An experimental analysis of DC excitation of transformers by geomagnetically induced currents. IEEE T Power Deliver 1994; 9: 1173–1183.

Kappenman

. Geomagnetic storms and their impact on power system. IEEE Power Eng Rev 1996; 16: 5–8.

Barnes

Van Dyke

. Economic consequences of geomagnetic storms (a summary). IEEE Power Eng Rev 1990; 10: 3–4.

Lesher

Porter

Byerly

. Sunburst—a network of GIC monitoring systems. IEEE T Power Deliver 1994; 9: 128–137.

Wen

Liu

Xiang

et al . Influences of geomagnetic induced currents on security and stability of power systems. Power Syst Tech 2010; 34: 24–30.

Pirjola

. Geomagnetically induced currents during magnetic storms. IEEE T Plasma Sci 2000; 28: 1867–1873.

Berge

Varma

Marti

. Laboratory validation of the relationship between geomagnetically induced current (GIC) and transformer absorbed reactive power. In: Electrical power and energy conference, Winnipeg, MB, Canada, 3–5 October 2011, pp.491–495. New York: IEEE.

Marti

Berge

Varma

. Determination of geomagnetically induced current flow in a transformer from reactive power absorption. IEEE Trans Power Deliv 2013; 28: 1280–1288.

Liu

Qin

. GIC reactive power loss calculation based on U-I curve for single-phase autotransformer. Electr Power Autom Equip 2015; 35: 55–60.

10.

Etemadi

Rezaei-Zare

. Optimal placement of GIC blocking devices for geomagnetic disturbance mitigation. IEEE T Power Syst 2014; 29: 2753–2762.

11.

Liu

. An optional control method for overcoming clustering transformer reactive power losses due to geomagnetic storm. Autom Electr Power Syst 2014; 38: 32–37.

12.

Berhardi

Cilliers

Gaunt

. Improvement in the modelling of geomagnetically induced currents in southern Africa. S Afr J Sci 2008; 104: 265–273.

13.

Trivedi

Vitorello

Kabata

et al . Geomagnetically induced currents in an electric power transmission system at low latitudes in Brazil: a case study. Space Weather 2007; 5: S04004.

14.

Watari

Kunitake

Kitamura

et al . Measurements of geomagnetically induced current in a power grid in Hokkaido, Japan. Space Weather 2009; 7: S03002.

15.

Liu

Xie

. Analysis of increase of noise of 500 kV transform. High Volt Eng 2005; 31: 85–87.

16.

Liu

Pirjola

. Geomagnetically induced currents in the high voltage power grid in China. IEEE T Power Deliver 2009; 24: 2368–2374.

17.

Liu

Niu

. Disastrous space weather risk on large-scale power grid. In: Proceedings of the 2010 5th international conference on critical infrastructure (CRIS), Beijing, China, 20–22 September 2010, pp.1–5. New York: IEEE.

18.

Liu

. Influence and hazard of disastrous space weather on power grid in China. Eng Sci 2010; 12: 29–33.

19.

Kappenman

Albertson

Damsk

et al . Solar wind monitor satellite. IEEE Power Eng Rev 1990; 10: 4–8.

20.

Kappenman

Radasky

Gilbert

et al . Advanced geomagnetic storm forecasting: a risk management tool for electric power operations. IEEE T Plasma Sci 2000; 28: 2114–2121.

21.

Zheng

Liu

David

et al . Modelling geomagnetically induced currents in multiple voltage levels of a power system illustrated using GIC-Benchmark case. P CSEE 2013; 33: 179–187.

22.

Rezaei-Zare

. Behavior of single-phase transformers under geomagnetically induced current conditions. IEEE T Power Deliver 2014; 29: 916–925.

23.

Jiang

Zhou

. Optimal probabilistic collocation points for stochastic response surface method. Chin J Comput Mech 2012; 29: 345–351.

24.

Luo

Yang

et al . Layout problem of Multi-component systems arising for improving maintainability. J Cent South Univ T 2014; 21: 1833–1841.

25.

Zhu

. Cultural based multi-objective particle swarm optimization algorithm using crowding distance sorting method. Control Decis 2012; 27: 1406–1410.

26.

Pirjola

. Review on the calculation of surface electric and magnetic fields and of geomagnetically induced currents in ground-based technological systems. Surv Geophys 2002; 23: 71–90.

27.

Liu

. Security risk assessment ideals and theoretical framework for power system considering geomagnetic storm. P CSEE 2014; 34: 1583–1591.

28.

Liu

Zhang

et al . Strong magnetic storm’s influence on China’s Guangdong power grid. Chin J Geophys 2008; 51: 976–981.