Static voltage stability analysis of integrated smart energy systems

The analysis of multi-energy flows forms the cornerstone for the study of state estimation, safety assessment, and optimization in integrated smart energy systems (ISES). The interactions between various energy flows and the inherent variability of renewable energy sources often lead to significant challenges to the static stability of ISES. This paper investigates the static voltage stability of ISES under multi-energy coupling conditions through multi-energy flow analysis in electric-gas-thermal energy systems. First, a steady-state model of the ISES is constructed by representing the interconnected energy subsystems as equivalent sources and loads in the power grid. Subsequently, through the coupling elements of ISES, the power flows of the natural gas system (NGS) and the district heating system (DHS) are converted into active power in the electrical power system (EPS), resulting in an equivalent power flow equation for the ISES. Then, using Brouwer’s fixed-point theorem, the analytical sufficient conditions for solving the equivalent power flow equation are derived. Finally, a simulation model based on MATLAB/Simulink is established. The steady-state criterion for ISES is obtained in this paper, and the correctness and effectiveness of the proposed conclusions are verified by the simulation results.

1 Introduction

With the depletion of fossil fuels and the increase in renewable energy sources, Integrated Smart Energy Systems (ISES) have gained widespread attention (Li et al., 2020; Liu et al., 2022; Peng et al., 2021; Ma et al., 2024; Wang et al., 2024; Yang et al., 2020). ISES refers to the coordinated optimization and stable operation of various energy supply systems within a certain area through the integration of multiple energy sources such as electricity, heat, cooling, and gas. By leveraging advanced information processing technologies and innovative management models, ISES aims to form a multi-energy collaborative management system that can meet diverse energy needs, aligning with future trends of openness, low carbon, and sustainability in the energy sector. Compared to traditional energy systems, ISES breaks the independence of individual energy supply networks, considering the complementary and cascading utilization relationships between various heterogeneous energy sources, thereby promoting the upgrading of the energy system structure (Huang et al., 2022).

As the proportion of renewable energy increases and the need for energy interconnection grows, the scale and complexity of ISES are also continuously increasing. However, one of the ensuing problems is static voltage stability, which refers to the ability of the system to maintain voltage levels under steady-state conditions (Liu et al., 2020; Nguyen et al., 2018; Liu et al., 2018). As the main energy carrier, electrical power system (EPS)’s coupling relationship with NGS and DHS is becoming increasingly complex. ISES can become unstable through various mechanisms, such as the loss of long-term equilibrium due to the increasingly complex coupling relationships among EPS, NGS, and DHS. Particularly, the uncertainty in power flow injections can push the equilibrium point of ISES out of the feasible region, making the system’s steady-state equilibrium non-existent (Dvijotham et al., 2018). Therefore, accurate static voltage stability analysis of ISES is crucial for subsequent state estimation, safety analysis, and optimal control.

Currently, researchers worldwide have conducted extensive studies and analyses on ISES, focusing mainly on ISES modeling, power flow calculation, steady-state analysis, and multi-energy system security. For example, in Shabanpour-Haghighi and Seifi (2015), a steady-state model of an electric-gas-thermal coupling system is established, considering the valve point effect of gas turbines, part-load characteristics of combined heat and power (CHP) units, and the impact of gas boilers on the power flow distribution of the energy system. In Zhang et al. (2021), an effective probabilistic multi-energy flow (PMEF) calculation method is proposed to obtain the probabilistic information and energy flow distribution of hydrogen-injected ISES through probabilistic flow analysis. In Massrur et al. (2018), an adaptive discretization method is proposed to achieve satisfactory accuracy with low computational burden. Regarding multi-energy flow and thermal inertia (Huang et al., 2023), proposes a sequential simulation reliability assessment method. Considering the time characteristics of renewable energy (Wu and Wang, 2022), uses a sequential Monte Carlo simulation method to evaluate the reliability of ISES. In Sun et al. (2019), the droop control and equal incremental rate standard of ISES are considered.

In summary, existing studies on ISES tend to focus on the efficiency of multi-energy flow calculation and steady-state calculation, without the ability to determine whether a steady state exists in ISES quickly and accurately. Additionally, although static voltage stability has been widely studied in traditional power systems, the specific characteristics of ISES, particularly the coordinated operation and coupling of multiple energy sources, have not been fully addressed. This paper aims to analyze the static voltage stability of ISES under multi-energy flow coupling conditions. The main contributions can be summarized as follows: (1) Derivation of static voltage stability criteria for ISES under multi-energy flow coupling, ensuring the existence of equilibrium points in ISES. (2) Proposal of a multi-energy flow equation iteration algorithm based on the fixed-point method, including the construction of iterative expressions and the selection of initial values.

2 Mathematical model of ISES

2.1 Symbols and definitions

Definition 1. ${\mathbb{C}}^m$, ${\mathbb{C}}_{+}^m$ and ${\mathbb{C}}^{m\times m}$ are the complex $m$-dimensional vector, the positive complex $m$-dimensional vector and complex $m\times m$ matrix, respectively. Define $1_m={[111]}^{\text{T}}$, $0_m={[00\dots 0]}^{\text{T}}$, and $I_m$ is an m-dimensional unit matrix. Let $x={[x_1{x}_2{x}_m]}^{\text{T}}$, and define $[[x]]=diag\left\{x\right\}$. For $x\in {\mathbb{C}}^n$, $\bar{x}$: conjugate of $x\in {\mathbb{C}}^n$. Let $x\in {\mathbb{C}}^n$ be a vector with $x_i\ne 0$, then $\frac{1}{x}$ denotes the vector $\left[\begin{array}{cccc}\hfill \frac{1}{x_1}\hfill & \hfill \frac{1}{x_2}\hfill & \hfill \dots \hfill & \hfill \frac{1}{x_n}\hfill \end{array}\right]$.

Definition 2. Define ${\Vert x\Vert }_{\infty }=\underset{1\le i\le m}{\max }\left\{|a_{ij}|\right\}$, for $x\in {\mathbb{C}}^n$, ${\Vert A\Vert }_{\infty }=\underset{1\le i\le m}{\max }\left\{{\sum }_{j=1}^m|a_{ij}|\right\}$.

Lemma 1. ($Brouwe{r}^{\prime }sFixed-PointTheorem$ (Khalil, 2001)). Let $f:{\mathcal{U}}^n\mapsto {\mathcal{U}}^n$ be a continuous mapping and $D\subset {\mathcal{U}}^n$ be a compact convex set. Then, if $f(x)$ is a self-mapping (i.e., $f(x)\in D$ for any $x\in D$), then there is a $x^{*}\in D$ such that $f(x^{*})=x^{*}$.

Lemma 2. Let $M_n$ denote the set of all $n\times n$ matrices. If for all $A,B\in M_n$, the following two properties are satisfied (Molitierno, 2016): (1)$‖A‖\ge 0$ and $‖A‖=0$ if and only if $A=0$; (2)$‖AB‖\le ‖A‖‖B‖$, then the function $‖\cdot ‖$: $M_n\to \mathbb{C}$ is a matrix norm.

Based on the characteristics of power flow distribution in the EPS model, NGS model, and DHS model in ISES, the power flow equation of ISES is determined. First, the mathematical model of ISES is shown in Figure 1.

Figure 1

Figure 1. Schematic diagram of ISES structure.

2.2 EPS model

In the EPS model used in this paper, there are n + 1 buses, including l balanced bus, g distributed generation buses, and l load buses, corresponding to l balanced node, g generator nodes, and l load nodes, with l = n − g. The EPS is set as a voltage regulation source, and the phase angle of the balance bus is fixed. Without loss of generality, it is assumed that the balance bus has a voltage of $V_S=1.05\angle 0^{°}$.The distributed generation buses are numbered $G=\left\{\begin{array}{ccc}\hfill 1,\hfill & \hfill \dots \hfill & \hfill ,g\hfill \end{array}\right\}$, and the load buses are numbered $L=\left\{\begin{array}{ccc}\hfill g+1,\hfill & \hfill \dots \hfill & \hfill ,n\hfill \end{array}\right\}$,with l, n and g all being integers. The topology of the EPS is shown in Figure 2.

Figure 2

Figure 2. The topology of EPS.

For PQ buses, the injected power is $s_i=P_i+j{Q}_i$, $\forall i\in \left\{\begin{array}{cc}\hfill G,\hfill & \hfill L\hfill \end{array}\right\}$.Considering the classification of nodes, the admittance matrix is introduced as shown in Equation 1 (Wang et al., 2017).

$\left[\begin{array}{cc}\hfill Y_{SS}\hfill & \hfill Y_{SL}\hfill \\ \hfill Y_{LS}\hfill & \hfill Y_{LL}\hfill \end{array}\right]\cdot \left[\begin{array}{c}\hfill V_S\hfill \\ \hfill V_L\hfill \end{array}\right]=\left[\begin{array}{c}\hfill I_S\hfill \\ \hfill I_L\hfill \end{array}\right](1)$

where $I_S\in {\mathbb{C}}^1$ represents the current vector of the balance bus, $I=I_D+j{I}_Q\in {\mathbb{C}}^n$represents the current vectors of the generator buses and load buses, and $V\in {\mathbb{C}}^n$ is the vector of generator and load voltages. Therefore, the power flow equation of the EPS can be expressed as Equation 2.

$s={\displaystyle \sum _{i=1}^n}\overline{Y_i}{V}_i\overline{V_i},\forall i\in \left\{\begin{array}{cc}\hfill G,\hfill & \hfill L\hfill \end{array}\right\}(2)$

Further, Equation 2 can be written in a compact form, resulting in Equation 3.

$\left[\left[V\right]\right]\overline{YV}=-s(3)$

The reference direction of the current is opposite to the reference direction of the voltage, hence the negative sign on the right side of Equation 3.

2.3 NGS model

According to Ma et al. (2024), for an N-node NGS, the network equation describing the gas flow is given by Equation 4.

$Y_g{p}_g=G_g(4)$

where $Y_g$ is the N-order generalized node admittance matrix, representing the network structure and parameters;$p_g$ is the N-dimensional node pressure vector; $G_g$ is the N-dimensional generalized node injection vector. Similar to the boundary conditions given for power flow calculations in electric power systems, all nodes are divided into pressure-defined nodes and injection-defined nodes as shown in Table 1, thus providing N state variables to solve the remaining N state variables. Based on this classification, the network equation describing the gas flow for an N-node NGS is given by Equation 5.

$\left[\begin{array}{cc}\hfill Y_{g,gg}\hfill & \hfill Y_{g,gp}\hfill \\ \hfill Y_{g,pg}\hfill & \hfill Y_{g,pp}\hfill \end{array}\right]\cdot \left[\begin{array}{c}\hfill p_{g,g}\hfill \\ \hfill p_{g,p}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill G_{g,g}\hfill \\ \hfill G_{g,p}\hfill \end{array}\right](5)$

where $G_{g,g}$ is the injection of injection-defined nodes;$p_{g,p}$ is the pressure of pressure-defined nodes, both are known boundary conditions;$p_{g,g}$ is the pressure of injection-defined nodes; $G_{g,p}$ is the injection of pressure-defined nodes, both are unknown state variables; Y is the generalized node admittance matrix. According to Equation 4, the pressure of injection-defined nodes and the injection of pressure-defined nodes can be solved as shown in Equation 6.

$\left\{\begin{array}{c}\hfill p_{g,g}=Y_{g,gg}^{-1}\left(G_{g,g}-Y_{g,gp}{p}_{g,p}\right)\hfill \\ \hfill G_{g,p}=Y_{g,pg}{p}_{g,g}+Y_{g,pp}{p}_{g,p}\hfill \end{array}\right.(6)$

Table 1

Table 1. Classification of nodes in NGS.

When there are pressure-defined nodes in the natural gas network, the node admittance matrix $Y_{g,gg}$ corresponding to injection-defined nodes is non-singular, i.e.,$Y_{g,gg}\ne 0$.

2.4 DHS model

The steady-state calculation of DHS is divided into hydraulic steady-state calculation and thermal steady-state calculation, where the hydraulic steady-state power flow calculation is completely analogous to the steady-state power flow calculation of NGS.

2.4.1 Hydraulic steady-state model

According to Ma et al. (2024), as shown in Equation 7, this equation is consistent with the mathematical form of the power grid equation and is suitable for solving hydraulic flow in radial or looped networks.

$Y_h{p}_h=G_h(7)$

where $Y_h$ is the N-dimensional hydraulic node admittance matrix, $p_h$pis the N-dimensional node pressure vector, and $G_h$ is the N-dimensional node injection vector. Similar to NGS, for convenience of solution, Ma et al. (2024) divides the DHS nodes into injection-defined nodes and pressure-defined nodes as shown in Table 2.

Table 2

Table 2. Classification of nodes in DHS (hydraulic energy flow).

Based on this classification, the network equation describing the hydraulic is shown by Equation 8.

$\left[\begin{array}{cc}\hfill Y_{h,gg}\hfill & \hfill Y_{h,gp}\hfill \\ \hfill Y_{h,pg}\hfill & \hfill Y_{h,pp}\hfill \end{array}\right]\cdot \left[\begin{array}{c}\hfill p_{h,g}\hfill \\ \hfill p_{h,p}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill G_{h,g}\hfill \\ \hfill G_{h,p}\hfill \end{array}\right](8)$

where $G_{h,g}$ is the injection of injection-defined nodes; $p_{h,p}$ is the pressure of pressure-defined nodes; $p_{h,g}$ is the pressure of injection-defined nodes; $G_{h,p}$ is the injection of pressure-defined nodes; Y is the generalized node admittance matrix. Therefore, the hydraulic state variables of DHS shown by Equation 9 can be solved from Equation 7.

$\left\{\begin{array}{c}\hfill p_{h,g}=Y_{h,gg}^{-1}\left(G_{h,g}-Y_{h,gp}{p}_{h,p}\right)\hfill \\ \hfill G_{h,p}=Y_{h,pg}{p}_{h,g}+Y_{h,pp}{p}_{h,p}\hfill \end{array}\right.(9)$

When there are pressure-defined nodes in DHS, the node admittance matrix $Y_{h,gg}$ corresponding to injection-defined nodes is non-singular. Finally, the branch state variables are calculated from the node state variables and the hydraulic branch equation, as shown in Equation 10.

$G_{h,b}=y_h\left(G_h^T{p}_h-E_{h,b}\right)(10)$

where $G_{h,b}$ is the branch injection vector; $Y_{h,b}$ is the branch admittance matrix; $a_h$ is the node-branch incidence matrix; $G{E}_{h,b}$ is the column vector composed of water pressure source parameters for each branch.

2.4.2 Thermal steady-state model

According to Liu et al. (2022), the network equation for the thermal steady-state of DHS is given by Equation 11.

$Y_t{T}_{t,t}=h_t(11)$

where $Y_t$ includes the heat transfer coefficient matrix and the node-branch incidence matrix, representing the network structure and parameters; $h_t$ includes the node injection temperature and the temperature difference of the heat exchanger working fluid, representing the network operation boundary conditions; $T_{t,t}$ represents the state variables of the network, i.e., the terminal temperature of the branches. Since $Y_t$ is non-singular, the terminal temperature of the branches $T_{t,t}$ can be obtained from Equation 11. After obtaining the terminal temperature of the branches, the initial temperature of the branches can be determined according to the thermal branch equation, as shown in Equation 12.

$T_{t,f}=K_t^{-1}{T}_{t,t}(12)$

2.4.3 Combined water-thermal network model

It should be noted that in practical DHS, the thermal load power is usually known, and fluctuations in the thermal load simultaneously affect both the hydraulic and thermal circuits. According to the power flow calculation method of electric power systems, the network equation of DHS in terms of node thermal power is obtained, as shown in Equation 13.

${\varphi }_i=c_p\left(T_{si}-T_{oi}\right)\left({\displaystyle \sum _{j=1}^{n_1}}{y}_{\mathit{hij}}{A}_{\mathit{hij}}{p}_{\mathit{hij}}-{\displaystyle \sum _{j=1}^{n_1}}{y}_{bj}{E}_{bj}\right)(13)$

where ${\varphi }_i$ is the injection thermal power (MW) of node i, $T_{si}$ and $T_{oi}$ are the supply temperature and load output temperature (°C) of node i respectively, $Y_{\mathit{hij}}$ is the element in the i-th row and j-th column of the generalized node admittance matrix of the supply network $(m\cdot s)$, $p_{hj}$ is the pressure of node j in the supply network (MPa),$y_{bj}$ is the admittance $(m\cdot s)$ of branch j, $E_{bj}$ is the water pressure source parameter (Pa) of branch j, $A_{\mathit{hij}}$ is the element in the i-th row and j-th column of the node-branch incidence matrix of the supply network.

2.5 Coupling units

2.5.1 CHP unit model

Combined heat and power (CHP) units are electric-thermal coupling elements, which are special thermal power units that can generate electricity and provide heat to users (e.g., urban centralized heating) by utilizing waste steam (steam that has completed its work in turbines). During operation, the units exhibit the characteristic of “determining electricity by heat,” i.e., adjusting the electricity output according to changes in heat output. The mathematical model is as follows shown by Equation 14:

$P_{\mathit{CPH}}={\mathrm{\Phi }}_{\mathit{CPH}}/c_m(14)$

where $P_{\mathit{CPH}}$ and ${\mathrm{\Phi }}_{\mathit{CPH}}$ represent the electric and thermal output power of the CHP unit, respectively; $c_m$ is the unit heat-to-electricity ratio.

2.5.2 Gas turbine model

Gas turbines are electric-gas coupling elements, usually acting as loads in the natural gas network to convert and output electricity, serving as sources in the power grid. The gas turbine model is as follows shown by Equation 15 (Liu et al., 2022):

$V_f^{GA}=a_1{P}_{GA}^2+a_2{P}_{GA}+a_3(15)$

where $a_1$, $a_2$, and $a_3$ are fitting coefficients; $P_{GA}$ is the electric power of the gas turbine (MW); $V_f^{GA}$ is the gas consumption of the gas turbine (m³/h). In unified energy flow theory, gas flow is generally measured in kg/s. When the density of natural gas is 0.7174 kg/m³, the unit conversion for gas flow is as follows shown by Equation 16:

$1kg/s=3600/0.7174m^3/h(16)$

Gas turbines act as sources in the power grid and loads in the gas network, connecting the PV nodes of the power grid with the injection-defined nodes of the gas network. Assuming that the heat source S1 in the thermal network operates under constant outlet pressure and constant outlet temperature, and the balance node is numbered $n_1$, the thermal power of the connected CHP unit is determined by the state variables of the thermal network (Liu et al., 2022), as follows shown by Equation 17:

${\varphi }_i=c_p\left(T_{si}-T_{oi}\right)\left({\displaystyle \sum _{j=1}^{n_1}}{y}_{\mathit{hij}}{A}_{\mathit{hij}}{p}_{\mathit{hij}}-{\displaystyle \sum _{j=1}^{n_1}}{y}_{bj}{E}_{bj}\right)(17)$

where $T_{s{n}_1}$ represents the supply temperature of the balance node (°C); $n_1$ and $m_1$ are the number of nodes and branches in the thermal network, respectively. The electric power of the unit is determined by the characteristic of “determining electricity by heat” as shown by Equation 18:

${\varphi }_i=c_p\left(T_{s{n}_1}-T_{r{n}_1}\right)\left({\displaystyle \sum _{j=1}^{n_1}}{y}_{\mathit{hij}}{A}_{\mathit{hij}}{p}_{\mathit{hij}}-{\displaystyle \sum _{j=1}^{m_1}}{y}_{bj}{E}_{bj}\right)/cn(18)$

2.6 Power flow equations of ISES

Combining the mathematical models of the various subsystems in ISES and the energy coupling relationships of the coupling elements, the power flow equations of ISES can be represented by Equation 19.

$F=\left\{\begin{array}{c}{s}_i-{\displaystyle \sum _{i=1}^n}{\bar{Y}}_i{V}_i\bar{V}\hfill \\ s_i+c_p\left(T_{s{n}_1}-T_{r{m}_1}\right)\left({\displaystyle \sum _{j=1}^{n_1}}{y}_{\mathit{hjj}}{A}_{\mathit{hij}}{p}_{hi}-{\displaystyle \sum _{j=1}^{m_1}}{y}_{bj}{E}_{bj}\right)/c_m-{\displaystyle \sum _{i=1}^n}{\bar{Y}}_i{V}_i\bar{V}\hfill \\ s_i+P_{GA}-{\displaystyle \sum _{i=1}^n}\overline{Y_i}{V}_i{\bar{V}}_i\hfill \\ G_{gi}-{\displaystyle \sum _{j=1}^{n_2}}{Y}_{\mathit{gij}}{p}_{gj}\hfill \\ G_{gi}-V_f^{GA}-{\displaystyle \sum _{j=1}^{n_2}}{Y}_{\mathit{gij}}{p}_{gj}\hfill \\ {\varphi }_i-c_p\left(T_{si}-T_{oi}\right)\left({\displaystyle \sum _{j=1}^{n_1}}{y}_{\mathit{hij}}{A}_{\mathit{hij}}{p}_{hi}-{\displaystyle \sum _{j=1}^{m_1}}{y}_{bj}{E}_{bj}\right)\hfill \\ G_{hi}-\left({\displaystyle \sum _{j=1}^{n_1}}{y}_{hj}{A}_{\mathit{hij}}{p}_{\mathit{hij}}-{\displaystyle \sum _{j=1}^{m_1}}{y}_{bj}{E}_{bj}\right)\hfill \end{array}\right.(19)$

Then, using the coupling elements as bridges, the network equations of NGS and DHS can be combined with the network equations of EPS. Their power flow can be equivalently represented as active power P substituted into Equation 2 to obtain the equivalent power flow equations of ISES.

$S=P_m^{\mathit{CHP}}+P_k^{GA}+s={\displaystyle \sum _{i=1}^n}\overline{Y_i}{V}_i\overline{V_i},\forall i\in \left\{\begin{array}{c}\hfill \left(1,\dots ,n\right)\hfill \end{array}\right\}(20)$

where, $P_m^{\mathit{CHP}}$ represents the equivalence of the natural gas flow of NGS to the active power of EPS at the gas turbine.$P_k^{GA}$ represents the equivalence of the thermal energy of DHS to the active power of EPS at the combined heat and power unit. $S$ represents the equivalent load power in ISES. The natural gas transported by NGS, after being processed through the gas turbine model, is equivalent to a power source in EPS, with the coupling node in the power system set as m. Therefore, in ISES, there is electrical coupling at node m, converting the gas consumption in NGS into generation power at the coupling node. Similarly, the combined heat and power unit adjusts the generation power based on changes in heating power, which is equivalent to the load of EPS. Therefore, in ISES, there is electrothermal coupling at node k, converting thermal load into electrical load.

3 Sufficient analytical conditions for the existence of ISES equilibrium

In this section, we will apply Brouwer’s fixed-point theorem to the power flow equations of ISES to derive the static voltage stability criterion of ISES. First, taking the conjugate of both sides of Equation 3 gives Equation 21.

$\left[\left[\overline{V_L}\right]\right]YV=-\bar{S}(21)$

Equation 21 can be expanded based on Equation 1 and written in the form of Equation 22.

$-{\left[\left[\overline{V_L}\right]\right]}^{-1}\bar{S}=Y_{LL}{V}_L+Y_{LS}{V}_S(22)$

Multiplying both sides of the equation by $[[\overline{V_L}]]$ gives Equation 23.

$\left[\left[\overline{V_L}\right]\right]Y_{LL}{V}_L+\left[\left[\overline{V_L}\right]\right]Y_{LS}{V}_S+\bar{S}=0(23)$

Note that the admittance matrix $Y_{LL}$ in Equation 23 is not the complete matrix constructed for all buses of the network, but a submatrix obtained after removing the slack bus. This submatrix is non-singular and therefore invertible Equation 10. By left-multiplying both sides by $Y_{LL}^{-1}{[[\overline{V_L}]]}^{-1}$, Equation 23 can be rewritten in the form of Equation 24, where $V^{*}=-Y_{LL}^{-1}{Y}_{LS}{V}_S$.

$V_L=V^{*}-Y_{LL}^{-1}{\left[\left[\overline{V_L}\right]\right]}^{-1}\bar{S}(24)$

For an n-dimensional quadratic equation with n unknowns (MDQE), solving it is quite challenging. The key lies in transforming the power flow equations into a form that constructs a mapping, and then obtaining the solution of the equations by analyzing the fixed points of the mapping.

To construct the mapping, we multiply both sides of Equation 24 by ${[[V^{*}]]}^{-1}$ and define $x={[[V^{*}]]}^{-1}{V}_L$, thus obtaining Equation 25.

$x=1_n-\frac{{\left[\left[V^{*}\right]\right]}^{-1}{Y}_{LL}^{-1}{\left[\left[V^{*}\right]\right]}^{-1}\left[\left[\bar{S}\right]\right]}{\bar{x}}(25)$

Let $A={[[V^{*}]]}^{-1}{V}_{LL}^{-1}{[[V^{*}]]}^{-1}[[\bar{S}]]$, where $S=P_m^{\mathit{CHP}}+P_n+Q$. Thus, the mapping form of the power flow equations of ISES can be obtained as follows shown by Equation 26.

$x=\gamma \left(x\right)=1_n-A\frac{1}{\bar{x}}(26)$

Therefore, the solvability of Equation 23 is equivalent to the existence of a fixed point of the function $\gamma (x)$. For this purpose, an iterative algorithm for the multi-energy flow equations is constructed. According to Lemma 1, what we need to prove is that $\gamma (x)$ is a self-mapping on its domain of definition, where the domain of $\gamma (x)$ is two-dimensional, i.e., $S_r(x)≔\{x|{\Vert x-1_n\Vert }_{\infty }\le r\}$. Since $S_r(x)\in {\mathbb{C}}^n,r>0$, we can substitute $x=a+jb$ into $S_r(x)$, obtaining ${\Vert a-1+jb\Vert }_{\infty }\le r$. Next, based on the properties of the norm, we get $\{x|{(a-1)}^2+b^2\le r^2\}$, which represents a circle with center (1,0) and radius r, as shown in Figure 3. We obtain $min\left|\overline{x_i}\right|=1-r$. On the other hand, according to Lemma 2, we get Equation 27.

${\Vert \frac{A}{\bar{x}}\Vert }_{\infty }\le {\Vert A\Vert }_{\infty }{\Vert \frac{1}{\bar{x}}\Vert }_{\infty }={\Vert A\Vert }_{\infty }\frac{1}{min\left|\overline{x_i}\right|}(27)$

Figure 3

Figure 3. The domain of $\gamma (x)$.

Therefore, we can express the range as Equation 28.

$S_r\left(\gamma \left(x\right)\right)=\left\{\gamma \left(x\right)|\left(1_n-{\Vert A\Vert }_{\infty }\frac{1}{1-r}\right)\le \gamma \left(x\right)\le 1_n\right.(28)$

According to Brouwer’s fixed-point theorem, if there exists a non-empty compact set $S_r$ such that $\gamma (x)$ is a self-mapping, i.e., $S_r(\gamma (x))\subseteq S_r(x)$, then there exists a unique solution $\gamma (x^{*})=x^{*}$. Based on this idea, a sufficient solvability condition for Equation 23 is obtained. For any $x\in S_r$, when$S_r\left(\gamma (x)\subseteq S_r(x)\right.$ is satisfied, we obtain Equation 29.

$\gamma \left(x\right)\ge 1_n-{\Vert A\Vert }_{\infty }\frac{1}{1-r}\ge 1-r(29)$

After rearranging, we can obtain the static stability criterion for ISES, which is Equation 30.

$4{\Vert A\Vert }_{\infty }\le 1(30)$

In fact, the static voltage stability of both is fundamentally based on the existence of solutions to the power flow equations or multi-energy flow equations. Due to the coupling effects of multi-energy flows, the multi-energy flow equations are more complex than the power system flow equations. In this study, we modeled the Combined Heat and Power (CHP) components as pivotal connections among the Electric Power System (EPS), Natural Gas System (NGS), and District Heating System (DHS). Specifically, we consider the natural gas flow from the NGS as an input, while the outputs are the electrical energy from the EPS and thermal energy from the DHS. This approach allows us to integrate parameters such as gas flow and thermal energy into the power flow equations of the EPS (as illustrated in Equations 19, 20), enabling the derivation of the power flow equations for the ISES. Consequently, the analytical static voltage stability criterion for the ISES is proposed in Equation 30.

4 Cost optimization

Under the conditions of static voltage equilibrium, an optimization problem is set up with the objective function as shown by Equation 31:

$MinimizeC=C_{\mathit{gen}}+C_{\mathit{trans}}+C_{\mathit{cons}}(31)$

Where $P_{\mathit{gen}}=P_{\mathit{load}}$ ensures that the power generated by the system matches the load demand, $Q_{\mathit{gen}}=Q_{\mathit{load}}$ ensures that the thermal energy generated by the system is consistent with the demand, and $V_{\mathit{min}}\le V_i\le V_{\mathit{max}}$ ensures that the voltage remains within specified limits, respectively. Considering the output volatility of renewable energy sources, setting minimum and maximum output ranges can be established.

5 Simulation and experimental study

Based on the MATLAB/Simulink simulation platform, an IEEE 30-bus system as shown in Figure 4 is constructed. The blue nodes represent the slack buses, the red nodes represent the PV buses, and the purple nodes represent the PQ buses. The black line segments between the nodes represent the cables. The relevant parameters are shown in Table 3 and Table 4. Node 2 realizes electrical energy conversion through a gas turbine, nodes 5 and 8 achieve electrothermal energy conversion through CHP coupling elements, and nodes 11 and 13 are EPS power nodes.

Figure 4

Figure 4. IEEE-30 node standard test case diagram.

Table 3

Table 3. EPS node parameters.

Table 4

Table 4. EPS line parameters.

By calculating the various power flow data of ISES using the thermal-to-electrical conversion method, including the load power of EPS, four test cases are designed to verify the static stability of the system. Let $4{\Vert A\Vert }_{\infty }=\tau$. In this simulation, the electrothermal coupling modules, i.e., nodes 5 and 8 of the CHP units operating in thermal-to-electrical mode, consume 0.4316 pu of thermal energy, equivalent to $P_l=0.2878pu$ according to Table 5. The electrical coupling module, the gas turbine, consumes 8,898.3 m³/h of natural gas, equivalent to $P_g=0.3586pu$ according to Table 5.

Table 5

Table 5. Fitting coefficient and thermoelectric ratio.

Case 1. Nodes 2, 11, and 13 are source nodes, and all other loads are EPS loads with no load nodes converted from thermal loads.

Case 2. Nodes 2, 11, and 13 are source nodes, and nodes 5 and 8 are load nodes converted from thermal loads. All other loads are EPS loads.

Example 1. In the per-unit system, $V_s$ is the reference voltage. Therefore, $V_1=1.05\angle 0^{°}$, and the line cable impedance is shown in Table 5. For Case 1: when the load power is scaled up to 1.5 times that of Table 3, ${\tau }_1=0.9759<1$; for Case 2: when the load power is scaled up to 2.5 times that of Table 3, ${\tau }_1=1.5865$.

Example 2. In the per-unit system, $V_s$ is the reference voltage. Therefore, $V_1=1.05\angle 0^{°}$, and the line cable impedance is shown in Table 5. For Case 1: when the load power is scaled up to 1.52 times that of Table 3, ${\tau }_1=0.9812<1$; for Case 2: when the load power is scaled up to 2.3 times that of Table 3, ${\tau }_1=1.4694$.

Example 3. The reference voltage $V^{*}$ and the line cable impedance are the same as in Case 1; only the complex power of the load is changed. For Case 1: when the complex power is scaled up to 1.6 times that of Table 3, ${\tau }_2=0.9954<1$. For Case 2: when the load power is scaled up to 2 times that of Table 3, $\tau =1.3972>1$.

The iterative process of the algorithm (Example 1) is shown in Figure 5. Figure 5 illustrates that $x_{(n+1)}=\gamma (x_n)$ has a solution when ${\tau }_2<1$; otherwise, it may not have a solution.

Figure 5

Figure 5. The iteration processes of the proposed algorithm: $x_{(n+1)}=\gamma (x_n)$. (A) Example 1 (1). (B) Example 1 (2). (C) Example 2 (1). (D) Example 2 (2). (E) Example 3 (1). (F) Example 3 (2). (G) Example 4. (H) Example 5.

Example 4. A Simulink model is established. In this case, the actual power of the load is involved in the calculation, set p = 0.18 MW, Q = 0 Var, $V^{*}=380V$, and the line cable impedance is set to $R=0.2\mathrm{\Omega },X=2^{-3}H$, ${\tau }_3=0.9972<1$.

Example 5. Increase the DHS load power, set to p = 0.20 MW, Q = 0Var, $V^{*}=380V$, and the line cable impedance is the same as in Example 4. ${\tau }_4=1.1082>1$.

The results of Figures 5A, C, E indicate that if the static voltage stability criterion, i.e., Equation 30, is satisfied, then the proposed iterative algorithm converges to the solution; Figures 5B, D, F show that if the static voltage stability criterion is not satisfied, the proposed algorithm may diverge, and the power flow equations may have no solution. Furthermore, by comparing Case 1 and 2 in Figure 5, when DHS nodes are converted to electrical load nodes and added to the system, the system load increases, and the system’s stability margin decreases, meaning the system is more likely to lose its equilibrium point. Figures 5G, H show the results of the fixed-point algorithm’s iteration for Examples 4 and 5, respectively. In Figure 5G, the fixed-point algorithm converges, and the system has a solution; in Figure 5H, the fixed-point algorithm diverges, and the system may have no solution.

Simulation Examples 4, 5: In the case of Example 4, the three-phase voltage waveform in Figure 6A varies between $\pm$ 311 V. When the load is 0.18 MW, it satisfies the static voltage stability criterion. The power waveform in Figure 6B does not exhibit distortion, indicating that the system has an equilibrium point. This verifies the correctness of the proposed solvability condition. In the case of Example 5, the voltage waveform in Figure 6C oscillates at t = 0.92 s, and the power waveform in Figure 6D is distorted. At this time, the system is overloaded and does not satisfy the static voltage stability criterion, resulting in system instability.

In the coupled network, the NGS is coupled with the EPS through the gas turbine at node 2, making node 2 another power source node for the EPS. The DHS is coupled with the EPS through the combined heat and power units at nodes 5 and 8, making nodes 5 and 8 two load nodes for the EPS. NGS and DHS, as sources or loads for EPS through coupling equipment, will certainly affect the original power flow distribution of the EPS. Therefore, we will analyze the changes in the static voltage stability criterion and the equilibrium point of the system by separately increasing the natural gas consumption, electrical load, and thermal load in the IEEE 30-bus system. Through this data, we will infer the impact of natural gas consumption, electrical load, and thermal load on the stability of the IEEE 30-bus system.

Figure 7A shows that increasing natural gas consumption can increase the system’s stability margin to a certain extent. When the gas consumption increases to $20000m^3\cdot h$, the static voltage stability criterion decreases by $15.8\%$. As shown in Figure 7B, when the load power of EPS is increased to S = 0.04 pu, the static voltage stability criterion rises significantly. When S = 0.07 pu, the system stability margin is only 0.15, and this data is derived solely under the condition of electrical load. Next, we analyze the variation of the static voltage stability criterion with the change in thermal load when thermal load is added to the system. The CHP at the DHS node acts as a load at EPS node 5, meaning that changes in thermal load also reflect changes in the grid load. As shown in Figure 7C, when the thermal load increment reaches 70 MW, the static voltage stability criterion reaches 0.92, greatly reducing the system stability margin. When the DHS load power increment reaches 90 MW, the system stability margin approaches 0. Furthermore, comparing the load variations in Figures 7B, C, it can be seen that electrical load, as the main load of the IEEE 30-bus system, has a greater impact on the system’s stability.

As shown in Figure 8A, when the gas consumption increases, the equivalent generator at node 2 can supply more active power to the system, and the node voltage can be maintained within an appropriate range, thus preserving the equilibrium point. Conversely, as shown in Figure 9, the voltage magnitude of the IEEE 30-bus system decreases after the addition of DHS load. This is because the DHS nodes are connected to the grid through CHP, and the thermal load adds an extra load to the system, increasing the demand for active power, which leads to a voltage drop. When the thermal load increment reaches 100 MW, the node voltage magnitude drops to 1.001 pu, as shown in Figure 8B, which is $4.67\%$ less than the equilibrium node voltage magnitude.

Figure 6

Figure 6. Simulation results of Case 3, 4. (A) Three-phase load voltage, for Example 3. (B) Load power, for Example 3. (C) Three-phase load voltage, for Example 4. Load power, for Example 4.

Figure 7

Figure 7. Static voltage stability criteria var ISES with the increase in gas consumption and load. (A) Variation of the static voltage stability criterion with increasing natural gas consumption. (B) Variation of the static voltage stability criterion with increasing electrical load. (C) Variation of the static voltage stability criterion with increasing thermal load.

Figure 8

Figure 8. Equilibrium varies with the increase in gas consumption and thermal load. (A) Variation of the static voltage stability criterion with increasing natural gas consumption. (B) Variation of the equilibrium point with increasing thermal load.

Figure 9

Figure 9. The node voltage of 30-bus system.

6 Conclusion

Based on the characteristics of the power flow distribution in the EPS, NGS, and DHS models within ISES, this paper determines the power flow equations of ISES as shown in Figure 10. Using Brouwer’s fixed-point theorem, the power flow equations are algebraically transformed to obtain the corresponding algebraic equation mapping expressions. By constructing a self-mapping for the algebraic equation mapping expressions, the steady-state criterion for ISES is obtained. According to the relevant parameters and the steady-state criterion of ISES, it is determined whether ISES has a steady state. This paper successfully converts the solvability problem of ISES power flow equations into the existence problem of mapping fixed points. The derived static voltage stability criterion can quickly and accurately determine whether ISES has a steady state, providing support for subsequent state estimation, security analysis, and optimization control of ISES.

Figure 10

Figure 10. Flow diagram.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

JH: Conceptualization, Funding acquisition, Writing–original draft, Writing–review and editing. XD: Investigation, Project administration, Writing–original draft, Writing–review and editing. DJ: Resources, Software, Validation, Writing–review and editing, Writing–original draft. ZL: Methodology, Resources, Validation, Writing–original draft, Conceptualization, Writing–review and editing.

Funding

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

Conflict of interest

Authors JH, XD, DJ, and ZL were employed by NARI Technology Nanjing Control Systems Co., Ltd.

The reviewer LY declared a shared affiliation with the author DJ to the handling editor at the time of review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Dvijotham, K., Nguyen, H., and Turitsyn, K. (2018). Solvability regions of affinely parameterized quadratic equations. IEEE Control Syst. Lett. 2, 25–30. doi:10.1109/LCSYS.2017.2721380

CrossRef Full Text | Google Scholar

Huang, Y., Sun, Q., Li, Y., Gao, W., and Gao, D. W. (2022). A multi-rate dynamic energy flow analysis method for integrated electricity-gas-heat system with different time-scale. IEEE Trans. Power Deliv. 38, 231–243. doi:10.1109/tpwrd.2022.3186762

CrossRef Full Text | Google Scholar

Huang, Y., Sun, Q., Li, Y., Zhang, H., and Chen, Z. (2023). Adaptive-discretization based dynamic optimal energy flow for the heat-electricity integrated energy systems with hybrid ac/dc power sources. IEEE Trans. Automation Sci. Eng. 20, 1864–1875. doi:10.1109/TASE.2022.3188277

CrossRef Full Text | Google Scholar

Khalil, H. K. (2001). Nonlinear systems. London: Pearson.

Google Scholar

Li, C., Yang, H., Shahidehpour, M., Xu, Z., Zhou, B., Cao, Y., et al. (2020). Optimal planning of islanded integrated energy system with solar-biogas energy supply. IEEE Trans. Sustain. Energy 11, 2437–2448. doi:10.1109/TSTE.2019.2958562

CrossRef Full Text | Google Scholar

Liu, N., Tan, L., Sun, H., Zhou, Z., and Guo, B. (2022). Bilevel heat–electricity energy sharing for integrated energy systems with energy hubs and prosumers. IEEE Trans. Industrial Inf. 18, 3754–3765. doi:10.1109/TII.2021.3112095

CrossRef Full Text | Google Scholar

Liu, Z., Liu, R., Zhang, X., Su, M., Sun, Y., Han, H., et al. (2020). Feasible power-flow solution analysis of dc microgrids under droop control. IEEE Trans. Smart Grid 11, 2771–2781. doi:10.1109/tsg.2020.2967353

CrossRef Full Text | Google Scholar

Liu, Z., Su, M., Sun, Y., Yuan, W., Han, H., and Feng, J. (2018). Existence and stability of equilibrium of dc microgrid with constant power loads. IEEE Trans. Power Syst. 33, 6999–7010. doi:10.1109/tpwrs.2018.2849974

CrossRef Full Text | Google Scholar

Ma, Z., Zhou, Y., Zheng, Y., Yang, L., and Wei, Z. (2024). Distributed robust optimal dispatch of regional integrated energy systems based on admm algorithm with adaptive step size. J. Mod. Power Syst. Clean Energy 12, 852–862. doi:10.35833/MPCE.2023.000204

CrossRef Full Text | Google Scholar

Massrur, H. R., Niknam, T., Aghaei, J., Shafie-Khah, M., and Catalão, J. P. (2018). Fast decomposed energy flow in large-scale integrated electricity–gas–heat energy systems. IEEE Trans. Sustain. Energy 9, 1565–1577. doi:10.1109/tste.2018.2795755

CrossRef Full Text | Google Scholar

Molitierno, J. J. (2016). Applications of combinatorial matrix theory to laplacian matrices of graphs. London: Chapman & Hall.

Google Scholar

Nguyen, H. D., Dvijotham, K., Yu, S., and Turitsyn, K. (2018). A framework for robust long-term voltage stability of distribution systems. IEEE Trans. Smart Grid 10, 4827–4837. doi:10.1109/tsg.2018.2869032

CrossRef Full Text | Google Scholar

Peng, Q., Wang, X., Kuang, Y., Wang, Y., Zhao, H., Wang, Z., et al. (2021). Hybrid energy sharing mechanism for integrated energy systems based on the stackelberg game. CSEE J. Power Energy Syst. 7, 911–921. doi:10.17775/CSEEJPES.2020.06500

CrossRef Full Text | Google Scholar

Shabanpour-Haghighi, A., and Seifi, A. R. (2015). An integrated steady-state operation assessment of electrical, natural gas, and district heating networks. IEEE Trans. Power Syst. 31, 3636–3647. doi:10.1109/tpwrs.2015.2486819

CrossRef Full Text | Google Scholar

Sun, Q., Zhang, N., You, S., and Wang, J. (2019). The dual control with consideration of security operation and economic efficiency for energy hub. IEEE Trans. Smart Grid 10, 5930–5941. doi:10.1109/TSG.2019.2893285

CrossRef Full Text | Google Scholar

Wang, D., Huang, D., Hu, Q., Jia, H., Liu, B., and Lei, Y. (2024). Electricity-heat-based integrated demand response considering double auction energy market with multi-energy storage for interconnected areas. CSEE J. Power Energy Syst. 10, 1688–1700. doi:10.17775/CSEEJPES.2022.02140

CrossRef Full Text | Google Scholar

Wang, Z., Cui, B., and Wang, J. (2017). A necessary condition for power flow insolvability in power distribution systems with distributed generators. IEEE Trans. Power Syst. 32, 1440–1450. doi:10.1109/TPWRS.2016.2588341

CrossRef Full Text | Google Scholar

Wu, T., and Wang, J. (2022). Reliability evaluation for integrated electricity-gas systems considering hydrogen. IEEE Trans. Sustain. Energy 14, 920–934. doi:10.1109/tste.2022.3229896

CrossRef Full Text | Google Scholar

Yang, H., Li, M., Jiang, Z., and Zhang, P. (2020). Multi-time scale optimal scheduling of regional integrated energy systems considering integrated demand response. IEEE Access 8, 5080–5090. doi:10.1109/ACCESS.2019.2963463

CrossRef Full Text | Google Scholar

Zhang, S., Wang, S., Zhang, Z., Lyu, J., Cheng, H., Huang, M., et al. (2021). Probabilistic multi-energy flow calculation of electricity–gas integrated energy systems with hydrogen injection. IEEE Trans. Industry Appl. 58, 2740–2750. doi:10.1109/tia.2021.3094487

CrossRef Full Text | Google Scholar