The Static Stability Region of an Integrated Electricity-Gas System Considering Voltage and Gas Pressure

Yunfei Mu; Zhibin Liu; Xiangwei Guo; Hongjie Jia; Kai Hou; Xiaodan Yu; Bofeng Luo; Hairun Li

doi:10.1016/j.eng.2023.12.009

Engineering ›› 2024, Vol. 36 ›› Issue (5) :139 -153. DOI: 10.1016/j.eng.2023.12.009

Research

Article

The Static Stability Region of an Integrated Electricity-Gas System Considering Voltage and Gas Pressure

Author information +

History +

PDF (3040KB)

Abstract

In an integrated electricity-gas system (IEGS), load fluctuations affect not only the voltage in the power system but also the gas pressure in the natural gas system. The static voltage stability region (SVSR) method is a tool for analyzing the overall static voltage stability in a power system. However, in an IEGS, the SVSR boundary may be overly optimistic because the gas pressure may collapse before the voltage collapses. Thus, the SVSR method cannot be directly applied to an IEGS. In this paper, the concept of the SVSR is extended to the IEGS-static stability region (SSR) while considering voltage and gas pressure. First, criteria for static gas pressure stability in a natural gas system are proposed, based on the static voltage stability criteria in a power system. Then, the IEGS-SSR is defined as a set of active power injections that satisfies multi-energy flow (MEF) equations and static voltage and gas pressure stability constraints in the active power injection space of natural gas-fired generator units (NGUs). To determine the IEGS-SSR, a continuation MEF (CMEF) method is employed to trace the boundary point in one specific NGU scheduling direction. A multidimensional hyperplane sampling method is also proposed to sample the NGU scheduling directions evenly. The obtained boundary points are further used to form the IEGS-SSR in three-dimensional (3D) space via a Delaunay triangulation hypersurface fitting method. Finally, the numerical results of typical case studies are presented to demonstrate that the proposed method can effectively form the IEGS-SSR, providing a tool for IEGS online monitoring and dispatching.

Graphical abstract

Keywords

Integrated electricity-gas system / Static stability region / Continuation multi-energy flow / Multidimensional hyperplane sampling

Cite this article

Download citation ▾

Yunfei Mu, Zhibin Liu, Xiangwei Guo, Hongjie Jia, Kai Hou, Xiaodan Yu, Bofeng Luo, Hairun Li. The Static Stability Region of an Integrated Electricity-Gas System Considering Voltage and Gas Pressure. Engineering, 2024, 36(5): 139-153 DOI:10.1016/j.eng.2023.12.009

登录浏览全文

4963

注册一个新账户忘记密码

1. Introduction

In the context of global carbon-reduction efforts, renewable energy generation is widely used in power systems [1], [2], [3]. The very uncertain power output of this type of energy generation necessitates generators that can quickly respond to severe power fluctuations [4]. Natural gas-fired generator units (NGUs) have become one of the primary choices for use with renewable energy generation in a power system due to their low pollution, fast response speed, and high efficiency [5], [6], which increase the integration of power systems and natural gas systems. The concept of the integrated electricity-gas system (IEGS) has been proposed for the unified analysis of power systems and natural gas systems.

With the expansion of the coupling between these two energy systems, security issues caused by interactions occur frequently and have become the main challenge to the security and stability operation of IEGSs [7], [8]. On the one hand, a short supply of natural gas directly affects the output of NGUs and may lead to power system blackouts. On the other hand, the load fluctuations in the power system affect the security of natural gas pipelines and nodal gas pressures [5]. For example, in recent years, several power outages have occurred in Taiwan Province in China [9], southern California in United States [10], and the United Kingdom [11], arousing public concern worldwide. All these accidents are attributed to the natural gas supply. Therefore, ensuring the stability of both systems of an IEGS at the same time is very significant.

In the field of steady-state stability analysis, studies on power system static stability analysis date back to the 20th century [12], [13]. The solution of power flow equations is generally acknowledged to be the approximate equilibrium point of a power system and, when there is no solution to the power flow equations, the loss of power system static voltage stability is assumed [14], [15]. Several indexes and methods have been proposed to measure the static voltage stability of power systems, such as the sensitivity index [16], singular value index [17], and load margin index [18]. Similar methods are further used in IEGSs for static stability analysis [19], [20]. Most of these studies are based on a “pointwise method,” which means that only the stability of one specific operating point can be estimated at a time.

The region-based theory is widely used in the power system analysis [21], especially in static voltage stability analysis to obtain the overall static stability of the system. The static voltage stability region (SVSR) concerning the static voltage stability and load margin under uncertain load growth is developed in Ref. [22], where the SVSR is defined as a region that is composed of all static voltage stable operating points under a certain network topology and parameters. The SVSR can be defined in the power injection space [23], the cut-set complex power space [24], and other parameter spaces [25]. The continuation power flow (CPF) method is still the most common method used to search for the SVSR and its boundary [26]. Based on the CPF method, several new methods using geometric parameter adjustment [27], arc-length parametrization, and piecewise approximants [28] have been proposed to speed up the calculation of the SVSR.

However, regarding IEGSs, very few studies have focused on a unified static stability analysis considering both the voltage and gas pressure. Owing to the significant differences in modeling power systems and natural gas systems, traditional static voltage stability analysis methods for power systems cannot be directly applied to natural gas systems and IEGSs. In fact, the natural gas system is considered to be the source of the power system with the coupling of NGUs. An increase in the power output of the NGUs further calls for an increase in gas demand, which affects the power and gas flow distribution in an IEGS [29], leading to a dramatic decrease in gas pressure [30]. The power output of the NGUs no longer increases due not only to the restrictions in the power system but also to those in the natural gas system—more specifically, regarding the gas pressure. A unified static stability analysis is needed to avoid overoptimization in the estimation and measurement of IEGS static stability while considering voltage and gas pressure. A clearer understanding of the restrictions set by static voltage and gas pressure stability can assist system operators in their decision-making for IEGS security monitoring.

Thus, in this paper, the static stability region (SSR) of an IEGS considering voltage and gas pressure (referred to as an IEGS-SSR) is proposed to ensure security monitoring with load fluctuations. The primary differences between the IEGS-SSR and existing works are encapsulated in two key aspects. The first aspect pertains to the differentiation between “region” and “margin” in the static stability analysis. In a static stability analysis, the “margin” signifies the distance between the current operating point and the collapse point of the IEGS [20], [31], [32]. The IEGS-SSR can determine a set of static stable operating points by considering voltage and gas pressure, providing a comprehensive perspective on the system’s static stability margin while taking voltage and gas pressure into account. The second difference is that existing static stability analyses in an IEGS primarily focus on the impact of other coupling systems on the static voltage stability within a power system [20], [32]. The IEGS-SSR can conduct a unified static stability analysis on both the power system and the natural gas system, while also investigating the reciprocal influence of the power system on the natural gas system.

The main contributions of this work are as follows:

(1) Based on a steady-state energy flow model of a natural gas system, the static gas pressure stability criterion of a natural gas system is developed. A concept and definition of the natural gas system static gas pressure stability region (SGSR) and its boundary are then proposed. Combined with the concept and restraints of power system static voltage stability, the IEGS-SSR and its boundary are further defined.

(2) A continuation multi-energy flow (CMEF) method extended from the CPF method is proposed in this paper to determine the boundary point of the IEGS under a specific NGU scheduling direction. Moreover, a multidimensional hyperplane sampling method is developed to obtain enough NGU scheduling directions. The IEGS-SSR is then constructed by sampling boundary points through hypersurface fitting. The proposed method can evenly form the IEGS-SSR with a low fitting error in engineering.

The remaining sections of this paper are organized as follows: Section 2 describes the steady-state multi-energy model of the IEGS, the corresponding multi-energy flow (MEF) calculation method, and the definition of the IEGS-SSR based on that model. Section 3 introduces the determination method of the IEGS-SSR. A case with several scenarios is studied in Section 4. Finally, Section 5 presents the main conclusions and discusses future work.

2. The IEGS-SSR

This section first presents the MEF model of the IEGS and the corresponding calculation method. The definition and constraints of the SVSR are then given and further extended to the SGSR in a natural gas system and the IEGS-SSR.

2.1. MEF model and calculation method

A static voltage stability analysis is used to describe the ability of a power system to maintain the voltage of all load buses without voltage collapse after disturbances such as load growth [33], [34]. When power load growth occurs in a power system, the power system operates from the initial equilibrium state. A static voltage stability analysis captures the power system’s operational state and evaluates whether this operational state after the power load growth is an equilibrium state. If the static voltage stability analysis shows no voltage collapse, the power system will finally reach a new equilibrium state (i.e., steady state) after the power load growth. During a static stability analysis, the dynamic process of power transmission to the load buses is overlooked. Although the operational time scale of a natural gas system is much longer than that of a power system, an IEGS static stability analysis—analogous to a static voltage stability analysis—captures the IEGS’s operational state and judges whether or not the IEGS will reach a new equilibrium state (i.e., steady state) after power or gas load growth. The dynamic process of gas flow rate transmission through the pipelines is overlooked as well. Thus, in the field of the static analysis in the IEGS, the dynamic process of the dynamic elements in the IEGS can be ignored [35].

A steady-state MEF model is presented in this section to calculate the state variables in an IEGS for a unified static stability analysis, including a power system model, a natural gas system model, and an NGU model.

2.1.1. Power system model

The power system model is described in Eqs. (1), (2).

(1)

P_{i}^{sp} - V_{i} \sum_{j \in i} V_{j} (G_{ij} \cos θ_{ij} + B_{ij} \sin θ_{ij}) = 0

(2)

Q_{i}^{sp} - V_{i} \sum_{j \in i} V_{j} (G_{ij} \sin θ_{ij} - B_{ij} \cos θ_{ij}) = 0

These equations represent the active and reactive power balance of each bus, where i, and j represent the buses in a power system; V_i and V_j are the voltage magnitude of buses i and j, respectively; G_ij and B_ij are the conductance and susceptance of the transmission line between buses i and j, respectively; θ_ij is the voltage phase angle difference between buses i and j;

P_{i}^{sp}

and

Q_{i}^{sp}

are the active and reactive power injection of bus i, respectively, as given in Eqs. (3), (4),

(3)

P_{i}^{sp} = P_{G, i} - P_{l, i}

(4)

Q_{i}^{sp} = Q_{G, i} - Q_{l, i}

where P _G _,i and Q _G _,i are the active and reactive power generation of bus i, respectively; P _l _,i and Q _l _,i are the active and reactive power load of bus i, respectively.

2.1.2. Natural gas system model

A model based on a node-branch method similar to the power system is applied in this part for the gas flow calculation. The variables of the natural gas system are the gas flow injection and the nodal gas pressure [36]. The nodes in the natural gas system can be divided into two types: slack nodes (also known as gas source nodes) and load nodes. The node classification and variables in the natural gas system are shown in Table 1.

Note: √ indicates that the variable is known at the node and × indicates that it is unknown at the node.

The natural gas flowing from node m to node n depends on the pressure difference between two nodes and the parameters of pipeline mn. The Weymouth, Panhandle A, and Panhandle B flow equations are commonly used to analyze flow through pipelines. Unifying these three equations, the gas flow through pipeline f_mn is calculated using Eqs. (5), (6) [37], [38], [39].

(5)

c_{mn} f_{mn}^{ξ} = s_{mn} \cdot (p_{m}^{2} - p_{n}^{2})

(6)

s_{mn} = \{\begin{matrix} + 1 p_{m} \geq p_{n} \\ - 1 p_{m} < p_{n} \end{matrix}

where c_mn is the resistance coefficient of the pipeline mn, which is related to the roughness, diameter, and pipeline length; and ξ refers to the flow index. Both c_mn and ξ are coefficients related to the gas pressure level of the system; p_m and p_n are the gas pressure of nodes m and n, respectively; function s_mn is the sign function.

At any node in the natural gas system, the inflow and outflow flows must be balanced, as shown in Eqs. (7), (8).

(7)

L_{m}^{sp} - \sum_{n \in Ω_{m}} f_{mn} = 0

(8)

L_{m}^{sp} = L_{s, m} - L_{l, m}

where

L_{m}^{sp}

is the natural gas injection at node m; Ω_m is the set of nodes connected with node m through pipelines; L _s _,m and L _l _,m are the natural gas supply and load at node m, respectively.

2.1.3. NGU model

NGUs generate electricity by consuming natural gas; for simplicity, the relationship between their natural gas consumption and the power output is expressed as a linear function, as shown in Eq. (9).

(9)

P_{G, r} = H η_{r} L_{G, r}

where P_G_,r is the active power generation of the rth NGU; H is the heat value of natural gas, which is 1.0833 × 10⁻² MW·h·m⁻³;

η_{r}

is the generating efficiency of the rth NGU, and the typical value is 0.55; and L _G _,r is the natural gas consumption of the rth NGU [37].

2.1.4. MEF calculation model

We define x_eg and y_eg as the state variables and injection variables of the IEGS. Composed of Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9), the MEF calculation model is expressed as Eqs. (10), (11), (12).

(10)

F_{eg} (x_{eg}, y_{eg}) = 0

(11)

x_{eg} = [V, θ, p]

(12)

y_{eg} = [P^{sp}, Q^{sp}, L^{sp}]

where function F _eg() represents the MEF equations in Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9); V, θ, and p are the vectors of voltage magnitude, voltage phase angle, and gas pressure, respectively; P ^sp, Q ^sp, and L ^sp are the vectors of active power injection, reactive power injection, and natural gas injection, respectively.

The calculation method can be divided simply into two types: the unified-solving method and the decoupling-solving method [37]. The unified-solving method means the MEF is calculated via the Newton–Raphson method with a unified Jacobian matrix J_eg. The decoupling-solving method decouples the IEGS from the coupling units and then calculates the MEF separately. In this paper, a decoupling-solving method is used for MEF calculation. The power systems and natural gas systems are decoupled from the NGUs. The output of the NGUs is added to the active power generation of the linked PV buses (the voltage magnitude and active power injection of the bus are given) in the power system, and the gas flow input is added to the load demand of the linked load node in the natural gas system. The MEF is then calculated using the alternative iteration method, as shown in Ref. [37].

2.2. Definition of the IEGS-SSR

The security region of an IEGS has been described in Refs. [40], [41]. It can be defined as a closed region in which the operating point satisfies the constraints of the MEF equations and other operational constraints. Depending on the specific operational constraints a security region must adhere to, the region can be categorized into different types, such as voltage security regions, pressure security regions, and so on. The security region is widely used in economic dispatch and optimal security control in an IEGS. However, during these dispatch and control processes, there may be instances when the system’s operational constraints are temporarily unmet. In such cases, it is crucial to maintain the solvability of the MEF and prevent voltage and gas pressure collapse within an IEGS model, similar to a power system [42]. Therefore, a new region that only considers MEF equation constraints (equivalent to the static stability constraints of the voltage and gas pressure in an IEGS) should be proposed to assess the static stability in an IEGS, while taking voltage and gas pressure into account. This region is referred to as the IEGS-SSR. Fig. 1 shows the relationship between the IEGS security region and the IEGS-SSR. The IEGS-SSR is a region in which the operating points satisfy the MEF equations and the static voltage and gas pressure stability constraints in an IEGS. In comparison with the IEGS-SSR, the IEGS security region is subject to additional operational constraints, such as voltage magnitude constraints and pressure security constraints [43]. Operating points within the IEGS-SSR may violate these additional operational constraints and extend beyond the boundary of the IEGS security region. Consequently, the IEGS security region is encompassed within the IEGS-SSR.

When the NGU injections ω _eg is determined, the state of the IEGS is uniquely determined. Thus, the IEGS-SSR Ω _eg can be described as the set of ω _eg that has both static voltage and pressure stability. To obtain the IEGS-SSR, the regions and boundary criteria of the SVSR and SGSR should be defined first. Then, the criteria of IEGS-SSR are determined.

The static voltage stability is used to describe the ability of the power system to maintain the voltage of the load bus near the equilibrium point after disturbances such as rapid load growth. The single-machine single-load power system illustrated in Fig. 2(a) is usually presented to explain how the voltage collapses. The P-V curve shown in Fig. 2(b) indicates the downtrend of V _end along with the load growth [20], [44], [45] and finally reaches its lower limit at point C, the nose point. The power flow at point C reaches its maximum [27]. Moreover, the power flow Jacobian matrix J _ee at point C is singular, which is described in Eq. (13).

(13)

\det (J_{ee}) = 0

Here, we define e_min as the minimum eigenvalue of J_ee, which is depicted by Eq. (14).

(14)

e_{min} = \min σ (J_{ee})

where function σ() is the spectral decomposition function.

The e_min monotonically approaches zero near nose point C with an increasing power load. In engineering practice, an e_min less than a threshold ε (where ε is usually set to 10^–5) is considered to represent a loss of the static voltage stability of a power system. Thus, the boundary criteria of the power system can be described as shown in Eq. (15) [46].

(15)

e_{min} \geq ε

The SVSR Ω_e is described as the set in Eq. (16).

(16)

Ω_{e =} \{\forall ω_{e} \in R^{1 \times N} | F_{e} (x_{e}, y_{e}) = 0, e_{min} \geq ε\}

The boundary is shown in Eq. (17).

(17)

\partial Ω_{e} = \{\forall ω_{e} \in R^{1 \times N} | F_{e} (x_{e}, y_{e}) = 0, e_{min} = ε\}

where ω_e is the NGU injections vector in a power system; R is the set of all real numbers; function F_e() represents the Eqs. (1), (2), (3), (4) in a power system; x_e and y_e are the state variables and injection variables of the power system, respectively.

Similar to the static voltage stability in a power system, the static gas pressure stability assesses the ability of a natural gas system to withstand gas pressure collapses following disturbances, such as a rapid growth in gas load. During IEGS operations, a sharp increase in power demand within the power system leads to a corresponding rise in power generation from the NGUs. The substantial gas demand from the NGUs may disrupt the gas flow balance equation in the natural gas system, potentially putting the natural gas system at risk of static gas pressure instability, as exemplified in the case of southern California in August 2020 [47], [48].

Based on the energy circuit theory of the IEGS proposed in Ref. [35], natural gas system pipelines can be equivalent to flow circuits, such as the electric circuit equivalent from power transmission lines. Using Thevenin’s theorem, for a gas flow load at a node, the “outside” natural gas system can be equivalent to one source with a pipeline connected to the load node. Thus, similar to a power system, a simplified single-pipeline single-load natural gas system model is illustrated in Fig. 2(c) to show how the gas pressure collapses through load growth.

Supporting p₀ is the pressure of the sending node, and the pressure of the receiving node p_end is derived using Eqs. (5), (6) and depicted as Eq. (18).

(18)

p_{end} = \sqrt{p_{0}^{2} - L_{end}^{ξ} \cdot c}

where L_end is the natural gas load at receiving node and c is the resistance coefficient.

The process of static gas pressure instability in a natural gas system can be defined similarly to static voltage instability in a power system [33]. When the natural gas load (L _end) increases, the gas pressure (p _end) begins to decrease from its initial equilibrium state, gradually reaching 0 bar (gauge pressure; 1 bar = 10⁵ Pa) [49]. At this point, the natural gas system reaches its maximum transmission gas flow rate and approaches the boundary of the static gas pressure stability operating condition. The gas pressure will no longer decrease when it reaches 0 bar because the reverse pressure differential between the node and the external environment prevents the node from supplying more natural gas to meet the load demand. The static gas load margin will no longer increase.

Determining the gas pressure involves a couple of positive and negative values, as shown in Eq. (18). A gas flow-pressure (L-p) curve of the natural gas system, analogous to the P-V curve of the power system, is illustrated in Fig. 2(d).

With an increasing natural gas load (L _end), the solution of the gas pressure approaches the inflection point C′. of the L-p curve. The gas pressure at C′. is zero. At that time, if the gas load continues to increase, the gas pressure (p _end) ultimately becomes an imaginary number, and the gas pressure solution no longer makes sense. Similar to the nose point C in the power system P-V curve, point C′ is the maximum gas flow (L _max) (also referred to as the critical flow in Ref. [50]) of the natural gas system.

Furthermore, for an unsimplified natural gas system, the natural gas system becomes unstable in terms of gas pressure when the minimum gas pressure p _min in the natural gas system decreases to zero, which is defined as gas pressure collapse in this paper.

As a result, whether or not the natural gas system remains static gas pressure stable is judged by the boundary criteria depicted in Eq. (19).

(19)

p_{min} \geq 0

Thus, the SGSR Ω_g is described by the set in Eq. (20).

(20)

Ω_{g} = \{\forall ω_{g} \in R^{1 \times N} | F_{g} (x_{g}, y_{g}) = 0, p_{min} \geq 0\}

The boundary is shown in Eq. (21).

(21)

\partial Ω_{g} = \{\forall ω_{g} \in R^{1 \times N} | F_{g} (x_{g}, y_{g}) = 0, p_{min} = 0\}

where ω_g is the natural gas injections vector in a natural gas system (converted into equal active power injections); function F_g() represents the Eqs. (5), (6), (7), (8) in a natural gas system; x_g and y_g are the state variables and injection variables of the natural gas system, respectively.

By a combined analysis of the IEGS, when voltage collapse and gas pressure collapse are taken as the criteria of the static stability of the IEGS, the IEGS-SSR can be defined as a set Ω_eg, as shown in Eq. (22).

(22)

Ω_{eg} = \{\forall ω_{eg} \in R^{1 \times N} | F_{eg} (x_{eg}, y_{eg}) = 0, e_{min} \geq ε, p_{min} \geq 0\}

Moreover, the boundary is defined as Eq. (23).

(23)

\partial Ω_{eg} = \{\forall ω_{eg} \in R^{1 \times N} | F_{eg} (x_{eg}, y_{eg}) = 0, e_{min} = ε or p_{min} = 0\}

When the operation point is within the boundary of the IEGS, the boundary criteria are satisfied through MEF equation calculations, and the system is considered to have static stability considering the voltage and gas pressure of the IEGS. In contrast, when the operation point is outside the boundary, the boundary criteria are no longer satisfied, and the system loses the static stability of the IEGS.

3. Determination of the IEGS-SSRs

With a specific network topology of an IEGS, the IEGS-SSR varies with respect to the load growth directions and the selection of NGUs for scheduling, similar to an SVSR in a power system [51], [52]. In this paper, the network topology and the direction of load growth are considered to be constants. Once the NGUs are chosen, the IEGS-SSR is determined by taking into account the variations in NGU scheduling directions. The determination of the IEGS-SSR depends on three issues. One issue is how to determine the boundary point on a given load growth direction, and another is how to determine all possible NGU scheduling directions. After the boundary points under these possible NGU scheduling directions are determined, the third issue is how to form the boundary hypersurface of the IEGS-SSR, as shown in Fig. 3. In this paper, the IEGS-SSR is defined in the active power injection space of the chosen NGUs that meet the load growth demand. The load growth in the natural gas system is converted into equal active power growth through the parameters of the NGUs to keep the results uniform.

In the following part of this section, a CMEF method and a multidimensional hyperplane sampling method are proposed to trace the boundary point in one specific NGU scheduling direction and to obtain numerous possible NGU scheduling directions. The results of the sampling boundary points are used to form the IEGS-SSR through a hypersurface fitting method based on the Delaunay triangulation method.

3.1. The CMEF method

For the power system static voltage stability analysis, the CPF method [26] is proposed to track the drop in voltage at the load bus along with the load growth. Analogous to this method, a CMEF method is proposed in this section that can calculate the varying p _min and e _min when tracking the solution manifold of MEF. A growth parameter, λ, is added to the MEF model to represent the growth ratio of the load or multi-energy generation. The prediction process and the correction process are used for MEF calculation while overcoming the difficulties in MEF convergence when the voltage or gas pressure nearly collapses. Using the criteria given by Eq. (22), the boundary point of the IEGS-SSR is determined.

3.1.1. Growth parameter λ

Adding the growth parameter λ into the MEF model, the load and energy generation variables in Eqs. (3), (4), and (7), (8), (9) are further modified, as depicted in the parts of Eqs. (24a), (24b), (24c), (24d).

(24a)

P_{l, i} (λ_{e}) = P_{l, i 0} + λ_{e} K_{l, i}^{P}

(24b)

Q_{l, i} (λ_{e}) = Q_{l, i 0} + λ_{e} K_{l, i}^{Q}

(24c)

L_{l, m} (λ_{g}) = L_{l, m 0} + λ_{g} K_{l, m}^{L}

(24d)

P_{G, r} (λ_{G}) = P_{G, r 0} + λ_{G} K_{G, r}^{P}

We assume that the active and reactive power demand increase at a fixed power factor at the same time. In Eqs. (24a), (24b), (24c), (24d), λ_e, λ_g, and λ_G are the growth parameters of the power system, natural gas system, and generators, respectively; P_l,_i₀ and Q_l,_i₀ are the initial active power load and initial reactive power load of bus i in the power system, respectively; L_l,_m₀ is the initial natural gas load at node m in the natural gas system; P_G,_r₀ is the initial active power generation of the rth NGU; while

K_{l, i}^{P}

K_{l, i}^{Q}

K_{l, m}^{L}

, and

K_{G, r}^{P}

are the active power load increment, reactive power load increment, gas load increment, and active power generation increment of the IEGS in one step length, respectively. The MEF calculation model of Eqs. (10), (11), (12) is further modified in the CMEF, as follows:

(25a)

H_{eg} (x_{eg}, y_{eg} (λ)) = 0

(25b)

x_{eg} = [V, θ, p]

(25c)

y_{eg} (λ) = [P^{sp}, Q^{sp}, L^{sp}]

where function H_eg() is the modified MEF equations.

To simplify the analysis, there are three modes that can be chosen to simulate the load growth in the system, as shown in Ref. [20]. In this paper, the mode in which the active power, reactive power, and/or gas flow of all load buses (or nodes) increase is selected as the load growth mode, and a fixed load growth direction K^L is given, as follows:

(26a)

K^{L} = [K_{1} / K_{\sum}, K_{2} / K_{\sum}, . . ., K_{n_{t}} / K_{\sum}]

(26b)

K_{\sum} = \sum_{i = 1}^{n_{t}} K_{i}

where n_t is the number of buses (or nodes) in the IEGS. To unify the form of the equations, K_i in Eq. (26b) represents

K_{l, i}^{P}

K_{l, i}^{Q}

, and

K_{l, m}^{L}

in active power load growth, reactive power load growth, and gas load growth, respectively, according to Eqs. (24a), (24b), (24c); K_Ʃ is the sum of K_i.

3.1.2. Prediction and correction processes

Prediction and correction processes are used to obtain the next-step MEF solution x_eg ⁽^k⁺¹⁾ based on the current solution x_eg ⁽^k⁾. In the prediction process, we set the step-k MEF solution x_eg ⁽^k⁾ as depicted in Eq. (27).

(27)

[x_{eg}^{(k)}, λ^{(k)}] = [V^{(k)}, θ^{(k)}, p^{(k)}, λ^{(k)}]

In the next step (named step-(k + 1)), the solution x_eg⁽^k⁺¹⁾ is calculated through the tangent method. First, the tangent vector of x_eg⁽^k⁾ is determined as the prediction direction for x_eg⁽^k⁺¹⁾. The tangent vector τ is expressed by Eq. (28).

(28)

τ = [\begin{matrix} d x_{eg} \\ d λ \end{matrix}] = {[\begin{matrix} \begin{matrix} \frac{\partial H_{eg}}{\partial x_{eg}} & \frac{\partial H_{eg}}{\partial λ} \end{matrix} \\ e_{h} \end{matrix}]}^{- 1} [\begin{matrix} 0 \\ 1 \end{matrix}]

where e_h is a row vector. In e_h, the hth number is set to 1 and the other number is 0 to ensure that the vector τ has a certain solution (because the number of variables to be determined is one greater than the number of MEF functions). In this paper, h is selected as the variable that changes fastest in the tangent direction of dV_i/V_i. The selected principle of h is expressed by Eq. (29).

(29)

|\frac{d V_{h}}{V_{h}}| = \max \{|\frac{d V_{1}}{V_{1}}|, |\frac{d V_{2}}{V_{2}}|, . . ., |\frac{d V_{n_{PQ}}}{V_{n_{PQ}}}|\}

where n_PQ is the number of PQ buses (the active power and reactive power injections of the bus are given) in the power system.

Thus, the prediction value of the step-(k + 1) solution x_eg ⁽^k^+1)* is calculated by Eq. (30) with a fixed step length σ.

(30)

[x_{eg}^{(k + 1) *}, λ^{(k + 1) *}] = [x_{eg}^{(k)}, λ^{(k)}] + σ τ

Then, in the correction process, to ensure that the MEF calculation has a certain solution, the CMEF calculation model is extended, as depicted in Eq. (31).

(31)

\{\begin{matrix} H_{eg} (x_{eg}) = 0 \\ V_{h} - V_{h}^{*} = 0 \end{matrix}

The prediction value x _eg ⁽ ^k ^+1)* is set to be the initial value of Eq. (31), and the solution x _eg ⁽ ^k ⁺¹⁾ is then obtained by the decoupling-solving MEF calculation method, as described in Section 2.1.4.

Next, we judge whether the result satisfies the IEGS-SSR boundary criteria in Eq. (22). If it does, we continue the next prediction and correction process until at least one of the static stability criteria of the voltage or gas pressure is not satisfied. At that time, the boundary point of the IEGS-SSR in a specific NGU scheduling direction is finally determined.

3.2. Sampling of possible NGU scheduling directions for hypersurface fitting

To obtain the whole IEGS-SSR, all the possible growth directions of NGUs in the n-dimensional space should be determined. However, the number of possible growth directions is theoretically infinite. To simplify the calculation, a multidimensional hyperplane sampling method is proposed in order to sample several possible NGU scheduling directions evenly and then find the exact IEGS-SSR boundary points on these sampling directions. These boundary points will then be used in hypersurface fitting in the next section.

The active power injection vector ω_eg of the r_s chosen NGUs in the system is given by Eq. (32).

(32)

ω_{eg} = [P_{G,1}, P_{G,2}, . . ., P_{G, r_{s}}]

According to Eq. (24d), with load growth, vector ω_eg can be further described as follows:

(33a)

ω_{eg} = ω_{eg,0} + λ_{G} K_{G, Σ}^{P} ω_{eg,p}

(33b)

K_{G, Σ}^{P} = \sum_{r = 1}^{r_{s}} K_{G, r}^{P}

(33c)

ω_{eg,p} = [α_{1}^{P}, α_{2}^{P}, \cdot \cdot \cdot, α_{r_{s}}^{P}]

(33d)

α_{r}^{P} = K_{G, r}^{P} / K_{G, Σ}^{P}, r = 1, 2, \cdot \cdot \cdot, r_{s}

where ω_eg,0 is the constant vector of the NGU initial power generation;

K_{G, Σ}^{P}

is the sum of active power generation of all chosen NGUs;

α_{r}^{P}

is the power generation percentage of NGU_r. Thus, ω_eg,p can be seen as an NGU scheduling direction vector. When enough NGU scheduling direction vectors are sampled, the boundary points on these directions can be determined by means of the CMEF method.

An IEGS with three NGUs is taken as an example in this section to show how to sample enough directions in three-dimensional (3D) space.

As shown in Fig. 4(a), the basis vectors of the 3D space are

\vec{OA}

\vec{OB}

, and

\vec{OC}

, and points A, B, and C form a triangle plane. For every point D_k in plane ABC, the vector

\vec{{OD}_{k}}

can be a possible scheduling direction vector. A ternary diagram is carried out to help visualize D_k in plane ABC, as shown in Fig. 4(b). We divide the three edges of the triangle into T sections (Fig. 4(b) shows T = 6). Through the divided points, we draw the parallel lines of the other two edges. There are a total of

C_{T + 2}^{2}

intersection points (including divided points and vertexes), and the coordinate of every intersection point in the ternary diagram can be seen as a sampling value of ω_eg,p.

To prove the uniformity of the sampling direction sampling, two uniformity indexes, ratio and var, are used in this paper [53], as shown in Eqs. (34a), (34b), (34c). Here, d_k refers to the distance between sampling point D_k and its closest sampling point; n_s is the number of the NGU growth directions sampled;

\bar{d}

is the average value of all d_k.

(34a)

ratio = \frac{\min d_{k}}{\max d_{k}}

(34b)

var = \frac{1}{n_{s}} \sum_{k = 1}^{n_{s}} {(d_{k} - \bar{d})}^{2}

(34c)

\bar{d} = \frac{1}{n_{s}} \sum_{k = 1}^{n_{s}} d_{k}

Three methods are then selected to compare the uniformity of the sampling points: ① the multidimensional hyperplane sampling method used in this paper, ② the Monte Carlo sampling method, and ③ the space segmentation algorithm method.

As shown in Fig. 5, the black points are the sampling points. In Fig. 5, 120 points are sampled using the multidimensional hyperplane sampling method (T = 14) and the Monte Carlo sampling method, and 124 points are sampled using the space segmentation algorithm method. Visually, the sampling uniformity of the multidimensional hyperplane sampling method is better than the others. Table 2 shows the specific index results of the three methods.

The result shows that the uniformity of the multidimensional hyperplane sampling method is the best of the three commonly used methods. This method can create sampling points evenly in the sample space and further avoid local overfitting in hypersurface fitting.

3.3. Hypersurface fitting of IEGS-SSR in 3D space

The IEGS-SSR in a 3D active power injection space is a 3D hypersurface. In this paper, the hypersurface is approximated through a fitting method using the sampling boundary points above in 3D space based on the Delaunay triangulation method.

In 3D space, every three points that are not collinear can form a unique plane (or hyperplane), as plane ABC shows in Fig. 6. The vertexes A, B, and C are the exact boundary points calculated through the CMEF method in three different directions (ignoring the systematic error caused by the increment of one step). The triangular plane ABC is fitted to approximate the exact IEGS-SSR within a certain area.

Using the Delaunay triangulation method, for a group of sampling points with the growth directions sampled as shown in Fig. 4(b), each point forms several triangular planes together with the adjacent points. The final IEGS-SSR is fitted as a hypersurface patched through all these triangular planes piecewise, like the surface of a football (or soccer ball in the United States).

3.4. Flowchart of IEGS-SSR determination

A flowchart of the IEGS-SSR determination is shown in Fig. 7. First, the initial operation point ω _eg,0 is obtained via MEF calculation. Practically, the direction of the load (both natural gas load and power load) growth can be predicted by the load forecasting method. For an IEGS containing R NGUs, the IEGS operator selects r _s (r _s ≤ R) NGUs to meet the load growth demand. Therefore, in the r _s-dimensional NGU active power injection space, n _s growth directions based on the multidimensional hyperplane sampling method are generated. Each of the n _s NGU growth directions is chosen in turn, and the boundary points are obtained through CMEF. Finally, n _s boundary points are obtained to generate the boundary of the IEGS-SSR through hypersurface fitting.

3.5. Fitting error analysis

To demonstrate the accuracy of the hypersurface fitting method, the fitting error e _f is obtained for further study. As shown in Fig. 6, point E is the boundary point on the fitted plane ABC and point F is the boundary point on the exact IEGS-SSR on the same growth direction of the NGUs. We calculate the fitting error with Eq. (35).

(35)

e_{f} = \frac{∥\vec{O^{'} F}∥ - ∥\vec{O^{'} E}∥}{∥\vec{O^{'} F}∥}

Two types of fitting triangular planes are generated through the sampling boundary points: In one type, the vertexes of the triangular planes are limited by the same kind of boundary criteria only, as shown in Fig. 8(a); in the other type, there are two types of vertexes, limited by the two different boundary criteria of the power system and natural gas system, as shown in Fig. 8(b).

The yellow planes are the final fitting planes of the two types. The SGSR and SVSR planes intersect in Type II, whereas they do not intersect in Type I. In Fig. 8(a), the plane

D^{'} E^{'} F^{'}

is “above” the other plane

A^{'} B^{'} C^{'}

. Thus, the fitting error of this type mainly depends on the accuracy of the boundary point calculation and the above hypersurface fitting method. In Fig. 8(b), there is an intersecting line, considering the intersection of the SGSR and SVSR in Type II. The points on the intersecting line are the boundary points limited by both the boundary criteria of the power system and the natural gas system, such as points M and N (where M and N are in the planes AOB and AOC, respectively). In this case, the maximum fitting error occurs around the intersecting line MN. In fact, the hypersurface can be divided into two parts: One part is the plane determined by the DMN and AMN, while the other is the plane determined by the MNEF and MNBC. These two parts can be seen as two Type-I planes, and the corresponding fitting planes will be planes AMN and MNEF. It is obvious that there is still a margin between the fitting plane AEF and the hypersurface AMN-MNEF for NGU dispatching. In other words, compared with the hypersurface AMN-MNEF, the fitting plane AEF is a more conservative fitting result for the IEGS-SSR. Thus, the fitting error of the plane AEF to the exact IEGS-SSR is larger, especially around MN. More boundary points should be sampled to reduce the fitting error in the planes as Type II.

4. Case studies

The proposed method is verified by a typical IEGS with five NGUs. The IEGS contains the New England (USA) 39-bus power systems [54] and a 22-node natural gas systems [49]. Topology diagrams of the power system and the natural gas system are introduced first; then, a case with four scenarios is studied to demonstrate the effectiveness of the proposed method.

4.1. Topology diagram of the IEGS and settings

A topology diagram of the IEGS is shown in Fig. 9. E_i stands for the ith bus of the power system, and G_m stands for the mth node of the natural gas system. There are ten generators, and bus 31 is set to be the slack bus of the power system. The initial gas pressure is 58 bar at gas source node G ₈. An IEGS with five NGUs is presented in a case study to form a region in a 3D active power injection space. The power factor of power load growth is assumed to be a constant at 0.8; the constant vector of the NGU initial power generation ω _eg,0 = [650,632,508,650,250], and

K_{G, Σ}^{P}

is set as 10 MW.

In this case, the IEGS-SSR is defined in the power injection space of three NGUs. The following four scenarios are studied to show the relationship between the IEGS-SSR, SVSR, and SGSR.

Scenario I: NGU₁, NGU₂, and NGU₃ (named NGU group O) are chosen to meet the power demand. Then, the SVSR and SGSR are formed separately, ignoring the system coupling. The load growth in the power system is converted to the gas consumption growth of the nodes connected to the chosen NGUs.

Scenario II: The IEGS-SSR is formed considering both the power system and the natural gas system boundary criteria.

Scenario III: The IEGS-SSR is formed considering both the power system and the natural gas system boundary criteria, and the initial natural gas loads of all gas nodes increase by 5%, 15%, 30%, and 50%.

Scenario IV: We form the new IEGS-SSRs in 3D active power injection spaces of other combinations of NGUs that meeting the power demand.

4.2. Case studies

In this section, the red and blue points in the figures represent the sampling boundary points constrained by the gas pressure or voltage. The black lines show the boundary of the fitted hyperplane.

4.2.1. Scenario I

The SGSR and SVSR of the IEGS are shown in Fig. 10. When T = 10, a total of 66 points are sampled for hypersurface fitting.

4.2.2. Scenario II

Taking both gas pressure and voltage restrictions into account, the IEGS-SSR is formed as follows (Fig. 11). When focusing on the IEGS-SSR results of Scenarios I and II, it is clear that the IEGS-SSR is the intersection of the SGSR and SVSR. Compared with the SVSR, the IEGS-SSR shrinks after natural gas system static gas pressure stability constraints are considered. It is evident that several operating points that were previously considered to be static voltage stable within the SVSR now exhibit static gas pressure instability and fall outside the boundary of the IEGS-SSR. This result indicates that the SVSR tends to be overly optimistic when characterizing the IEGS-SSR within an IEGS.

Here, we take two NGU growth directions as examples. Table 3 provides detailed information on the two directions.

As shown in Fig. 12(a), the gas pressure of G ₂₁ drops to zero when the power generation of NGU₂ increases to 1802 MW (where λ _G = 195), and the natural gas system gas pressure limits the further growth of the load in the power system. The IEGS static stability margin is much lower than the static voltage stability margin of the SVSR in direction A (where λ _G = 254). In Fig. 12(b), the voltages in the power system drop, and the Jacobian matrix is finally singular when the power generation of NGU₁ increases to 2222 MW (where λ _G = 262). This result means that the key constraints restricting the MEF are the constraints of the power system.

Furthermore, to demonstrate the inclusion relation between the IEGS-SSR and the IEGS security region, a comparison of the two regions is presented in Fig. 13. There is a significant margin between the boundaries of the IEGS-SSR and the IEGS security region. This result confirms that the IEGS security region is encompassed by the IEGS-SSR.

4.2.3. Scenario III

The results of Scenario III are illustrated in Fig. 14.

Comparing Fig. 14, Fig. 11, it shows the impact of the gas load growth in a natural gas system on the IEGS-SSR. The load growth in the natural gas system leads to a decrease in the IEGS static stability margin. The results of ω _eg,p = [0.7,0.3,0] are provided as an example. The active power generation of the three NGUs is [2288,1334,508] (where λ _G = 234) in the SVSR of Scenario I, which represents the boundary point considering only the power system static voltage stability constraints. The results of the growth parameter λ _G in the SGSR and IEGS-SSR with different gas load increases for the same ω _eg,p are shown in Table 4.

With natural gas load growth, the growth parameter λ _G in the SGSR decreases from 247 to 201, which is ultimately lower than that in the SVSR. The growth parameter λ _G in the IEGS-SSR begins to decrease when the gas load in the natural gas system increases by 15%. This result means that, in the IEGS-SSR, as the intersection of the SGSR and SVSR, the key constraints restricting the IEGS-SSR shift from the static voltage stability constraints in the power system to the static gas pressure stability constraints in the natural gas system in these growth directions. The gas load growth in a natural gas system will reduce the power generation capability of the NGUs, which will lead to a decrease in the IEGS static stability margin. With the ever-increasing natural gas loads in the natural gas system, the IEGS-SSR shrinks further.

4.2.4. Scenario IV

In the above case with five NGUs coupling the power system and the natural gas system, there are ten combinations of NGUs (in 3D space). In this scenario, six other groups of NGUs are selected to meet the power growth demand, as follows:

Group A: NGU₁, NGU₂, and NGU₄

Group B: NGU₁, NGU₂, and NGU₅

Group C: NGU₁, NGU₄, and NGU₅

Group D: NGU₁, NGU₃, and NGU₅

Group E: NGU₂, NGU₃, and NGU₄

Group F: NGU₂, NGU₄, and NGU₅

Fig. 15 shows the different IEGS-SSRs formed in different NGU-combination spaces.

Comparing the results in Fig. 11, Fig. 15, it reveals that, even in the same power-generation growth direction, the IEGS load margin varies when different NGUs are chosen to meet the power growth demand. Table 5 compares the detailed boundary information of NGU group O, group A, and group B in two NGU scheduling directions.

When NGU₁ and NGU₂ are kept unchanged and the third NGU is varied to meet the load demand in the IEGS, the growth parameter λ _G in the IEGS-SSR varies significantly, as shown in Table 5. In direction C, the λ _G of group O is larger than that of group B, while the situation is reversed in direction D. Thus, when dispatching, it is better to find the best groups of NGUs with the largest load margin to keep the IEGS in a secure operating state.

4.3. Computational accuracy and time analyses

To check the fitting error of the IEGS-SSR hypersurface, a pointwise method based on Monte Carlo sampling is carried out, sampling more than 20 000 boundary points to form the IEGS-SSR (considered the exact boundary). Then, we calculate the fitting error using Eq. (35). The fitting error ternary diagrams of groups O, B, and C are shown in Fig. 16. Detailed fitting error results are provided in Table 6.

As shown in Fig. 16, in most areas of the three IEGS-SSRs, the average positive fitting errors are within 0.32% of the exact boundary. High positive fitting errors mainly appear within the fitting triangular planes with different colored vertexes. This result is consistent with the discussion in Section 3.5. The fitting errors reveal that the IEGS-SSR can be accurately determined. However, the presence of negative fitting errors indicates that the fitted IEGS-SSR may slightly overestimate the evaluation of static stability in an IEGS. In practice, power systems in North America are typically required to maintain a static voltage stability margin of over 5% under normal conditions [55]. This criterion can also be extended to guide the operation of the IEGS. Even though the IEGS-SSR displays optimism with a minimum negative fitting error of approximately −0.31%, this should not significantly impact the operation and dispatch processes within the IEGS.

In a power system, the security region is used to assess the system’s security level, typically at intervals of half an hour or even shorter [56]. Similarly, the IEGS-SSR can be calculated based on the current operating point and used to evaluate the static stability of an IEGS within similar half-hour intervals. Therefore, computational efficiency is crucial in determining the IEGS-SSR. The computation time of the IEGS-SSR in Scenarios II and IV is shown in the Table 7.

The computation time of the IEGS-SSR is related to the number of sampled boundary points n _s in the active power injection space of the chosen NGUs. Thus, the calculation speed increases by at least 90% compared with that of the pointwise method, considering that the number of sampled boundary points decreases from 20 000 to 66. The results show that the proposed method can effectively form the IEGS-SSR.

5. Conclusions

In this paper, the concept of the power system SVSR is extended to IEGSs. First, an MEF model of an IEGS is presented. Then, a static stability analysis of the IEGS is proposed, and the IEGS-SSR and its boundary are defined based on the analysis. By tracing the boundary of the IEGS-SSR via the CMEF and the multidimensional hyperplane sampling method, the 3D IEGS-SSR is obtained. The conclusions are summarized as follows:

(1) The IEGS-SSR gives a global perspective of the system static stability margin. In view of the stability and security of the NGU energy supply, it is necessary to consider the constraints of the natural gas system in the maximum transmission power determination of the power system.

(2) The IEGS-SSR shrinks when the natural gas system is considered. When the gas load increases in the natural gas system, the IEGS-SSR shrinks further.

(3) The IEGS-SSRs vary according to the multidimensional spaces of NGUs chosen to meet the load growth. To ensure that the IEGS has enough stability and security margins when dispatching, more regions should be formed to guide the operators to find the best group of NGUs and the corresponding scheduling direction.

(4) The proposed method can accurately form the IEGS-SSR hypersurface with an average positive fitting error within 0.32%.

The IEGS-SSR provides an online security monitoring method for IEGS operators; it simplifies the security monitoring process, allowing operators to simply confirm whether the current operating point is within the IEGS-SSR in order to estimate the static stability of the IEGS. Moreover, the IEGS-SSR can be integrated as boundary constraints within the IEGS dispatch and system control. Further research will focus on the topological properties of the IEGS-SSR boundary and a fast IEGS-SSR boundary-determination method to determine the IEGS-SSR with high efficiency and to make the IEGS-SSR widely applicable in practice.

Acknowledgments

This work was financially funded by the National Natural Science Foundation of China (52222704 and 52177107).

Compliance with ethics guidelines

Yunfei Mu, Zhibin Liu, Xiangwei Guo, Hongjie Jia, Kai Hou, Xiaodan Yu, Bofeng Luo, and Hairun Li declare that they have no conflicts of interest or financial conflicts to disclose.

References

Publishing order | Descend order by publishing year | Descend order by cited within

[1]	Paul S, Dey T, Saha P, Dey S, Sen R. Review on the development scenario of renewable energy in different country. In: Proceedings of 2021 Innovations in Energy Management and Renewable Resources; 2021 Feb 5-7; Kolkata, India. IEEE; 2021. p. 1-2.

[2]	Youssef E, El Azab RM, Amin AM. Comparative study of voltage stability analysis for renewable energy grid-connected systems using PSS/E. In: Proceedings of SoutheastCon 2015; 2015 Apr 9-12; Fort Lauderdale, FL, USA. IEEE; 2015. p. 1-6.

[3]	J. Zhang, Y. Mu, Z. Wu, Z. Liu, Y. Gao, H. Jia, et al.. Optimal scheduling method of regenerative electric heating for emergency residential building heating: an affine arithmetic-based model predictive control approach. IET Energy Syst Integr, 5 (1) ( 2023), pp. 40-53

[4]	Y. Mi, C. Liu, J. Yang, H. Zhang, Q. Wu. Low-carbon generation expansion planning considering uncertainty of renewable energy at multi-time scales. Glob Energy Interconnect, 4 (3) ( 2021), pp. 261-272

[5]	D. Chen, Z. Bao, L. Wu. Integrated coordination scheduling framework of electricity-natural gas systems considering electricity transmission N - 1 contingencies and gas dynamics. J Mod Power Syst Clean Energy, 7 (6) ( 2019), pp. 1422-1433

[6]	V. Shabazbegian, H. Ameli, M.T. Ameli, G. Strbac.Stochastic optimization model for coordinated operation of natural gas and electricity networks. Comput Chem Eng, 142 ( 2020), p. 107060

[7]	B.C. Erdener, K.A. Pambour, R.B. Lavin, B. Dengiz. An integrated simulation model for analyzing electricity and gas systems. Int J Electr Power Energy Syst, 61 ( 2014), pp. 410-420

[8]	S. Badakhshan, M. Ehsan, M. Shahidehpour, N. Hajibandeh, M. Shafie-Khah, J.P.S. Catalão. Security-constrained unit commitment with natural gas pipeline transient constraints. IEEE Trans Smart Grid, 11 (1) ( 2020), pp. 118-128

[9]	Hui H, Ding Y, Luan K, Xu D. Analysis of “8.15” blackout in Taiwan and the improvement method of contingency reserve capacity through direct load control. In:Proceedings of 2018 IEEE Power & Energy Society General Meeting (PESGM); 2018 Aug 5-10; Portland, OR, USA. IEEE; 2018. p. 1-5.

[10]	P. Fairley. The troubled link between gas and electricity grids. IEEE Spectrum, 53 (6) ( 2016), pp. 11-12

[11]	J. Bialek.What does the GB power outage on 9 August 2019 tell us about the current state of decarbonized power systems?. Energy Policy, 146 ( 2020), Article 111821

[12]	F.J. Mercede, J.C. Chow, H.H. Yan, R. Fischl. A framework to predict voltage collapse in power systems. IEEE Trans Power Syst, 3 (4) ( 1988), pp. 1807-1813

[13]

Abed AM.WSCC voltage stability criteria, undervoltage load shedding strategy, and reactive power reserve monitoring methodology. In:Proceedings of 1999 IEEE Power Engineering Society Summer Meeting. Conference Proceedings (Cat. No.99CH36364); 1999 Jul 18- 22; Edmonton, AB, Canada. IEEE; 1999. p. 191-7.

[14]	H. Kwatny, A. Pasrija, L. Bahar. Static bifurcations in electric power networks: loss of steady-state stability and voltage collapse. IEEE Trans Circ Syst, 33 (10) ( 1986), pp. 981-991

[15]	V.A. Venikov, V.A. Stroev, V.I. Idelchick, V.I. Tarasov. Estimation of electrical power system steady-state stability in load flow calculations. IEEE Trans Power Apparatus Syst, 94 (3) ( 1975), pp. 1034-1041

[16]	K. Alzaareer, M. Saad, H. Mehrjerdi, Q. Salem, S. Harasis, A.M.I. Aldaoudeyeh, et al.. Sensitivity analysis for voltage stability considering voltage dependent characteristics of loads and DGs. IEEE Access, 9 ( 2021), pp. 156437-156450

[17]	B. Liu, Y. Wang, Z. Yang, Y.W. Li, H.Z. Xie, Y. Zhang, et al.. Evaluation method of voltage stability based on minimum singular value. Appl Mech Mater, 511-512 ( 2014), pp. 1128-1132

[18]	M.J. Vahid-Pakdel, H. Seyedi, B. Mohammadi-Ivatloo. Enhancement of power system voltage stability in multi-carrier energy systems. Int J Electr Power Energy Syst, 99 ( 2018), pp. 344-354

[19]	Zhao B, Luo B, Zhou J, Meng X, Li Z.Static sensitivity analysis method of integrated energy system based on power flow model. In:Proceedings of 2018 2nd IEEE Conference on Energy Internet and Energy System Integration (EI2); 2018 Oct 20- 22; Beijing, China. IEEE; 2018. p. 1-9.

[20]	Qiao Z, Guo Q, Sun H, Jia L, Zhao W.Static voltage stability margin considering the coupling of natural gas and power system. In:Proceedings of 2018 IEEE Power & Energy Society General Meeting (PESGM); 2018 Aug 5-10; Portland, OR, USA. IEEE; 2018. p. 1-5.

[21]	H.D. Chiang, Z.Y. Wang, L. Zeng. Dynamic relationship between KKT saddle solutions and optimal solutions in AC OPF problems. IEEE Trans Power Syst, 39 (1) ( 2024), pp. 1637-1646

[22]	Wei W, Jia H, Zhang P, Wang C, Wu J, Lee ST.Development of power system voltage stability region (PSVSR) for static voltage security assessment. In:Proceedings of 2006 International Conference on Power System Technology; 2006 Oct 22- 26; Chongqing, China. IEEE; 2006. p. 1-6.

[23]	Y. Mu, H. Jia. An approach to determining the local boundaries of voltage stability region with wind farms in power injection space. Sci China Technol Sci, 53 (12) ( 2010), pp. 3232-3240

[24]	Yu Y, Li H, Han Q, Lee ST, Zhang P.Practical boundary of static voltage stability region in cut-set complex power space of power systems. In:Proceedings of 2005 IEEE/PES Transmission & Distribution Conference & Exposition: Asia and Pacific; 2005 Aug 18- 18; Dalian, China. IEEE; 2005. p. 1-6.

[25]	Y. Qiu, H. Wu, Y. Song, J. Wang. Global approximation of static voltage stability region boundaries considering generator reactive power limits. IEEE Trans Power Syst, 33 (5) ( 2018), pp. 5682-5691

[26]	H.D. Chiang, A.J. Flueck, K.S. Shah, N. Balu. CPFLOW: a practical tool for tracing power system steady-state stationary behavior due to load and generation variations. IEEE Trans Power Syst, 10 (2) ( 1995), pp. 623-634

[27]	E. Kuroda, M. Watanabe, D. Kato, N. Saito, M. Yatsu. Fast computation method of static voltage stability using geometric parameter adjustment for the continuation power flow. Electr Eng Jpn, 214 (2) ( 2021), pp. 1-10

[28]	W. Zhang, T. Wang, H.D. Chiang. A novel FFHE-inspired method for large power system static stability computation. IEEE Trans Power Syst, 37 (1) ( 2022), pp. 726-737

[29]	W.H. Tang, W.W. Li, J.H. Zheng, C.Q. Wu, L.X. Wang, Q.L. Wei, et al.. A composite voltage stability index for integrated energy systems based on L-index and the minimum eigenvalue of reduced Jacobian matrix. Int J Electr Power Energy Syst, 141 ( 2022), Article 108136

[30]	Liu Z, Mu Y, Guo X, Jia H, Wu G, Zu G.Static stability analysis of integrated electricity and gas system considering voltage and pressure. In:Proceedings of 2021 IEEE 5th Conference on Energy Internet and Energy System Integration (EI2); 2021 Oct 22- 24; Taiyuan, China. IEEE; 2021, p. 1750-4.

[31]	W. Xu, H. Lu, G. He, Y. Chen, D. Yang. Optimal schedule strategy of integrated energy system considering static voltage stability. Lect Notes Electr Eng, 718 ( 2021), pp. 660-675

[32]	S. Lu, Y. Xu, W. Gu, X. Fang, Z. Dong. On thermal dynamics embedded static voltage stability margin. IEEE Trans Power Syst, 38 (3) ( 2023), pp. 2982-2985

[33]	Tang Y. Static analysis of voltage stability. In: Voltage Stability Analysis of Power System. Singapore: Springer Singapore; 2021.

[34]	B. Rajani, P.S. Raju. Comparision of facts controllers for improvement of voltage/line stability in transmission system using SSSC & STATCOM. J Theor Appl Inf Technol, 41 ( 2012), pp. 192-200

[35]	B. Chen, Q. Guo, G. Yin, B. Wang, Z. Pan, Y. Chen, et al.. Energy-circuit-based integrated energy management system: theory, implementation, and application. Proc IEEE, 110 (12) ( 2022), pp. 1897-1926

[36]	A. Martinez-Mares, C.R. Fuerte-Esquivel. A unified gas and power flow analysis in natural gas and electricity coupled networks. IEEE Trans Power Syst, 27 (4) ( 2012), pp. 2156-2166

[37]	Yuan K, Sun C, Song Y, Li H.A multi-energy flow calculation method considering multiple energy coupling operation modes. In:Proceedings of 2020 IEEE 4th Conference on Energy Internet and Energy System Integration (EI2); 2020 Oct 30-Nov 1; Wuhan, China. IEEE; 2020. p. 2661-6.

[38]	N.J. Labiba. Fluid flow analysis of a transmission line of Jalalabad gas transmission and distribution system limited. Int J Oil Gas Coal Eng, 3 (2) ( 2015), pp. 24-32

[39]	Y.Q. Bao, M. Wu, X. Zhou, X. Tang. Piecewise linear approximation of gas flow function for the optimization of integrated electricity and natural gas system. IEEE Access, 7 ( 2019), pp. 91819-91826

[40]	T. Jiang, R. Zhang, X. Li, H. Chen, G. Li.Integrated energy system security region: concepts, methods, and implementations. Appl Energy, 283 ( 2021), p. 116124

[41]	Chen S, Wei Z, Sun G, Sun Y, Lu N. Steady-state security regions of electricity-gas integrated energy systems. In: Proceedings of 2016 IEEE Power and Energy Society General Meeting (PESGM); 2016 Jul 17-21; Boston, MA, USA. IEEE; 2016. p. 1-5.

[42]	Li S, Xia W, Liu Y.A sensitivity-based exchange model between L index and eigenvalues to quantify power system static voltage stability. In:Proceedings of 2023 IEEE 6th International Electrical and Energy Conference (CIEEC); 2023 May 12- 14; Hefei, China. IEEE; 2023. p. 397-401.

[43]	X. Li, G. Tian, Q. Shi, T. Jiang, F. Li, H. Jia.Security region of natural gas network in electricity-gas integrated energy system. Int J Electr Power Energy Syst, 117 ( 2020), p. 105601

[44]	Y.J. Qin, J.H. Zheng, Z. Li, Q.H. Wu.Reduction of non-linear many objectives for coordinated operation of integrated energy systems. Int J Electr Power Energy Syst, 117 ( 2020), p. 105657

[45]	M. Varadarajan, K.S. Swarup. Solving multi-objective optimal power flow using differential evolution. IET Gener Transm Distrib Inst Eng Technol, 2 (5) ( 2008), pp. 720-730

[46]	B. Gao, G.K. Morison, P. Kundur. Voltage stability evaluation using modal analysis. IEEE Trans Power Syst, 7 (4) ( 1992), pp. 1529-1542

[47]

Li P, Li Q, Bai J, Wang S, Huang Y.Accommodation ability evaluation of high renewable energy penetrated power system with large-scale concentrating solar power. In:Proceedings of 2021 6th Asia Conference on Power and Electrical Engineering (ACPEE); 2021 Apr 8- 11; Chongqing, China. IEEE; 2021, p. 263-8.

[48]	Que L, Jiang Z, Lai H, Xie Z, Li S.A reliability evaluation method for transmission grid based on risk theory. In:Proceedings of Third International Conference on Mechanical Design and Simulation (MDS 2023); 2023 Mar 3- 5; Xi’an, China. 2023.

[49]	Osiadacz AJ.Simulation and analysis of gas networks. Oak Ridge: US Department of Energy; 1987.

[50]	J. Xiao, Y. Qu, B. She, L. Lv, C. Song, Z. Sun, et al.. Complete limits of flow network based on critical flow: concept, model, algorithm, visualization, and applications. Arab J Sci Eng, 48 (1) ( 2023), pp. 985-1007

[51]	Y. Qiu, H. Wu, Y. Zhou, Y. Song. Global parametric polynomial approximation of static voltage stability region boundaries. IEEE Trans Power Syst, 32 (3) ( 2017), pp. 2362-2371

[52]	Hu L, Zhao J, Hongjie J, Guo X, Zeng Q.Computing the boundary of simple power system feasible region based on the optimization method. In: Proceedings of IEEE PES Innovative Smart Grid Technologies; 2012 May 21-24; Tianjin, China. IEEE; 2012. p. 1-6.

[53]	A. Bishnu, S. Desai, A. Ghosh, M. Goswami, S. Paul. Uniformity of point samples in metric spaces using gap ratio. SIAM J Discrete Math, 31 (3) ( 2017), pp. 2138-2171

[54]	Qiu C, Cheng J, Chen Q, Jiang P.An optimal method for load restoration in AC-DC system. In:Proceedings of 2017 China International Electrical and Energy Conference (CIEEC); 2017 Oct 25- 27; Beijing, China. IEEE; 2018. p. 107-11.

[55]	wecc. org. Internet. System performance criterion under normal conditions, following loss of a single BES element, and following extreme BES events. Report. Salt Lake City: Western Electricity Coordinating Council; 2012 Apr.

[56]	Y. Yu, Y. Liu, C. Qin, T. Yang. Theory and method of power system integrated security region irrelevant to operation states: an introduction. Engineering, 6 (7) ( 2020), pp. 754-777