Front. Energy Res. Frontiers in Energy Research Front. Energy Res. 2296-598X Frontiers Media S.A. 1363822 10.3389/fenrg.2024.1363822 Energy Research Original Research Station-network cooperative planning method of urban integrated energy system based on energy flow model Fu et al. 10.3389/fenrg.2024.1363822 Fu Xin 1 Yu Shunjiang 2 * He Qibo 1 Wang Long 1 Niu Chengeng 2 Yang Li 2 1 State Grid Wuxi Power Supply Company of Jiangsu Electric Power Co., Ltd., Wuxi, China 2 College of Electrical Engineering, Zhejiang University, Hangzhou, China

Edited by: Shengyuan Liu, State Grid Zhejiang Electric Power Co., Ltd., China

Reviewed by: Zhaobin Du, South China University of Technology, China

Xun Dou, Nanjing Tech University, China

 *Correspondence: Shunjiang Yu, yushunjiang@zju.edu.cn 22 05 2024 2024 12 1363822 31 12 2023 25 04 2024 Copyright © 2024 Fu, Yu, He, Wang, Niu and Yang. 2024 Fu, Yu, He, Wang, Niu and Yang

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Coordinated siting and sizing for energy stations and supply networks in urban integrated energy system (UIES) is significant for economic improvement and carbon emissions reduction. A station-network cooperative planning method of UIES based on energy flow model is proposed. First, an operation model of heat network based on energy flow theory is proposed, which solves the problem that the temperature mixing equation in the traditional operation model cannot be applied to heat network planning. On this basis, a bi-level model for station-network cooperative planning of UIES is constructed, in which the upper level optimizes the siting and sizing of the energy station and the topology of the supply network, and the lower level optimizes the operation of the UIES and feeds back the operation cost of the UIES to the upper level. Finally, a solution method of cooperative planning based on the Karush-Kuhn-Tucher condition is proposed, to transform the bi-level nonlinear optimization model into a single-level linear optimization model for efficient solution. Case studies on the 55-node and 77-road urban topology show that the proposed method can perform an effective planning on energy supply network topology and rationally configure the capacity of various devices in the energy station.

 energy flow model urban integrated energy system bi-level planning model station-network cooperative optimization heat network section-at-acceptance Smart Grids 1 Introduction

With the continuous depletion of traditional fossil fuels and the increasingly serious problem of environmental pollution, there is an urgent need to form a new highly efficient, clean and sustainable way of utilizing multiple energy sources in a complementary manner (Strezoski et al., 2022). As an advanced energy utilization concept emerging in recent years, urban integrated energy system (UIES) can realize the coordinated complementary and efficient utilization of various types of energy sources, thereby increasing the rate of renewable energy consumption and reducing carbon emissions (Li et al., 2021). Therefore, it is of great significance to study the optimal configuration and operation optimization strategy of UIES under the background of dual carbon target (Chen et al., 2022).

Currently, global scholars have conducted a large number of studies on the planning of UIES. In Farrokhifar et al. (2020), a comprehensive summary of the planning research of UIES is provided from the perspectives of mathematical modeling, operation constraints, optimization objectives and methods of handling uncertainty. In Mu et al. (2020), a planning model for energy conversion and storage devices in UIES considering dynamic energy con-version efficiency coefficients is developed. In order to achieve carbon emission reduction, a device capacity planning strategy in UIES based on various device investment constraints is proposed in Wang et al. (2019). In Dong et al. (2023), an optimal expansion planning model for an integrated energy system consisting of a power grid, a gas network and multiple energy hubs is proposed. The references mentioned above focus on the planning problem of the capacity of the device in energy stations when the siting of the station and energy supply network are predetermined. However, the siting and sizing of the energy station and the energy supply network are interacting with each other, and the global optimal solution of UIES planning cannot be obtained when the siting of the station and the supply network are predetermined (Xiao et al., 2018). Therefore, further research is needed to investigate the cooperative planning of the siting and sizing of energy station with energy supply network.

In this context, a coordinated siting and sizing method for PIES that considers load complementary characteristics is proposed in Liu. (2020a). In Liu (2020b), an alternating optimization method based on cellular network theory is proposed for UIES’s energy station site selection, energy supply area division, and pipe network planning, and the Kruskal algorithm is used to determine the pipe network topology planning strategy. In Zhang et al. (2015), a long-term expansion planning methodology for energy hubs containing electric, gas and heat to optimally determine the least-cost planning schedules for generating units, transmission lines, gas boiler (GB) and combined heat and power (CHP) generation is proposed. However, these references do not use an exact model that includes variables such as thermal mass flow and temperature when modeling heat networks, but simply de-scribes them through thermal power balance constraints only, which will lead to inability to obtain accurate planning results.

To address this problem, a planning method of UIES consisting of a simulation model and an optimization model is proposed in Hong et al. (2018), where the simulation model adopts a two-layer bus structure, and considers both the energy flow in the external bus and the detailed network in the internal bus to achieve an accurate modeling of the thermal network in mass regulation mode. In Chen et al. (2024), a two-stage planning model of UIES for coordinating the location and size of the energy station as well as the supply net-work is proposed, where the optimal location of the energy station and the topology of the supply network are optimized to obtain the optimal location of the energy station in the first stage, and the size of the energy station and the pipeline network selection are optimized to obtain the optimal location of the energy station and the supply net-work topology in the second stage. However, all of the above references are planning with the topology of the heat network determined and does not consider planning for the heat supply pipeline. This is because the temperature mixing equation is only valid if the inlet and outlet of the heat pipeline are correctly determined, whereas the ther-mal mass flow direction cannot be predetermined when considering the planning of the heat pipeline, and therefore the model needs to be improved.

In this regard, by introducing an auxiliary heat variable and approximating the calculation equation of heat loss, the heat network model is transformed into a linear energy flow model in Xue et al. (2021), which can be applied to the planning problem of heat net-work to significantly improve the computational efficiency. Based on this model, for a mesh multi-energy distribution network (MMDN) containing electric, gas and heat, a resilience-oriented extended planning and reinforcement model is proposed in Li et al. (2023), to enhance the resilience of the MMDN through multi-energy coupling support, reinforcement, and reconfiguration strategies to resist low-probability high-damage un-foreseen events. In order to increase the penetration rate of wind power, a coordinated planning model for electric and heat based on the energy flow model and considering seasonal heat network reconfiguration is developed in Du et al. (2024), which is formulated as a mixed integer linear programming model. However, none of the energy flow-based heat network operation models in the above reference take into account the transmission delay characteristic of thermal mass, which cannot fully utilize the heat storage capacity of the heat network.

In summary, it can be seen that how to apply the energy flow model to the station-network cooperative planning of UIES and consider the transmission delay of thermal mass still needs to be investigated. Aiming at this problem, a station-network cooperative planning method of UIES based on energy flow model is proposed, which realizes the cooperative optimization of the siting and sizing of energy station and the topology of energy supply network. The main contributions are as follows:

(1) By introducing energy flow variables to replace the product of flow and temperature, and establishing the transmission delay characteristic equation, an energy flow-based heat network operation model is proposed, which solve the problem that the existing models do not take into account the transmission delay characteristic of thermal mass.

(2) A bi-level model for station-network cooperative planning of UIES is constructed, in which the upper level optimizes the siting and sizing of energy station and the topology of energy supply network with the objective of minimizing the total cost and passes it to the lower level, and the lower level optimizes the operation of the UIES according to the planning scheme and feeds back the total operation cost of the UIES to the upper level.

(3) A solution method of cooperative planning based on the Karush-Kuhn-Tucher (KKT) condition is proposed. Firstly, the nonlinear bi-level cooperative planning model is transformed into a bi-level linear model based on mathematical linearization theory to improve the solution efficiency, and then the bi-level linear model is transformed into a mixed-integer linear planning model which is easy to be solved based on the KKT condition.

The remainder is organized as follows. In Section 2, the operation model of heat network based on energy flow theory is proposed. In Section 3, the bi-level model for station-network cooperative planning of UIES is constructed. In Section 4, the solution method of cooperative planning based on the KKT condition is proposed. In Section 5, case studies are conducted to demonstrate the Effectiveness and efficiency of the proposed method. Conclusion is finally given in Section 6.

 2 Operation model of heat network based on energy flow theory 2.1 Exact operation model of heat network

In the heat network, the return water is heated by the heat source, and enters the supply pipeline after the temperature rises, and is transmitted to each secondary heat exchanger station through the supply pipeline, and flows into the return pipeline again after heat exchange. The exact operation model of heat network is as follows (Wu et al., 2018; Ha et al., 2022): H i , t , s N = c w m i T i , t , s S − T i , t , s R ∑ k i m k i S = ∑ i j m i j S ∑ k i m k i R = ∑ i j m i j R T i , t , s S ∑ k i m k i S = ∑ k i m k i S T k i , t S , out T i , t , s R ∑ k i m k i R = ∑ k i m k i R T k i , t , s R , out T i , t , s S = T i j , t , s S , in T i , t , s R = T i j , t , s R , out T i j , t , s S , out = T t , s am + T i j , t − τ i j , s S , in − T t , s am e − λ i j L i j m i j S A i j ρ w c w T i j , t , s R , out = T t , s am + T i j , t − τ i j , s R , in − T t , s am e − λ i j L i j m i j R A i j ρ w c w T min ⁡ , i S ≤ T i , t , s S ≤ T max ⁡ , i S T min ⁡ , i R ≤ T i , t , s R ≤ T max ⁡ , i R m min ⁡ , i j S ≤ m i j S ≤ m max ⁡ , i j S m min ⁡ , i j R ≤ m i j R ≤ m max ⁡ , i j R where H N i,t,s is the heat power of heat network node i at period t under scenario s; c w is the specific heat capacity of water; m i is the mass flow rate of heat exchanger at node i; T S i,t,s and T R i,t,s are the temperature of water flow at node i in supply/return network at period t under scenario s, respectively; m S ij and m R ij is the mass flow rate in pipe ij of supply and return networks, respectively; λ ij is the thermal conductivity of pipe ij; L ij is the length of the road between nodes i and j; A ij is the cross-sectional area of pipe ij; ρ w is the density of water; T S, in ij,t,s and T ij , t , s R,in are the inlet temperatures of water flow in pipe ij of supply and return network at period t under scenario s, respectively; T ij , t , s S,out and T ij , t , s R,out are the outlet temperatures of water flow in pipeline ij of supply and return network at period t under scenario s, respectively; Tam t,s is the ambient temperature at period t under scenario s; τ ij is the transmission delay of pipeline ij; T S max,i and T S min,i are the maximum and minimum temperatures of water in supply network for node i, respectively; T max , i R and T min , i R are the maximum and minimum temperatures of water in return network for node i, respectively; m S max,ij and m min , ij S are the maximum and minimum flow rates of the water in pipe ij of supply network, respectively; m R max,ij and m R min,ij are the maximum and minimum flow rates of the water in pipe ij of return network, respectively.

The above model cannot be directly applied to the station-network cooperative planning of UIES mainly because of the following two problems:

(1) In the exact operation model, the temperature mixing equations, i.e., Eqs. 4, 5, are only valid if the inlet and outlet of the heat pipe are correctly determined, whereas the direction of mass flow cannot be predetermined when considering the planning for heat network, and therefore the model cannot be used.

(2) The product term of mass flow rate and temperature exists in the exact model, which, if applied to the station-network cooperative planning of UIES, will make the planning model become a nonlinear optimization model, which not only has a poor solution efficiency, but may not be able to achieve the solution when the size of network is large.

Therefore, a heat network operation model based on theory (Xue et al., 2021) is constructed in this paper, which effectively solves the above two problems.

 2.2 Energy flow-based heat network operation model

Based on the energy flow theory, utilizing an auxiliary variable instead of the nonlinear term (i.e., the product of temperature and mass flow rate) is considered in this paper, and defines the available heat quantity at the inlet and outlet of the pipe ij at period t under scenario s as h in ij,t,s and h out ij,t,s , respectively, as shown in the following equation: h i j , t , s in = c w m i j S T i j , t , s S , in − T i j , t , s R , in h i j , t , s out = c w m i j S T i j , t , s S , out − T i j , t , s R , out

Then Eqs. 4, 5 are transformed into the following equation: ∑ k i h k i , t , s out + H i , t , s N = ∑ i j h i j , t , s in

Due to the flow direction of heat pipeline in the UIES planning problem varies with the planning scheme, the energy flow-based heat network operation model in Xue et al. (2021) is unable to account for transmission delay characteristics of thermal mass. To address this issue, this paper depicts the flow direction of heat pipeline based on the spanning tree variables in the spanning tree theory, and rewrites the delay characteristic equation Eqs. 8, 9 into the following form: − 1 − b i j H M ≤ h i j , t , s in − h i j , t − τ i j , s out − h i j , t , s loss ≤ 1 − b i j H M − 1 − b j i H M ≤ h i j , t , s out − h i j , t − τ i j , s in − h i j , t , s loss ≤ 1 − b j i H M h i j , t , s loss = c w m i j S T i j , t , s S , in + T i j , t , s R , in − 2 T t , s am 1 − e − λ i j L i j m i j S A i j ρ w c w where h ij , t , s loss is the loss of available heat quantity of the pipe ij at period t under scenario s; b ij H and b ji H are binary spanning tree variables. If node i is the parent of node j in the heat network, b H ij = 1 and b H ji = 0; if node j is the parent of node i, b H ji = 1 and b H ij = 0. When b H ij = 1, the flow direction of pipeline ij is from node i to node j, therefore Eq. 17 transform to h i j , t , s in − h i j , t − τ i j , s out − h i j , t , s loss = 0 to depict the transmission delay characteristic of the available heat quantity of the pipe ij. Eq. 18 is the same as Eq. 17.

Based on Eqs 10–13, the upper and lower constraints on the available heat quantity are respectively established as shown in Eqs 20–23 and give the calculation of the upper and lower limits of the available heat quantity. h min ⁡ , i j in ≤ h i j , t , s in ≤ h max ⁡ , i j in h min ⁡ , i j out ≤ h i j , t , s out ≤ h max ⁡ , i j out h max ⁡ , i j in = c w m max ⁡ , i j S T max S − T min R h max ⁡ , i j out = c w m max ⁡ , i j S T max S − T min R h min ⁡ , i j in = 0 h min ⁡ , i j out = 0 where h max , ij in and h min , ij in are the upper and lower limits of the available heat quantity at the inlet of the pipe ij, respectively; h out max, ij and h out min, ij are the upper and lower limits of the available heat quantity at the outlet of the pipe ij, respectively; T S max and T S min are the maximum and minimum temperatures of water in supply network, respectively; T R max and T R min are the maximum and minimum temperatures of water in return network, respectively.

If the heat loss h ij , t , s loss is determined, the operation model of heat network based on energy flow theory is independent of temperature and mass flow rate and can be solved directly without considering the direction of water flow. The nonlinear term λ ij l ij / m i j S A i j ρ w c w in Eq. 19 is very close to 0, usually less than 0.001, so considering 1-e-x≈x, the loss of available heat quantity can be approximated as (Dorfner and Hamacher, 2014): h i j , t , s loss ≈ c w m i j S T i j , t , s S , in + T i j , t , s R , in − 2 T t , s am λ i j L i j m i j S A i j ρ w c w = λ i j L i j A i j ρ w T i j , t , s S , in + T i j , t , s R , in − 2 T t , s am .

Since the only controllable variable affecting the loss of available heat quantity is temperature, the temperature of pipe should be kept lower to minimize the heat loss. In calculating the loss of available heat quantity, the inlet temperature of pipe in Eq. 21 is set as its lower limit, therefore: h i j , t , s loss = λ i j L i j A i j ρ w T min ⁡ , i j S , in + T min ⁡ , i j R , in − 2 T t , s am where T S, in min, ij and T R, in min, ij are the lowest outlet temperatures of water flow in pipe ij of supply and return network, respectively.

 3 Bi-level planning model of station-network cooperative optimization of UIES 3.1 Objective of the upper level planning model

The objective of the upper level is to minimize the total cost C all, including the investment cost C inv and the operation cost C ope as shown in the following equation: min ⁡ C all = C inv + C ope

The investment cost C inv includes the device investment cost C dv and the pipeline investment cost C pl, as shown in the following equation: C inv = C line + C device C device = ∑ i ∈ N P ∑ d ∈ Ω D c inv d S i d R ir d C line = ∑ i j ∈ L ∑ d ∈ Ω L c inv d l i j y i j d R ir d R ir d = r 1 + r τ d 1 + r τ d − 1 where N P is the set of candidate nodes of energy station; L is the set of roads; ΩD is the set of the types of devices, including CHP, GB, electrical heater (EH), energy storage (ES), photovoltaic (PV) and wind turbine (WT); ΩL is the set of the types of pipelines, including electric lines and heat pipes; cd inv is the investment cost factor for the d-category device/pipeline; Rd ir is the payback factor for d-category device/pipeline, which spread the investment cost of the device over the planning period equally over the years of the life cycle; τ d is the life of d-category devices/pipelines; Sd i is the planning capacity of d-category device at node i; yd ij indicates whether d-category pipeline is installed between nodes i and j. If it is installed, yd ij,k = 1, otherwise, yd ij ,k = 0.

 3.2 Constraints of the upper level planning model 3.2.1 Constraints of planning for energy station

This paper considers that there are multiple candidate nodes of energy station, and defines the binary variable y i S to indicate whether energy station is planned at node i. If it is planned, y S i = 1; otherwise, y S i = 0. The number of energy stations in the UIES needs to meet the following constraints: N min S ≤ ∑ i y i S ≤ N max S where N S min and N S max are the lower and upper limits of the number of energy stations, respectively.

The CHP units, GB, EH, ES, PV and WT can be configured only when the node is planning an energy station, so there: y i S S min ⁡ , i d ≤ S i d ≤ y i S S max ⁡ , i d , d ∈ CHP , GB , EH , ES , WT , PV where Sd i denotes the installed capacity of d-category device at node i; Sd min,i and Sd max,i are the minimum and maximum installed capacity of d-category device at node i, respectively.

The total capacity of each type of equipment in the UIES needs to meet the following constraints: S min d ≤ ∑ i S i d ≤ S max d , d ∈ CHP , GB , EH , ES , WT , PV where Sd min and Sd max are the lower and upper limits of the total installed capacity of d-category equipment in the UIES.

 3.2.2 Constraints of planning for energy supply network

In this paper, the planning of electric lines and heat pipes is considered, in order to ensure the radial structure of the electric network, the following constraints are in place. Eqs 31, 32 and give the calculation of the number of parent node of a node based on the spanning tree theory. Eq. 36 indicates that if a candidate node is planned as an energy station, it must be the root node and has no parent node. Eq. 37 indicates that each load node has a parent node. Eq. 38 indicates that nodes other than candidate nodes and load nodes have at most one parent node. y i j E = b i j E + b j i E c i E = ∑ j i b j i E c i E ≤ 1 − y i S , i ∈ N P c i E = 1 , i ∈ N L c i E ≤ 1 , i ∈ N \ N P ∪ N L where N, N P and N L are the sets of all nodes, candidate energy station nodes and load nodes, respectively; y E ij indicates whether the electric line is planned between nodes i and j. If it is planned, y E ij = 1; otherwise, y E ij = 0; b E ij and b E ji are binary spanning tree variables. If node i is the parent of node j in the electric network, b E ij = 1 and b E ji = 0; if node j is the parent of node i, b E ji = 1 and b E ij = 0; c E i denotes the number of parent nodes of node i in the electric network.

Similarly, to ensure the radial structure of the heat network, the following constraints are imposed based on Eqs 34–38: y i j H = b i j H + b j i H c i H = ∑ j i b j i H c i H ≤ 1 − y i S , i ∈ N P c i H = 1 , i ∈ N L c i H ≤ 1 , i ∈ N \ N P ∪ N L where y H ij indicates whether the heat pipe is planned between nodes i and j. If it is planned, y ij H = 1; otherwise, y ij H = 0; c i H denotes the number of parent nodes of node i in the heat network.

 3.3 Objective of the lower level planning model

The objective of the lower level is to minimize the operating cost C ope, which includes purchasing electric cost C e buy , purchasing gas cost C g buy and maintenance cost C main, as shown in the following equation: min ⁡ C ope = C e buy + C g buy + C main C e buy = ∑ s ∑ t ∑ i α s sce P i , t , s buy c t E C g buy = ∑ s ∑ t ∑ i α s sce G i , t , s buy c t G C main = ∑ i ∈ N P ∑ d c m d S i d + ∑ i j ∈ L l i j y i j E c m E + y i j H c m H where P i , t , s buy is the power of purchasing electric at node i at period t under scenario s; G i , t , s buy is the power of purchasing gas at node i at period t under scenario s; P d i,t,s is the output of d-category device at node i at period t under scenario s; c m d is the unit maintenance cost of d-category device; c E m and c H m are the unit maintenance cost of electric line and heat pipe, respectively.

 3.4 Objective of the lower level planning model 3.4.1 Constraints of device

(1) Constraints of CHP.

The constraints of CHP are established as follows (Chen et al., 2023): H i , t , s CHP = K e CHP P i , t , s CHP H i , t , s CHP = K h CHP G i , t , s CHP P min CHP S i CHP ≤ P i , t , s CHP ≤ P max CHP S i CHP − r max CHP S i CHP Δ t ≤ P i , t , s CHP − P i , t − 1 , s CHP ≤ r max CHP S i CHP Δ t where P CHP i,t,s, G CHP i,t,s and H CHP i,t,s are the electric power, gas consumption power and heat power of the CHP at node i at period t under scenario s, respectively; K e CHP is the heat-to-electric ratio coefficient of the CHP; K h CHP is the gas-heat conversion efficiency of the CHP; P max CHP and P min CHP are the upper and lower limits of the generation power of the CHP, respectively; r max CHP is the maximum ramp rate of the CHP.

(2) Constraints of GB (Zhang et al., 2023).

H i , t , s GB = K GB G i , t , s GB H min GB S i GB ≤ H i , t , s GB ≤ H max GB S i GB where G i , t , s GB and H GB i,t,s are the gas consumption power and heat power of the GB at node i at period t under scenario s, respectively; K GB is the heat-gas conversion factor of the GB; H max GB and H min GB are the upper and lower limits of the heat power of the GB, respectively.

(3) Constraints of EH (Chen et al., 2023).

 H i , t , s EH = K EH P i , t , s EH H min EH S i EH ≤ H i , t , s EH ≤ H max EH S i EH where P EH i,t,s and H EH i,t,s are electric power and heat power of the EH at node i at period t under scenario s, respectively; K EH is the heat-electric conversion factor of the EH; HEH max and HEH min are the upper and lower limits of the heat power of the EH, respectively.

(4) Constraints of ES (Wang et al., 2021).

 x i , t , s ES , c + x i , t , s ES , f ≤ 1 0 ≤ P i , t , s ES , c ≤ x i , t , s ES , c S i ES 0 ≤ P i , t , s ES , f ≤ x i , t , s ES , f S i ES P i , t , s ES = P i , t , s ES , c − P i , t , s ES , f E i , t , s ES = E i , t − 1 , s ES + η ES , c P i , t , s ES , c − P i , t , s ES , f η ES , f 10 % E i ES , p ≤ E i , t , s ES ≤ 90 % E i ES , p where x ES,c i,t,s and x ES,f i,t,s are the charging and discharging state of the ES at node i at period t under scenario s, respectively; P ES,c i,t,s and P ES,f i,t,s are the charging and discharging power of the ES at node i at period t under scenario s, respectively; E ES i,t,s is the capacity of the ES at node i at period t under scenario s; η ES,c and η ES,f are the charging and discharging efficiencies of the ES, respectively; E i ES , p is the planning capacity of the ES at node i.

(5) Constraints of PV (Chen et al., 2023).

 0 ≤ P i , t , s PV ≤ α t , s PV , pre S i PV where α t , s PV , pre is the predicted output coefficients for PV at period t under scenario s; P PV i,t,s is the power outputs of PV at node i at period t under scenario s.

(6) Constraints of WT.

 0 ≤ P i , t , s WT ≤ α t , s WT , pre S i WT where α t , s WT , pre is the predicted output coefficients for WT at period t under scenario s; P i , t , s WT is the power outputs of WT at node i at period t under scenario s.

(7) Constraints of Power of Purchasing Gas.

 G i , t , s buy = G i , t , s CHP + G i , t , s GB

3.4.2 Constraints of electric network

The power flow model of the electric network is established based on the LinDistFlow model as follows (Baran and Wu, 1989; Su et al., 2023): P i , t , s N = P i , t , s buy + P i , t , s CHP + P i , t , s PV + P i , t , s WT + P i , t , s ES − P i , t , s EH − P i , t , s load Q i , t , s N = Q i , t , s buy + Q i , t , s CHP + Q i , t , s PV + Q i , t , s WT − Q i , t , s EH − Q i , t , s load ∑ k i ∈ L P k i , t , s L + P i , t , s N = ∑ i j ∈ L P i j , t , s L ∑ k i ∈ L Q k i , t , s L + Q i , t , s N = ∑ i j ∈ L Q i j , t , s L − M 1 − y i j E ≤ U i , t , s 2 − U j , t , s 2 − 2 P i j , t , s L R i j − 2 Q i j , t , s L X i j ≤ M 1 − y i j E U min 2 ≤ U i , t , s 2 ≤ U max 2 P i j , t , s L 2 + Q i j , t , s L 2 ≤ y i j E S i j 2 where P N i,t,s and Q N i,t,s are the active and reactive power at node i at period t under scenario s, respectively; P L ij,t,s and Q L ij,t,s are the active and reactive power transmitted by line ij at period t under scenario s, respectively; U 2 i,t,s is the squared value of the voltage at node i at period t under scenario s; R ij and X ij are the resistance and reactance of line ij, respectively; M comes from the big-M method and is a sufficiently positive number, so that for any unplanned line ij, U 2 i,t,s and U 2 j,t,s do not directly influence each other; U min 2 and U max 2 are the squares of the lower and upper limits of voltage, respectively; S ij is the capacity of line ij.

R ij and X ij in Eq. 67 require the substitutions shown in Eq. 70: R i j = y i j E l i j R X i j = y i j E l i j X where R and X are the resistance and reactance per unit length of the electric line, respectively.

Then Eq. 67 is transformed into Eq. 71: − M 1 − y i j E ≤ U i , t , s 2 − U j , t , s 2 − 2 P i j , t , s L y i j E l i j R − 2 Q i j , t , s L y i j E l i j X ≤ M 1 − y i j E

3.4.3 Constraints of heat network

In addition to the energy flow-based heat network operation model developed in Section 2, the constraint shown in Eq. 72 is included. H i , t , s N = H i , t , s CHP + H i , t , s GB + H i , t , s EH − H i , t , s load

In summary, the constraints of heat network include Eqs 16, 20, 21, 25, 72.

 4 Solution method of station-network cooperative optimization of UIES

Firstly, the nonlinear terms in the cooperative optimization model are linearized based on mathematical theory, and the bi-level nonlinear optimization model is transformed into a bi-level linear optimization model. Then based on the Lagrangian function and KKT condition, the lower level model is transformed into the additional constraints of upper level model, and the bi-level linear optimization model is transformed into a single-level linear optimization model, which can be efficiently solved by commercial solvers.

4.1 Linearization of the cooperative planning model

The nonlinear terms P L ij,t, syE ijl ij R and Q L ij,t, syE ijl ij R appear in Eq. 71, which need to be linearized to improve the efficiency of the solution. Define the auxiliary variable Z P ij,t,s = P L ij,t,s y E ij , Z Q ij,t,s = Q L ij,t,s y E ij , where y E ij is a 0–1 variable, so there are the constraints shown in Eqs 73–76. When y E ij = 1, Eqs 73, 74 are translated into Z P ij,k,t,s = P L ij,t,s ; when y E ij = 0, Eqs 73, 74 are translated into Z P ij,k,t,s = 0. Eqs 75, 76 are identical to the above. P i j , t , s L − M 1 − y i j E ≤ Z i j , t , s P ≤ P i j , t , s L + M 1 − y i j E − M y i j E ≤ Z i j , t , s P ≤ M y i j E Q i j , t , s L − M 1 − y i j E ≤ Z i j , t , s Q ≤ Q i j , t , s Q + M 1 − y i j E − M y i j E ≤ Z i j , t , s Q ≤ M y i j E

Then Eq. 71 is transformed into: − M 1 − y i j E ≤ U i , t , s 2 − U j , t , s 2 − 2 Z i j , t , s P l i j R − 2 Z i j , t , s Q l i j X ≤ M 1 − y i j E .

Further, Eq. 69 is linearized based on quadratic constraint linearization method, and two square constraints Eq. 78 are employed to substitute for Eq. 69 (Chen et al., 2016). − y i j E S i j ≤ P i j , t , s L ≤ y i j E S i j − y i j E S i j ≤ Q i j , t , s L ≤ y i j E S i j − 2 y i j E S i j ≤ P i j , t , s L + Q i j , t , s L ≤ 2 y i j E S i j − 2 y i j E S i j ≤ P i j , t , s L − Q i j , t , s L ≤ 2 y i j E S i j

At this point, the bi-level nonlinear optimization model is transformed into a bi-level linear optimization model.

 4.2 Transformation of the cooperative planning model based on KKT conditions

The Lagrangian function of the lower level model is shown in Eq. 79. L q u , q l , λ , μ = f q u , q l + ∑ i = 1 I λ i b i q u , q l + ∑ j = 1 J μ j d j q u , q l where f (q u , q l ) is the objective function of the lower level model; b (q u , q l ) and d (q u , q l ) are the set of inequality constraints and the set of equation constraints of the lower level model, respectively; q u and q l are the decision variables of the upper level model and lower level model, respectively; λ i ∈ λ and μ j ∈ μ are the Lagrange multipliers of the inequality constraints and equation constraints in the lower level model, respectively; I and J are the number of inequality constraints and equation constraints in the lower level model, respectively.

Based on Eq. 79, the Lagrangian function of the lower level model is constructed as in Eq. 80. L LM = ∑ s ∑ t ∑ i α s sce P i , t , s buy c t E + ∑ s ∑ t ∑ i α s sce G i , t , s buy c t G + ∑ i ∈ N P ∑ d c m d S i d + ∑ i j ∈ L l i j y i j E c m E + y i j H c m H + λ i , t , s 1 H i , t , s CHP − K e CHP P i , t , s CHP + λ i , t , s 2 H i , t , s CHP − K h CHP G i , t , s CHP + λ i , t , s 3 H i , t , s GB − K GB G i , t , s GB + λ i , t , s 4 H i , t , s EH − K EH P i , t , s EH + λ i , t , s 5 P i , t , s ES − P i , t , s ES , c + P i , t , s ES , f + λ i , t , s 6 E i , t , s ES − E i , t − 1 , s ES − η ES , c P i , t , s ES , c + P i , t , s ES , f / η ES , f + λ i , t , s 7 G i , t , s buy − G i , t , s CHP − G i , t , s GB + λ i , t , s 8 P i , t , s N − P i , t , s buy − P i , t , s CHP − P i , t , s PV − P i , t , s WT − P i , t , s ES + P i , t , s EH + P i , t , s load + λ i , t , s 9 Q i , t , s N − Q i , t , s buy − Q i , t , s CHP − Q i , t , s PV − Q i , t , s WT + Q i , t , s EH + Q i , t , s load + λ i , t , s 10 ∑ k i ∈ L P k i , t , s L + P i , t , s N − ∑ i j ∈ L P i j , t , s L + λ i , t , s 11 ∑ k i ∈ L Q k i , t , s L + Q i , t , s N − ∑ i j ∈ L Q i j , t , s L + λ i , t , s 12 ∑ k i h k i , t , s out + H i , t , s N − ∑ i j h i j , t , s in + λ i , t , s 13 H i , t , s N − H i , t , s CHP − H i , t , s GP − H i , t , s EH − H i , t , s load + u i , t , s 1 , ⁡ min P min CHP S i CHP − P i , t , s CHP + u i , t , s 1 , ⁡ max P i , t , s CHP − P max CHP S i CHP + u i , t , s 2 , ⁡ min − r max CHP S i CHP Δ t − P i , t , s CHP + P i , t − 1 , s CHP + u i , t , s 2 , ⁡ max P i , t , s CHP − P i , t − 1 , s CHP − r max CHP S i CHP Δ t + u i , t , s 3 , ⁡ min H min GB S i GB − H i , t , s GB + u i , t , s 3 , ⁡ max   H i , t , s GB − H max GB S i GB + u i , t , s 4 , ⁡ min H min EH S i EH − H i , t , s HP + u i , t , s 4 , ⁡ max H i , t , s EH − H max EH S i EH − u i , t , s 5 , ⁡ min x i , t , s ES , c + x i , t , s ES , f + u i , t , s 5 , ⁡ max x i , t , s ES , c + x i , t , s ES , f − 1 − u i , t , s 6 , ⁡ min P i , t , s ES , c + u i , t , s 6 , ⁡ max P i , t , s ES , c − x i , t , s ES , c S i ES − u i , t , s 7 , ⁡ min P i , t , s ES , f + u i , t , s 7 , ⁡ max P i , t , s ES , f − x i , t , s ES , f S i ES + u i , t , s 8 , ⁡ min 10 % E i ES , p − E i , t , s ES + u i , t , s 8 , ⁡ max E i , t , s ES − 90 % E i ES , p − u i , t , s 9 , ⁡ min P i , t , s PV + u i , t , s 9 , ⁡ max P i , t , s PV − α t , s PV , pre S i PV − u i , t , s 10 , ⁡ min P i , t , s WT + u i , t , s 10 , ⁡ max P i , t , s WT − α t , s WT , pre S i WT + u i j , t , s 11 , ⁡ min − M 1 − y i j E − U i , t , s 2 + U j , t , s 2 + 2 Z i j , t , s P l i j R + 2 Z i j , t , s Q l i j X + u i j , t , s 11 , ⁡ max U i , t , s 2 − U j , t , s 2 − 2 Z i j , t , s P l i j R − 2 Z i j , t , s Q l i j X − M 1 − y i j E + u i j , t , s 12 , ⁡ min P i j , t , s L − M 1 − y i j E − Z i j , t , s P + u i j , t , s 12 , ⁡ max Z i j , t , s P − P i j , t , s L − M 1 − y i j E + u i j , t , s 13 , ⁡ min Q i j , t , s L − M 1 − y i j E − Z i j , t , s Q + u i j , t , s 13 , ⁡ max Z i j , t , s Q − Q i j , t , s Q − M 1 − y i j E + u i j , t , s 14 , ⁡ min − M y i j E − Z i j , t , s P + u i j , t , s 14 , ⁡ max Z i j , t , s P − M y i j E + u i j , t , s 15 , ⁡ min − M y i j E − Z i j , t , s Q + u i j , t , s 15 , ⁡ max Z i j , t , s Q − M y i j E + u i j , t , s 16 , ⁡ min − y i j E S i j − P i j , t , s L + u i j , t , s 16 , ⁡ max P i j , t , s L − y i j E S i j + u i j , t , s 17 , ⁡ min − y i j E S i j − Q i j , t , s L + u i j , t , s 17 , ⁡ max Q i j , t , s L − y i j E S i j + u i j , t , s 18 , ⁡ min − 2 y i j E S i j − P i j , t , s L − Q i j , t , s L + u i j , t , s 18 , ⁡ max P i j , t , s L + Q i j , t , s L − 2 y i j E S i j + u i j , t , s 19 , ⁡ min − 2 y i j E S i j − P i j , t , s L + Q i j , t , s L + u i j , t , s 19 , ⁡ max P i j , t , s L − Q i j , t , s L − 2 y i j E S i j + u i j , t , s 20 , ⁡ min − 1 − b i j H M − h i j , t , s in + h i j , t − τ i j , s out + h i j , t , s loss + u i j , t , s 20 , ⁡ max h i j , t , s in − h i j , t − τ i j , s out − h i j , t , s loss − 1 − b i j H M + u i j , t , s 21 , ⁡ min − 1 − b j i H M − h i j , t , s out + h i j , t − τ i j , s in + h i j , t , s loss + u i j , t , s 21 , ⁡ max h i j , t , s out − h i j , t − τ i j , s in − h i j , t , s loss − 1 − b j i H M + u i j , t , s 22 , ⁡ min h min ⁡ , i j in − h i j , t , s in + u i j , t , s 22 , ⁡ max h i j , t , s in − h max ⁡ , i j in + u i j , t , s 23 , ⁡ min h min ⁡ , i j out − h i j , t , s out + u i j , t , s 23 , ⁡ max h i j , t , s out − h max ⁡ , i j out

According to the constructed Lagrangian function and KKT complementary relaxation conditions, the lower level model can be transformed into the additional constraints of the upper level model, thus transforming the bi-level linear optimization model into a single-level linear optimization model, and the simplified form of the model is as follows: min q u ∈ Q U , q l ∈ Q L F q u , q l B q u , q l ≤ 0 , D q u , q l = 0 b q u , q l ≤ 0 , d q u , q l = 0 ∇ q l L q u , q l , λ , μ = 0 λ i ≥ 0 , i = 1 , 2 , … , I λ i b i q u , q l = 0 , i = 1 , 2 , … , I where F (q u , q l ) is the objective function of the upper level model; B (q u , q l ) and D (q u , q l ) are the set of inequality constraints and the set of equation constraints of the upper level model, respectively.

The objective function corresponding to F (q u , q l ) is shown in Eq. 87. B (q u , q l ) includes the constraints in Eq. 88. D (q u , q l ) includes the constraints in Eq. 89. b (q u , q l ) includes the constraints in Eq. 90. d (q u , q l ) includes the constraints in Eq. 91. The constraint shown in Eq. 84 can be easily obtained by calculating the partial derivatives of the decision variables of the lower level model in Eq. 80. min ⁡ C all = C inv + C ope Equations  31 − 33 , 36 , 38 , 41 , Equations  34 − 35 , 37 , 39 − 40 , Equations  50 − 51 , 53 , 55 , 56 − 58 , 61 − 63 , 68 , 73 − 78 Equations  48 − 49 , 52 , 54 , 59 , 60 , 64 − 66 ,

At this point, the cooperative optimization model is transformed into a single-level linear optimization model that can be efficiently solved by the mathematical solvers.

 5 Case studies

In this section, the proposed method is demonstrated on a 55-node and 77-road urban topology. A computer with an Intel i9-13900HX CPU and 32 GB memory is used, and the MILP problems are modelled in MATLAB R2021a with the YALMIP package and solved by GUROBI 10.0.0 with the parameter MIPGap set as 0.01%.

5.1 Test system

The urban road topology shown in Figure 1 is used to test the proposed method, containing 55 nodes and 77 roads. There are six load centers in the region, and the maximum electric and heat load demand during the planning period is shown in Table 1. Typical operation scenarios are shown in Figures 2–5. The technical parameters of the proposed planning device are shown in Table 2.

Urban road topology.

Maximum electrical and heat loads during the planning period.

Node Electrical load (kW) Heat load (kW)
10 4180 4130
13 1720 2660
32 860 2210
45 4670 4330
49 1430 2380

The curve of electrical load for typical operating scenarios.

The curve of heat load for typical operating scenarios.

The curve of PV for typical operating scenarios.

The curve of WT for typical operating scenarios.

Parameters of devices to be planned.

Device Parameters Value Device Parameters Value
CHP c i n v CHP 7000 yuan/kW GB c i n v GB 2450 yuan/kW
c m CHP 0.05 yuan/kW c m GB 0.03 yuan/kW
τ CHP 25 years τ GB 25 years
K h CHP / K e CHP 0.3/0.45 K h GB 0.95
EH c i n v EH 2800 yuan/kW ES c i n v ES 1960 yuan/kWh, 560 yuan/kW
c m EH 0.04 yuan/kW c m ES 0.043 yuan/kW
τ EH 25 years τ ES 12 years
K e EH 0.95 η ES,c/η ES,f 0.95
PV c i n v PV 5000 yuan/kW WT c i n v WT 7500 yuan/kW
c m PV 0.04 yuan/kW c m WT 0.05 yuan/kW
τ PV 30 years τ WT 30 years
5.2 Analysis of the results of the proposed method

Figure 6 shows the results of energy station siting and energy supply network planning, Table 3 shows the results of device capacity allocation in the energy station, and Table 4 shows the various costs of the UIES.

Results of energy station siting and energy supply network planning.

Results of device capacity allocation in energy station.

CHP (MW) GB (MW) EH (MW) PV (MW) WT (MW) ES
4.35 7.33 2.23 12.60 10.00 2.40 MW/4.80 MWh

Various costs of UIES.

C inv C ope
C dv (104 yuan) C pl (104 yuan) C e buy (104 yuan) C g buy (104 yuan) C main (104 yuan)
1406.6 1399.9 76.7 1322.9 162.3

As can be seen in Figure 1, the proposed method chooses to construct energy station at node 14 and constructs 15 power supply lines/heat pipes. The planned energy supply network does not have a ring network and satisfies the requirement of radial constraint. In order to further verify the effectiveness of the proposed method, the planning scheme and the total system cost are solved one by one when the energy station is constructed at other nodes, and the result reveals that the total cost of the planning scheme of the proposed method is one of the smallest. In summary, the proposed method can effectively plan the energy supply network.

As can be seen from Tables 3, 4, the proposed method in configures larger capacity PV and WT, which in conjunction with CHP can basically satisfy the electric load demand during the planning period, so the cost of purchasing electric in UIES is smaller. In addition, the proposed method configures a larger capacity of GB and a smaller capacity of EH, which means that the heat load is mainly satisfied by GB consuming gas, and is assisted by EH when more electric is generated in the region. This shows that the proposed method can reasonably plan the capacity of energy devices including CHP, GB, EH, PV, WT, and ES.

 5.3 Analysis of the efficiency of the proposed heat network operation model

In order to analyze the efficiency of the proposed heat network operation model, the following scenario is set up as a comparison with the proposed method.

Case 1–1: Replace of the heat network operation model of the proposed method with the energy flow-based heat network operation model from Xue et al. (2021).

Table 6 demonstrates the comparison of UIES’s costs of the proposed method and Case 1–1. The results of energy station siting and energy supply network planning of Case 1–1 are the same as Figure 6. Figures 7, 8 show the total heat power of the heat equipment in the system of the proposed method and Case 1–1.

Total heat power of the heat equipment in the system of the proposed method.

Total heat power of the heat equipment in the system of Case 1–1.

The comparison in Table 5 shows that compared to Case 1–1, the planning cost of proposed method is reduced by 5.79% and the operating cost is reduced by 3.98%, which in turn lead to a 5.15% reduction in the total cost. As can be seen from the comparison of Figures 7, 8, this is because considering the transmission delay can utilize the virtual heat storage characteristics of the heat network, which in turn smooths out the fluctuations of the heat power of the system and reduces the maximum power of heat load, thus reducing the required capacity of the heat equipment and lowering the investment cost and operation cost.

Comparison of UIES’s costs of the proposed method and Case 1–1.

C all (104 yuan) C inv (104 yuan) C ope (104 yuan)
Proposed method 4368.4 2806.5 1561.8
Case 1–1 4605.5 2978.9 1626.6
5.4 Analysis of the impact of considering different devices on planning results

In order to analyze the impact of considering different devices on planning results, the following scenarios are set up as a comparison with the proposed method. Table 6 shows the comparison of the planning results.

Comparison of system costs in five scenarios.

Cost (104 yuan) Proposed method Case 2–1 Case 2–2 Case 2–3 Case 2–4
C dv 1406.6 1081.8 1202.4 1587.2 1266.7
C pl 1399.9 1399.9 1399.9 1399.9 1399.9
C buy e 76.7 376.1 362.8 376.3 46.6
C buy g 1322.9 2252.9 1419.0 1140.1 1633.1
C main 162.3 162.2 162.2 162.3 162.2
C all 4368.4 5272.9 4546.5 4665.8 4508.6

Case 2–1: Planning for PV and WT is not considered.

Case 2–2: Planning for ES is not considered.

Case 2–3: Planning for GB is not considered.

Case 2–4: Planning for EH is not considered.

The comparison between the proposed method and Case 2–1 shows that in Case 2–1, the device investment cost for UIES is reduced because the planning for PV and WT is not considered, but UIES not only needs to purchase more electric from the superior grid, but also needs to purchase more natural gas for CHP generation, so the cost of purchasing electric and gas is significantly increased. Thus, in Case 2–1, although the investment cost for UIES is reduced by 11.6%, the total cost is increased by 20.7%.

The comparison between the proposed method and Case 2–2 shows that in Case 2–2, the device investment cost for UIES is reduced because the planning for ES is not considered, but the purchasing electric cost for UIES is significantly increased due to the absence of ES to provide peak shaving and valley filling. Thus, in Case 2–2, although the investment cost for UIES is reduced by 7.3%, the total cost is increased by 4.1%.

The comparison between the proposed method and Case 2–3 shows that in Case 2–3, not considering the planning for GB makes the purchasing gas cost for UIES decrease significantly, but the capacity of CHP and EH is increased in order to meet the heat load demand, which leads to an increase of the device investment cost and purchasing electric cost for UIES on the contrary, and consequently the total cost of the system increases by 6.8%.

The comparison between the proposed method and Case 2–4 shows that in Case 2–4, the purchasing electric cost for UIES is reduced due to the fact that EH is not considered in the planning, but the capacity of CHP and GB is increased in order to meet the heat load demand, which leads to an increase in the purchasing gas cost for UIES instead, and consequently the total cost of the system is increased by 3.2%.

In summary, it can be seen that the consideration of PV, WT, ES, EH, GB in station-network cooperative planning of UIES are all conducive to reducing the total cost of the system and improving the economics of the planning scheme.

 6 Conclusion

In this paper, a station-network cooperative planning method of UIES based on energy flow model is proposed to realize the cooperative optimization of the siting and sizing of energy station and the topology of energy supply network. And the following conclusions are obtained through arithmetic example analysis: 1) The proposed method can effectively plan the energy supply network and reasonably configure the capacity of energy devices including CHP, GB, EH, PV, WT and ES.

2) Compared with the existing energy flow-based heat network operation model, the proposed method takes into account the transmission delay of the heat network, which can more fully utilize the heat storage capacity and thus improve the economic efficiency of the planning scheme.

3) The consideration of PV, WT, ES, EH, GB in station-network cooperative planning of UIES are all conducive to reducing the total cost of the system and improving the economics of the planning scheme.

 Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

 Author contributions

XF: Conceptualization, Methodology, Project administration, Resources, Software, Writing–original draft. SY: Conceptualization, Methodology, Software, Validation, Writing–original draft, Writing–review and editing. QH: Funding acquisition, Methodology, Software, Writing–original draft. LW: Data curation, Funding acquisition, Investigation, Visualization, Writing–original draft. CN: Investigation, Visualization, Writing–original draft. LY: Conceptualization, Project administration, Resources, Supervision, Validation, Writing–review and editing.

 Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was funded by the Science and Technology Project of State Grid Jiangsu Electric Power Co., Ltd. (No. J2023083). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

 Conflict of interest

Authors XF, QH, and LW were employed by the State Grid Wuxi Power Supply Company of Jiangsu Electric Power Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 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.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2024.1363822/full#supplementary-material

References Baran M. Wu F. (1989). Network reconfiguration in distribution systems for loss reduction and load balancing. IEEE Trans. Power Deliv. 4 (2), 1401–1407. 10.1109/61.25627 Chen C. Li Y. Qiu W. Liu C. Zhang Q. Li Z. (2022). Cooperative-game-based Day-Ahead scheduling of local integrated energy systems with Shared energy storage. IEEE Trans. Sustain. Energy 13 (4), 1994–2011. 10.1109/TSTE.2022.3176613 Chen C. Zhang T. Wang Y. Sun G. Wu C. Lin Z. (2024). Coordinated siting and sizing for integrated energy system considering generalized energy storage. Int. J. Electr. Power & Energy Syst. 155, 109619. 10.1016/j.ijepes.2023.109619 Chen C. Zhu Y. Zhang T. Li Q. Li Z. Liang H. (2023a). Two-stage multiple cooperative Games-based Joint planning for Shared energy storage provider and local integrated energy systems. Energy 284, 129114. 10.1016/j.energy.2023.129114 Chen S. Su W. Wu B. (2023b). Two stage robust planning of park integrated energy system considering low carbon. Front. Energy Res. 10, 1100089. 10.3389/fevo.2022.1100089 Chen X. Wu W. Zhang B. (2016). Robust Restoration method for active distribution networks. IEEE Trans. Power Syst. 31 (5), 4005–4015. 10.1109/TPWRS.2015.2503426 Chen Y. Lao W. Qi D. Hui H. Yang S. Yan Y. (2023c). Distributed Self-triggered Control for Frequency Restoration and active power Sharing in Islanded Microgrids. IEEE Trans. Industrial Inf. 19 (10), 10635–10646. 10.1109/TII.2023.3240738 Dong W. Lu Z. He L. Zhang J. Ma T. Cao X. (2023). Optimal expansion planning model for integrated energy system considering integrated demand Response and Bidirectional energy Exchange. CSEE J. Power Energy Syst. 9 (4), 1449–1459. 10.17775/CSEEJPES.2021.09220 Dorfner J. Hamacher T. (2014). Large-scale District heating network optimization. IEEE Trans. Smart Grid 5 (4), 1884–1891. 10.1109/TSG.2013.2295856 Du Y. Xue Y. Wu W. Shahidehpour M. Shen X. Wang B. (2024). Coordinated planning of integrated electric and heating system considering the optimal reconfiguration of District heating network. IEEE Trans. Power Syst. 39 (1), 794–808. 10.1109/TPWRS.2023.3242652 Farrokhifar M. Nie Y. Pozo D. (2020). Energy systems planning: a Survey on models for integrated power and natural gas networks coordination. Appl. Energy 262, 114567. 10.1016/j.apenergy.2020.114567 Ha T. Xue Y. Lin K. Zhang Y. Thang V. V. Nguyen T. (2022). Optimal operation of energy hub based Micro-energy network with integration of renewables and energy storages. J. Mod. Power Syst. Clean Energy 10 (1), 100–108. 10.35833/MPCE.2020.000186 Hong B. Chen J. Zhang W. Shi Z. Li J. Miao W. (2018). Integrated energy system planning at Modular Regional-User level based on A two-layer bus structure. CSEE J. Power Energy Syst. 4 (No.2), 188–196. 10.17775/CSEEJPES.2018.00110 Li C. Li P. Yu H. Li H. Zhao J. Li S. (2021). Optimal planning of Community integrated energy station considering Frequency regulation Service. J. Mod. Power Syst. Clean Energy 9 (2), 264–273. 10.35833/MPCE.2019.000056 Li T. Han X. Wu W. Sun H. (2023). Robust expansion planning and Hardening strategy of meshed multi-energy distribution networks for resilience Enhancement. Appl. Energy 341, 121066. 10.1016/j.apenergy.2023.121066 Liu X. (2020a). Energy stations and pipe network collaborative planning of integrated energy system based on load complementary characteristics. Sustain. Energy Grids Netw. 23, 100374. 10.1016/j.segan.2020.100374 Liu X. (2020b). Pipeline network layout design of integrated energy system based on energy station site selection and load complementary characteristics. IEEE Access 8, 1–92790. 10.1109/ACCESS.2020.2991239 Mu Y. Chen W. Yu X. Jia H. Hou K. Wang C. (2020). A double-layer planning method for integrated community energy systems with varying energy conversion efficiencies. Appl. Energy 279, 115700. 10.1016/j.apenergy.2020.115700 Strezoski L. Padullaparti H. Ding F. Baggu M. (2022). Integration of utility distributed energy resource management system and aggregators for evolving distribution system operators. J. Mod. Power Syst. Clean Energy 10 (2), 277–285. 10.35833/MPCE.2021.000667 Su Y. Liu F. Kang K. Wang Z. (2023). Implicit sparsity of LinDistFlow model. IEEE Trans. Power Syst. 38 (5), 4966–4969. 10.1109/TPWRS.2023.3291214 Wang J. Du W. Yang D. (2021). Integrated energy system planning based on life cycle and emergy theory. Front. Energy Res. 9, 713245. 10.3389/fenrg.2021.713245 Wang Y. Li R. Dong H. Ma Y. Yang J. Zhang F. (2019). Capacity planning and optimization of business park-level integrated energy system based on investment constraints. Energy 189, 116345. 10.1016/j.energy.2019.116345 Wu C. Gu W. Jiang P. Li Z. Cai H. Li B. (2018). Combined economic dispatch considering the time-delay of district heating network and multi-regional indoor temperature control. IEEE Trans. Sustain. Energy 9 (1), 118–127. 10.1109/TSTE.2017.2718031 Xiao H. Pei W. Dong Z. Kong L. (2018). Bi-Level planning for integrated energy systems incorporating demand response and energy storage under uncertain environments using novel metamodel. CSEE J. Power Energy Syst. 4 (2), 155–167. 10.17775/CSEEJPES.2017.01260 Xue Y. Shahidehpour M. Pan Z. Wang B. Zhou Q. Guo Q. (2021). Reconfiguration of district heating network for operational flexibility enhancement in power system unit commitment. IEEE Trans. Sustain. Energy 12 (2), 1161–1173. 10.1109/TSTE.2020.3036887 Zhang X. Shahidehpour M. Alabdulwahab A. Abusorrah A. (2015). Optimal expansion planning of energy hub with multiple energy infrastructures. IEEE Trans. Smart Grid 6 (5), 2302–2311. 10.1109/TSG.2015.2390640 Zhang Y. Li J. Ji X. Ye P. Yu D. Zhang B. (2023). Optimal Dispatching of electric-heat-Hydrogen integrated energy system based on Stackelberg game. Energy Convers. Econ. 4, 267–275. 10.1049/enc2.12094