Figure 1 illustrates the basic framework based on security regions, considering the collaborative support of multiple city systems. It is specifically divided into two parts: characterizing the low-carbon feasible space of each city system based on security regions, and the IEHGS scheduling model that accounts for the collaborative support of multiple city systems. Details are provided as follows:

To fully utilize energy, city systems are coupled with distributed renewable energy sources, combined heat and power units and micro gas turbines fueled by Hydrogen Enriched Compressed Natural Gas (HCNG), grid-connected electric vehicles, and various types of energy storage devices. Considering the low-carbon operation constraints and renewable energy uncertainties of city systems, a security region model is constructed to effectively quantify the low-carbon feasible space of city systems.

The security region model constructed above contains random variables, making it an uncertainty model that is difficult to solve directly. This section uses the distribution-ally robust chance-constrained method to handle the model, transforming the original uncertain security region model into a deterministic one. On this basis, based on the boundary points of the security region, the low-carbon feasible space of the city system is analytically represented using a multi-dimensional piecewise linear approximation method.

The model aims to minimize the total operating cost C of the provincial power grid IEHGS. The total cost C includes: the handling cost of renewable energy uncertainty C Uncer deal , electricity purchase cost C Power buy , compensation cost for city system support C LCFS supp , and carbon emission cost C CO 2 emi . The details are as follows: min   C = C Uncer deal + C Power buy + C LCFS supp + C CO 2 emi (21) Where each cost can be calculated using Equations 21–25. C Uncer deal = ∑ t ∈ Ω T ∑ m ∈ Ω GS a GS U F m . t GS . U + a GS D F m . t GS . D + ∑ t ∈ Ω T ∑ w ∈ Ω GFU a GFU U P w . t GFU . U + a GFU D P w . t GFU . D (22) Where Ω T represents the set of scheduling periods, Ω G S and Ω G F U refer to the sets of gas sources and gas turbines, respectively. m denotes the gas node, F m . t GS . U and F m . t GS . D represent the upper and lower control quantities of the gas source, respectively. a GS U and a GS D are the unit control costs, while www is the equipment index. P w . t GFU . U and P w . t GFU . D represent the upper and lower control quantities of the gas turbine, with a GFU U and a GFU D being the control costs. C Power buy = ∑ w ∈ Ω Rene ∑ t ∈ Ω T a Rene buy P w . t Re + ξ w . t Re + ∑ i ∈ Ω slack ∑ t ∈ Ω T a UPG buy P i . t UPG + − a UPG sell P i . t UPG − (23) Where Ω Rene and Ω s l a c k represent the sets of renewable energy stations and the slack node, respectively. P w . t Re and ξ w . t Re represent the predicted output and prediction error of the renewable energy stations, respectively. P i . t UPG + and P i . t UPG − represent the power purchased from and sold to the distribution network root node i and the upper-level grid, respectively. a Rene buy , a UPG buy , and a UPG sell are the feed-in tariffs for renewable energy stations, and the purchase and sale of electricity prices for the distribution network and upper-level grid, respectively. C LCFS supp = ∑ k ∈ Ω Park a Ele sup   E Ele . k sup + a H 2 sup   E H 2 . k sup + a NG sup   E NG . k sup (24) Where Ω P a r k represents the set of urban systems. a Ele sup   , a H 2 sup   , and a NG sup   are the compensation prices for the energy support provided by the urban systems in terms of electricity, hydrogen, and natural gas, respectively. E Ele . k sup , E H 2 . k sup , and E NG . k sup represent the electricity, hydrogen, and natural gas support amounts provided by the urban systems, respectively. C CO 2 emi = a CO 2 ce E car total − E car base (25) Where a CO 2 ce represents the carbon trading price, E car total and E car base are the total carbon emissions and the free carbon emission allowance for the regional IEHGS, respectively. E car base can be further expressed as: E car total = ∑ t ∈ Ω T ∑ m ∈ Ω GS e car Gas F m . t GS . E + η m − i GFU ξ i . t Re + ∑ i ∈ Ω slack e car UPG P i . t UPG + − ∑ i ∈ Ω slack e car UPG P i . t UPG − (26) Where η m − i GFU ξ i . t Re represents the additional output from the gas source to cope with the uncertainty of renewable energy, while e car Gas and e car UPG are the carbon emission intensities for the gas source and the upper-level grid, respectively.

Based on the analytical expressions of the low-carbon feasible space of urban systems in Section 2.3, further consideration is given to the low-carbon safe operation constraints of the provincial-level IEHGS, as detailed below.

Gas turbines and gas sources can serve as flexible resources to address the uncertainty of renewable energy. Considering the impact of renewable energy uncertainty, the distributional robust chance constraints for the gas turbine power output and reserve capacity are constructed as Formulas 27, 28; similarly, the distributional robust chance constraints for the gas source’s gas output and reserve capacity are constructed as Formulas 29, 30. Pr ∼ P Re P i . lo GFU ≤ P i . t GFU . E − ξ i . t Re ≤ P i . up GFU ≥ 1 − β Re , ∀ P Re ∈ R Re (27) Pr ∼ P Re − P i . t GFU . D ≤ − ξ i . t Re ≤ P i . t GFU . U ≥ 1 − β Re , ∀ P Re ∈ R Re (28) Pr ∼ P Re F m . lo GS ≤ F m . t GS . E + η m − i GFU ξ i . t Re ≤ F m . up GS ≥ 1 − β Re , ∀ P Re ∈ R Re (29) Pr ∼ P Re − F m . t GS . D ≤ η m − i GFU ξ i . t Re ≤ F m . t GS . U ≥ 1 − β Re , ∀ P Re ∈ R Re (30) Where P Re is the probability distribution function of the renewable energy uncertainty random variable ξ i . t Re , R Re is the fuzzy set of P Re , 1 − β Re is the confidence level, P i . t GFU . E and η m − i GFU are the planned active power and energy conversion coefficient of the gas turbine, respectively. P i . t GFU . E − ξ i . t Re and F m . t GS . E + η m − i GFU ξ i . t Re represent the actual outputs of the gas turbine and gas source point under the impact of uncertainty. P i . up GFU and P i . lo GFU are the power limits of the gas turbine, and F m . up GS and F m . lo GS are the gas output limits of the gas source point.

The node voltage constraint is shown in Equation 31, the upper grid interaction power constraint is shown in Equation 32, and the power flow and branch capacity constraints are shown in Equations 33–37. Considering the mitigation effect of gas turbines and gas sources on renewable energy uncertainty, the node power balance constraints are shown in Equations 38, 39. U i . ⁡ min ≤ U i . t ≤ U i . ⁡ max (31) − P i . up grid ⩽ P i . t grid ⩽ P i . up grid − Q i . up grid ⩽ Q i . t grid ⩽ Q i . up grid (32) P i . t Grid + P i . t GFG + P i . t DG − P i . t Load − P i . t P 2 G = ∑ j ∈ i P i j . t (33) Q i . t Grid + Q i . t GFG + Q i . t DG − Q i . t Load = ∑ j ∈ i Q i j . t (34) P i j . t = U i . t − U j . t r i j / r i j 2 + x i j 2 + θ i . t − θ j . t x i j / r i j 2 + x i j 2 (35) Q i j . t = U i . t − U j . t x i j / r i j 2 + x i j 2 − θ i . t − θ j . t r i j / r i j 2 + x i j 2 (36) P i j . t 2 + Q i j . t 2 ≤ S i j cap 2 (37) P i . t UPG + P i . t Re + P i . t GFU . E + ∑ k ∈ Ω i . Park P ELE . t . k LINE − P i . t P 2 G − P i . t Load = ∑ j ∈ Ω bus . i P i j . t (38) Q i . t UPG + Q i . t Re + Q i . t GFU − Q i . t Load = ∑ j ∈ Ω bus . i Q i j . t (39) Where Ω b u s . i and Ω i . P a r k represent the set of distribution network nodes connected to node i and the set of city systems, respectively. P i j . t and Q i j . t are the active and reactive power transmitted through branch ij, respectively. P i . t Load and Q i . t Load are the active and reactive loads at the node, respectively. P i . t UPG and Q i . t UPG are the active and reactive power exchanged between the upper grid and node i, respectively. P i . t P 2 G represents the electricity consumption of P2G. Q i . t Re and Q i . t GFU are the reactive power of renewable energy and gas turbines, respectively.

The IEHGS scheduling model established above can be written in the following form: min x ∈ X F x , ξ s . t .   Pr ∼ p Re f 1 x , ξ ≤ 0 ≥ 1 − β Re , ∀ P Re ∈ R Re ⏟ D i s t r i b u t i o n a l l y   r o b u s t   c h a n c e   c o n s t r a i n t s f 2 x ≤ 0 ⏟ A n a l y t i c a l   e x p r e s s i o n   o f   t h e   l o w − c a r b o n   f e a s i b l e   s p a c e   o f   u r b a n   s y s t e m s , f 3 x ≤ 0 ⏟ D e t e r m i n i s t i c   c o n s t r a i n t s (40) Where x represents all decision variables, and X is the feasible region of x . ξ is a random variable. F x , ξ is the unified expression of the objective function Formulas 21–26; the constraints involving random variables in the scheduling model (i.e., Formulas 27–30) are uniformly represented as distributionally robust chance constraints; the analytical expression of the low-carbon feasible space of the urban system in Section 2.3 is uniformly represented as f 2 x ≤ 0 ; other remaining constraints (such as Formulas 38, 39) are uniformly represented as f 3 x ≤ 0 . F x , ξ , f 1 x , ξ , and f 3 x ≤ 0 contain random variables or nonlinear constraints, which increases the difficulty of solving the model. To address this, the following transformation method is proposed.

1) Deterministic transformation of the objective function min x ∈ X F x , ξ

Based on distributionally robust chance constraints, the objective function in Equation 40 can be equivalently transformed into: min x ∈ X F x , ξ = min x ∈ X max P Re ∈ R Re E P F x , ξ (41)

Based on Theorem 4 from reference (Hu and Hong, 2012), Equation 41 can be equivalently transformed into as shown in Equation 42: min x ∈ X , ν Re ≥ 0 ν Re InE P Re . 0 e F x , ξ / ν Re + ν Re κ Re (42) Where P Re . 0 is the reference distribution of the random variable, κ Re is the KL divergence tolerance, and ν Re is the auxiliary variable.

Further, F x , ξ can be expressed as F x , ξ = d T x + h T x + g , where d and h are the coefficient vectors of x and ξ in Equations 21–26, and g is the constant term in Equations 21–26. Therefore, the following expression holds: min x ∈ X , ν Re ≥ 0 d T x + g + ∑ i ∈ Ω h ∑ t ∈ Ω T h i 2 σ i . t Re 2 2 ν Re + h i μ i . t Re + ν Re κ Re (43)

For the variable ν Re in Equation 43, when ν Re = ∑ i ∈ Ω h ∑ t ∈ Ω T h i 2 σ i . t Re 2 2 κ Re (44)

Then, the terms containing ν Re in Equation 43 can attain their minimum value by Equation 44, which can be expressed as shown in Equation 45: min ν Re ≥ 0 ∑ i ∈ Ω h ∑ t ∈ Ω T h i 2 σ i . t Re 2 2 ν Re + ν Re κ Re = 2 κ Re ∑ i ∈ Ω h ∑ t ∈ Ω T h i 2 σ i . t Re 2 (45)

Based on this, Equation 43 can be further transformed into Equation 46: min x ∈ X d T x + g + ∑ i ∈ Ω h ∑ t ∈ Ω T h i μ i . t Re + 2 κ Re ∑ i ∈ Ω h ∑ t ∈ Ω T h i 2 σ i . t Re 2 (46)

2) Deterministic transformation of the distributionally robust chance constraint Pr ∼ P Re f 1 x , ξ ≤ 0 ≥ 1 − β Re , ∀ P Re ∈ R Re

The constraints represented by the distributionally robust chance constraint in Equation 40 can be further expressed as: Pr ∼ P Re a x Re + b ξ i . t Re + c ≤ 0 ≥ 1 − β Re , ∀ P Re ∈ R Re (47) Where x Re a, b, and c represent the decision variables, coefficients of the decision variables, coefficients of the random variables, and constant terms in Equations 27–30, respectively.

Based on the equivalent transformation method for fuzzy chance-constrained problems, adjusting the risk parameter β Re in Equation 47 allows the fuzzy chance constraint to be converted into the following chance constraint: Pr ∼ P Re . 0 a x Re + b ξ i . t Re + c ≤ 0 ≥ 1 − β ¯ Re (48) Where β ¯ Re is the adjusted risk parameter, which is associated with β Re and κ Re and can be obtained using a binary search algorithm.

Based on the Bernstein approximation method, Equation 48 can be further converted to Equation 49: inf α Re > 0 α Re In E P Re . 0 e a x Re + b ξ i . t Re + c / α Re − α Re In β ¯ Re ≤ 0 (49) Where α Re is an auxiliary variable.

3) Linearization of nonlinear constraint f 3 x ≤ 0

The nonlinearity in f 3 x ≤ 0 arises from branch capacity constraints, which can be handled using a piecewise linearization method to reduce the computational burden.

At this point, the original non-convex, nonlinear scheduling model with uncertainty has been transformed into a mixed-integer linear optimization model, which can be solved directly using commercial solvers.

To investigate the impact of the collaborative support function of multi-city systems on the low-carbon operation of the IEHGS, the following two comparative cases are set up to analyze the collaborative support function of multiple city systems in IEHGS_1.

Case 1:

The IEHGS scheduling model proposed in this paper, formulates a low-carbon operation strategy for IEHGS_1 based on security regions and considers the collaborative support function of multi-city interconnected power systems.

Case 2:

The low-carbon operation strategy for IEHGS_1 without considering the support function of city systems (Hu et al., 2021).

The operating costs of IEHGS_1 in Cases 1, 2 are shown in Table 2. Compared to Case 2, the total operating cost of IEHGS_1 in Case 1 decreases by 7.21%. Specifically, the electricity purchase cost of IEHGS_1 is reduced by 60.78%, and no carbon emission costs are incurred, allowing the sale of surplus carbon allowances for profit. This result demonstrates that the proposed IEHGS low-carbon operation strategy (Case 1), which considers the collaborative support of multiple city systems, effectively reduces electricity purchase costs and carbon emission expenses, thereby lowering the total operating cost of IEHGS. This outcome validates the effectiveness of incorporating city system collaboration into IEHGS dispatching as proposed in this paper.

The optimal support points provided by each city system in Case 1 are located within their own low-carbon feasible space, as shown in Figure 4. This result indicates that each city can collaboratively utilize different low-carbon resources and system flexibility resources, effectively assisting in the heterogeneous energy conversion of IEHGS_1. This provides heterogeneous energy support for the low-carbon operation of IEHGS_1, thereby reducing its energy procurement costs.

In addition, the provision of clean, heterogeneous energy from multiple cities to IEHGS_1 reduces its consumption of fossil fuel (natural gas), thereby contributing to lower carbon emissions and reducing the need for carbon allowance trading. As shown in Table 3, the carbon emissions of IEHGS_1 in Case 1 are reduced by 24.70%. Since there is no requirement for additional carbon allowances, the surplus allowances can be sold to make a profit. Therefore, considering the collaborative support of multiple city systems also effectively reduces the carbon emission costs of IEHGS_1.

Based on this, the various city systems of IEHGS_1 can be characterized by the low-carbon feasible space defined by security regions. While exploring the low-carbon operating potential of the city systems, they provide collaborative support for the low-carbon operation of IEHGS_1, enhancing its low-carbon performance and improving the supply guarantee capability of the multi-city interconnected power system.

To verify the effectiveness of the security region in characterizing the low-carbon feasible space of city systems and providing support for the low-carbon operation of the provincial grid IEHGS, the following two case studies are set up for comparative analysis.

Case 3:

The low-carbon operation strategy of IEHGS_1 is developed based on the collaborative support role of the multi-city interconnected power system considered by the security region, as proposed in this paper.

Case 4:

The low-carbon operation strategy of IEHGS_1 is developed based on the collaborative support role of the multi-city system, considering the interconnection line energy exchange capacity, as proposed in reference (Yang and Liu, 2023).

The optimal support points provided by each city system in Case 4 are shown in Figure 5. It can be observed that the optimal support points in Case 4 are located outside the low-carbon feasible space, meaning these support points are infeasible. This is because Case 4 only considers the transmission capacity constraints of the interconnection lines and does not account for the internal security operation constraints of the city systems, resulting in the optimal support points being unsuitable for the city systems. However, the support points of each park system in Case 3 are all located within the low-carbon feasible space (as shown in Figure 4), which verifies the effectiveness of considering the collaborative support role of multi-city systems based on the security region.