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.

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. V L ̄ Y V = − S ̄ (21)

Equation 21 can be expanded based on Equation 1 and written in the form of Equation 22. − V L ̄ − 1 S ̄ = Y L L V L + Y L S V S (22)

Multiplying both sides of the equation by [ [ V L ̄ ] ] gives Equation 23. V L ̄ Y L L V L + V L ̄ Y L S V S + S ̄ = 0 (23)

Note that the admittance matrix Y L L 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 L L − 1 [ [ V L ̄ ] ] − 1 , Equation 23 can be rewritten in the form of Equation 24, where V * = − Y L L − 1 Y L S V S . V L = V * − Y L L − 1 V L ̄ − 1 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 − V * − 1 Y L L − 1 V * − 1 S ̄ x ̄ (25)

Let A = [ [ V * ] ] − 1 V L L − 1 [ [ V * ] ] − 1 [ [ S ̄ ] ] , where S = P m 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 = γ x = 1 n − A 1 x ̄ (26)

Therefore, the solvability of Equation 23 is equivalent to the existence of a fixed point of the function γ ( 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 γ ( x ) is a self-mapping on its domain of definition, where the domain of γ ( x ) is two-dimensional, i.e., S r ( x ) ≔ { x | ‖ x − 1 n ‖ ∞ ≤ r } . Since S r ( x ) ∈ C n , r > 0 , we can substitute x = a + j b into S r ( x ) , obtaining ‖ a − 1 + j b ‖ ∞ ≤ r . Next, based on the properties of the norm, we get { x | ( a − 1 ) 2 + b 2 ≤ r 2 } , which represents a circle with center (1,0) and radius r, as shown in Figure 3. We obtain m i n x i ̄ = 1 − r . On the other hand, according to Lemma 2, we get Equation 27. ‖ A x ̄ ‖ ∞ ≤ ‖ A ‖ ∞ ‖ 1 x ̄ ‖ ∞ = ‖ A ‖ ∞ 1 m i n x i ̄ (27)

Therefore, we can express the range as Equation 28. S r γ x = γ x | 1 n − ‖ A ‖ ∞ 1 1 − r ≤ γ x ≤ 1 n (28)

According to Brouwer’s fixed-point theorem, if there exists a non-empty compact set S r such that γ ( x ) is a self-mapping, i.e., S r ( γ ( x ) ) ⊆ S r ( x ) , then there exists a unique solution γ ( x * ) = x * . Based on this idea, a sufficient solvability condition for Equation 23 is obtained. For any x ∈ S r , when S r γ ( x ) ⊆ S r ( x ) is satisfied, we obtain Equation 29. γ x ≥ 1 n − ‖ A ‖ ∞ 1 1 − r ≥ 1 − r (29)

After rearranging, we can obtain the static stability criterion for ISES, which is Equation 30. 4 ‖ A ‖ ∞ ≤ 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.

Under the conditions of static voltage equilibrium, an optimization problem is set up with the objective function as shown by Equation 31: M i n i m i z e C = C gen + C trans + C cons (31) Where P gen = P load ensures that the power generated by the system matches the load demand, Q gen = Q load ensures that the thermal energy generated by the system is consistent with the demand, and V min ≤ V i ≤ V 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.

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.

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 ‖ A ‖ ∞ = τ . 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.2878 p u 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.3586 p u according to Table 5.

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 ∠ 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, τ 1 = 0.9759 < 1 ; for Case 2: when the load power is scaled up to 2.5 times that of Table 3, τ 1 = 1.5865 .

Example 2

In the per-unit system, V s is the reference voltage. Therefore, V 1 = 1.05 ∠ 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, τ 1 = 0.9812 < 1 ; for Case 2: when the load power is scaled up to 2.3 times that of Table 3, τ 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, τ 2 = 0.9954 < 1 . For Case 2: when the load power is scaled up to 2 times that of Table 3, τ = 1.3972 > 1 .

The iterative process of the algorithm (Example 1) is shown in Figure 5. Figure 5 illustrates that x ( n + 1 ) = γ ( x n ) has a solution when τ 2 < 1 ; otherwise, it may not have a solution.

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 * = 380 V , and the line cable impedance is set to R = 0.2 Ω , X = 2 − 3 H , τ 3 = 0.9972 < 1 .

Example 5

Increase the DHS load power, set to p = 0.20 MW, Q = 0Var, V * = 380 V , and the line cable impedance is the same as in Example 4. τ 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 ± 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 20000 m 3 ⋅ 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.

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.