Optimal Scheduling of Electricity–Gas–Heating System Based on Improved Alternating Direction Method of Multipliers

: With the joint optimization of the electricity–gas–heating system (EGHS) attracting more and more attention, a distributed optimized scheduling framework for EGHS based on an improved alternating direction method of multipliers (ADMM) algorithm is put forward in this paper. The framework of the proposed algorithm is a co-ordinated convex distribution framework with inner and outer layers. The outer layer is a penalty convex–concave procedure (PCCP), the inner layer is an ADMM-FE (forced equality) procedure. In this framework, the outer layer optimization uses the convex and concave procedure to turn the non-convex airﬂow equation into a second-order cone constraint with successive iterations, and the inner layer ADMM-FE algorithm solves the convex model to obtain a convergent solution. In the end, we compare the established algorithm with the traditional ADMM algorithm and the centralized optimization algorithm through example simulation analysis, and the results verify the e ﬀ ectiveness of the proposed model and optimization algorithm framework.


Introduction
With the further advancement of the industrialization reform of the energy industry, the energy internet has received increasing attention and research. Integrated energy systems (IES) that include gas units, gas to power (P2G) equipment, and combined heat and power (CHP) equipment break the barriers to exchange between different forms of energy. Therefore, the integrated energy system can fully tap the potential of each energy to complement each other efficiently, and meet various requirements of economic and environmental protection [1][2][3][4]. In addition, in the open electricity market environment, the electricity supplier, gas supplier, and heating supplier in the electric heat interconnection system generally belong to different suppliers, it is important to ensure the privacy of the information in the co-ordinated operation, so data needs to be protected [5,6].
Many research works are reported, focusing on optimization technologies of IES which combine two types of energy flow. Some researchers established an electricity-gas system considering the influence of temperature, and convexified the power grid [7]. Finally, the optimal power flow was studied on this model. A robust co-optimization scheduling model was proposed to study the co-ordinated optimal operation of electricity-gas systems, and the optimization problem of the model was solved by the alternating direction method of multipliers [8]. The research proposed the piecewise linear approximation methods to solve the non-convex and non-linear Weymouth equation in gas networks [9]. However, the above studies are mostly based on a single electricity-gas system or system or electricity-heating system. There are few studies that concentrate on the electricity-gasheating system and most are lacking consideration of network dynamic characteristics. Based on this, this paper builds a day-ahead scheduling optimization model of the electricity-gas-heating system (EGHS), considering the dynamic characteristics of the network in order to tap the potential of complementary advantages between different energy flows.
On the other hand, in the existing research, the optimization methods of the IES system mostly adopt the centralized optimization method. A convex relaxation approach was proposed to achieve the co-ordinated operation of electricity and natural gas systems [10]. Some researchers used a genetic algorithm to achieve the co-ordinated optimization of the district electricity and heating system [11]. However, the centralized solution used in the above research requires an upper-level co-ordination center to collect the operating information of each subsystem to solve the optimal strategy, which ignores the actual needs of each subsystem to protect its privacy. Distributed computing methods can solve this problem. The distributed algorithm mainly includes the alternating direction multiplier method [12][13][14] and the Benders decomposition algorithm [15][16][17]. This kind of distributed algorithm is only for the 2-Block problem, meaning it is only used to solve the electricity-gas system or the electricity-heating system. The EGHS contains three subsystems, so it belongs to the 3-Block problem in distributed optimization. Therefore, conventional distributed algorithms cannot ensure the convergence of calculation results [18].
Aiming at the above problems, this paper builds a day-ahead scheduling model of the EGHS considering the dynamic characteristics of the network, then uses inner and outer algorithms to solve the model. The outer penalty convex-concave procedure (PCCP) is responsible for the model convexification and the inner alternating direction method of multipliers (ADMM)-FE (forced equality) gets a convergent solution.
The rest of the paper is structured as follows: Section 2 establishes a day-ahead scheduling model of the EGHS. Section 3 introduces the PCCP-ADMM-FE algorithm for the day-ahead scheduling. Numeral testing results are provided in Section 4. Finally, conclusions are summarized in Section 5.

The Structure of EGHS
The EGHS studied in this paper consists of a power network, a natural gas network, and a heating network. Among them, the grid converts electric energy to natural gas through P2G equipment. The gas network converts natural gas to heat energy and electric energy through CHP equipment, and the gas network also converts natural gas to electric energy through the gas turbine. The coupling relationship between the systems is shown in Figure 1.

Power Network Model
In order to simplify the power network model and ensure the convexity of the model, this paper uses the DC power flow model to simulate the power network constraints: In the above, Ω EN is the grid node set; A G , A P2G , and A CHP are node-unit association matrix, node-P2G association matrix, and node-CHP association matrix; P D,t is the node load matrix; B is the imaginary part of the node admittance matrix; θ t is the node voltage phase angle vector; P max 1,ij and x ij are the maximum transmission power and reactance of the line respectively; P G,t is the unit active output vector; P P2G,t is the P2G device consumes active power vector; P CHP,t is the CHP device active input power vector; a d and a u are the upper and lower limits of the unit climbing rate constraint;θ t slack is the balance node phase angle.

Natural Gas Network Model
The transmission speed of natural gas in pipelines is related to various factors such as pressure and temperature. The natural gas steady-state model ignores the transmission delay of natural gas in the pipeline, and simply equates the pipeline injection gas flow with the outflow gas flow, resulting in an increase in calculation error. This paper uses gas pipeline partial differential equations to describe the dynamic characteristics of natural gas networks: In the above, b t ij,x , f t ij,x are the gas pressure and gas flow at pipe x at time t, respectively. M 1 and M 2 are pipe transmission constants and the specific expression of M1, M2 are provided in Appendix B; Equations (2)-(3) are partial differential equations, which are difficult to solve directly. Therefore, they are transformed into algebraic equations by Euler's finite difference method. The pipe ij is divided into N segments, each segment length is ∆x ij = l ij /N, and the equation is expressed as follows: In the above, ∆t is the time step and ∆x ij is the space step. There is a delay effect in the transmission of natural gas in the pipeline, and some natural gas is retained in the pipeline, which is called "pipe storage". After completing a scheduling period, the storage should meet the following constraint: In the above, Ω GP is a natural gas pipeline collection. The pressurization station must meet the following constraints: In the above, b t i , b t ij,0 is the pressure before and after boosting; k min c , k max c are the lower and upper limits of the career ratio; f t i is the pipe airflow and f max in is the pipe airflow limit. Node air pressure and air source output constraints are as follows: In the above, Ω g and Ω GB are gas source sets and gas node sets, respectively. Gas tank constraints are as follows: In the above, Ω s is the collection of gas storage tanks; Q in S,n,t and Q out S,n,t are the natural gas input and output of the gas storage tank n at time t, respectively; µ in n , µ out n are the inflation and deflation efficiency of the gas storage tank n at time t; S cap n is the rated gas storage capacity of the gas storage tank n.
The node traffic balancing constraint is as follows: In the above, B g , B P2G , B S , B CHP , B GT , and A g are node-air source correlation matrix, node-P2G correlation matrix, node-gas tank association matrix, node-CHP correlation matrix, node-gas unit association matrix, and node-Pipeline association matrix, respectively. f g,t , f P2G,t , f CHP,t , f GT,t , and f D,t are the gas source output vector, P2G natural gas injection vector, CHP natural gas injection vector, gas unit gas injection vector, and natural gas load vector, respectively.
The dynamic model of the natural gas network can accurately describe the transmission characteristics of natural gas pipelines. The natural gas network is coupled to the power network and the heating network through a gas turbine, CHP. Accurately establishing the natural gas network model can fully explore the potential of EGHS's energy complementarity.

Heating Network Model
In the heating network, the water flow is mostly used as the heat medium. The dynamic characteristics of the heat network are mainly reflected in the heat transfer delay and temperature loss. Therefore, in the heating network model, the water flow rate is regarded as constant, and the heating network dynamic characteristics are described as follows: In the above, T x k,t and M x k,t are the water temperature and water flow at the distance from the nozzle x at time t of the k-pipe. µ k is the heat loss factor, and c w is the specific heat capacity of the water. ρ k is the density of water and R k is the radius of the pipe. T w k,t is the ambient temperature of the k pipe at time t.
Solving Equation (11): In the above, T in k,t and T out k,t+∆τ k , are the water inflow and outflow temperatures of the kth pipe at time t. ∆τ k is the time it takes for water to pass through the k-pipe.
The water flow can be regarded as a constant, so the above construction can be converted into the following form: In the above, L k is the length of the k pipe and M k is the water flow of the k pipe. The electrical and heating output relationship of the CHP unit can be expressed as follows: In the above, k is the water flow. T s c,t is the water supply temperature. T r c,t is the return water temperature. Heat exchange station constraints are as follows: Heating network node water temperature constraints: In the above, Ω pipe,out,k and Ω pipe,in,k are the collections of pipes with node k at the outflow point and the inflow point, respectively. T mix,k is the mixing temperature of node k.
The water temperature of the water supply network and the backwater network meet the following constraints: In the above, T max s and T min s are the upper and lower limits of the water temperature of the water supply network, respectively. T max r and T min r are the upper and lower limits of the water temperature of the backwater network.
The dynamic model of the heating network can be used to describe the operating characteristics more accurately. The heating energy stored in the pipe can be used as the energy buffer of EGHS, which increases the flexibility of scheduling [19].

Coupling Device Constraint
Power networks and natural gas networks are coupled with gas turbines and P2G devices and need to meet the following constraints: In the above, Ω GT , Ω P2G are gas turbines and P2G sets, respectively. P GT,i,t , f GT,i,t are the active power of the gas turbine and the corresponding natural gas consumption at time t. P P2G,j,t , f P2G,j,t are the active power of P2G device and the corresponding natural gas conversion amount at time t. h 2,i , h 1,i , h 0,i , and h P2G,i are coefficients.
In addition, the natural gas network is coupled to the heating network and the power network through the CHP device. The working mode of the CHP is to determine the electrical output through the heat output, and the constraints are as follows: In the above, Ω CHP is a collection of CHP. H CHP,c,t , P CHP,c,t are CHP heating power and electric power respectively. K CHP,c is the ratio of CHP output power and thermal energy. f CHP,c,t is the natural gas consumption of CHP unit. η CHP is CHP unit conversion efficiency. H GV is the calorific value of natural gas.

Objective Function
The objective function is the total operation cost: In the above, T is the scheduling period. Ω G is the collection of coal-fired units. Ω g is the collection of wind turbines. Ω W is the collection of wind turbines. Ω s is the collection of gas cans. Ω h is the collection of heat sources. P G,i,t is the output of coal-fired unit i at time t. a i , b i , and c i are the unit cost coefficients. f g,j,t is the gas output of the gas source j at time t. ξ j is the cost coefficient of purchasing Appl. Sci. 2020, 10, 1214 7 of 18 gas. P W,k,max is the upper limit of the wind turbine output. P W,i,t is the output of wind turbine i at time t. δ k is the penalty coefficient for discarding wind energy. C S,n,t is the cost coefficient of the gas storage tank. H h,m,t is the output of heat source m at time t. C h,m,t is the cost of purchasing heating energy.
The electric-gas-heating system's daily dispatching model aims at minimizing the operating cost of the system, that is, minimizing the system's energy supply cost through optimized scheduling.

Distributed Optimized Scheduling of EGHS
In order to protect the privacy of each system operation, this paper uses an improved ADMM algorithm to solve the day-ahead scheduling problem. The application of the ADMM algorithm requires that the model is convex, while Equations (5) and (21) are non-convex constraints. Therefore, this paper proposes a convex distributed optimization framework for internal and external co-ordination.

Outer Layer PCCP Algorithm
The constraint in Equation (5) is processed as follows [20]: In the above, Equation (27) is a convex constraint. The relaxation of Equation (28) is as follows: In the above, b t,(k−1) ij,d+1 and f t,(k−1) ij,d+1 are the pressure value and airflow of the pipeline position d + 1 at time t after the k − 1 iteration; δ t ij,d is a relaxation variable with a value greater than or equal to 0. The constraint in Equation (21) is processed as follows: Since unnecessary gas consumption will increase the cost of energy supply to the system, the constraint in Equation (29) is always tight. In order to keep the relaxation domain tight, the relaxation variables are added to the inner optimization objective function: In the above, W (k) is the objective function value including the penalty function obtained by the kth optimization; ρ (k) is the penalty factor for the kth iteration of the outer PCCP; Ω d is the set of relaxation variables.
The convergence conditions of PCCP are as follows: In the above, ε 1 and ε 2 are convergence determination parameters; δ is the relaxation variable matrix; If the convergence conditions are not met, the penalty factor is updated as follows: In the above, ρ max is the maximum penalty factor; v c is the dynamic adjustment coefficient.

Inner Layer ADMM-FE Algorithm
Generally, ADMM is only used to solve two-operator problems. The electricity-gas-heating system model constructed in this paper is a three-operator problem. If the ADMM is directly extended by the three operators, the calculation result does not necessarily converge. The electricity-gas-heating system optimization problem is as follows: In the above, x, y, and z are separable three operators. A, B, and C are the constraint coefficient vectors related to the grid, gas network, and heat grid, respectively. D is the parameter vector of the correlation equation.
Augmented Lagrangian function: In the above, L 3 β (x, y, z, λ) is an augmented Lagrangian function that introduces Lagrangian multiplier λ and penalty factor η.
The standard ADMM iteration format is as follows: The above formula is a standard ADMM iterative format, and its convergence is guaranteed, but y k+1 and z k+1 are updated synchronously, and distributed computing is not implemented. In order to achieve fully distributed computing, ADMM is directly promoted as follow: Comparing Equation (35) with Equation (36), Equation (35) uses (Ax k+1 , By k , Cz k ,λ k ) to update y k+1 , but Equation (29) uses (Ax k+1 , By k , Cz k+1 , λ k ) to update y k+1 . The equal status of y and z is destroyed, resulting in the algorithm may not converge. Therefore, modify the iteration format as follows: The above formula updates y k+1 and z k+1 with the same information, guarantees the equal and parallel processing of y and z, and adds regular items when updating y and z. The meaning is to add a self-constraining mechanism to the variable to correct the variables. The improved ADMM algorithm guarantees the convergence of the algorithm by forcing equal updates of the variables. The improved algorithm is ADMM-FE.

PCCP-ADMM-FE Algorithm Steps
Based on the above algorithms, a convex distributed optimization framework for internal and external co-ordination is established in this paper. The specific steps are as follows: (1) Initialize the parameters of the EGHS system, ignore Equation (28), and use the second-order cone program model to obtain a good initial iteration point. Initialize PCCP parameters and enter the inner layer calculation.
(2) Relax the equality constraint to the objective function, then yield an augmented Lagrange function: In the above, X P2G , X G2P , and X CHP are the P2G, G2P, and CHP variables in the power grid, respectively. Y P2G , Y G2P , and Y CHP are P2G, G2P, and CHP variables in the gas network, respectively. Z P2G is the CHP variable in the heating network.
The iteration format of ADMM-FE is as follows: In the above, C e , C g , and C h are the power grid, gas network, and heat grid costs separated in the objective function, respectively. λ xy,k is the Lagrange multiplier of the grid and gas network coupling variables. λ zy,k is the Lagrange multiplier of the coupled variable of the heating network and the gas network. λ xz,k is the Lagrange multiplier of the coupled variable of the grid and the heating network.
In the above, ρ, τ, and β are penalty parameters. k 1 , k 2 , k 3 are the conversion coefficients of the coupling variables of the CHP equipment between the grid, the gas network, and the heating network. λ P2G,k and λ G2P,k are the Lagrange multipliers of the grid and gas network P2G coupling equipment and G2P coupling equipment respectively. λ CHPxy,k , λ CHPzy,k , and λ CHPxz,k are Lagrange multipliers of CHP coupling devices.
(3) The inner ADMM-FE convergence decision is as follows: In the above, ε x , ε y , and ε z are the absolute values of the maximum difference between the boundary variables of each subsystem after the kth iteration. ε k is the maximum of ε k x , ε k y , and ε k z . When the difference between the boundary variables calculated by the kth and the k+1th time is less than a certain value, and the rate of change of the objective function values between the continuous two-time calculation is less than a certain value, the algorithm is considered to be convergent. The convergence is determined as follows: In the above, a and b are constants used to determine whether or not to converge.

Case Studies
Based on the proposed EGHS scheduling model and solving algorithm, an example is established for analysis. The Example System 1 is as shown in Appendix A Figure A1. It consists of a modified IEEE39 node network, a Belgian 20-node natural gas network, and an 8-node heating network. The Example System 2 is as shown in Appendix A Figure A2. In this process, the heat network of Example System 1 is changed to a 20-node heat network. All tests are performed on a computer with four processors running at 2.30 GHz with 64 GB of memory. Programs are coded using MATLAB R2018a and GAMS.

Example Description
The power network is a modified IEEE39 node power system. There are 46 transmission lines in the system. The original system contains 10 coal-fired units. Set G 1 , G 7 , and G 8 as gas units, which are respectively connected to the 4, 10, and 12 nodes of the natural gas network. Set G 4 and G 5 as the wind turbine connected to the 33 and 34 nodes. In order to maximize wind power, the P2G input is connected to nodes 33 and 34. The system load is set to 85% of the original load, and the cost of wind abandonment is $50/(MW.h). The remaining generator set parameters are recorded in Appendix A Tables A1 and A2.
The natural gas network is a modified Belgian 20-node natural gas network. There are 24 gas pipelines in the system. The system consists of two gas sources, three pressurization stations, and four gas storages. P2G output is connected at 13 and 14. The system load is 20% of the original load, and the gas source energy supply cost is $ 0.3/m 3 . The remaining parameters are recorded in Appendix A  Table A3.
In the example of System 1, the heating network is an 8-node heating network. The system contains 14 water pipes and 2 heat sources. Nodes (1~4) and nodes (5~8) respectively form independent systems, and the two systems are not directly connected. CHP is connected to 1 node and 5 nodes respectively. Heating network parameters are recorded in Appendix A Figure A3.
In the example of System 2, the heating network is a 20-node heating network: the system has a total of 19 water pipelines; CHP is connected at 1 and 8 nodes respectively.

Comparison of PCCP-ADMM-FE and Centralized Algorithm Results
In order to compare the optimization results and convergence speed of the proposed algorithm with the centralized algorithm, the optimization results of the example system under the ADMM-FE algorithm are compared with the centralized optimization results. Set ρ 0 = 10 −3 , the optimization results are shown in Table 1. The analysis of Table 1 is as follows: (1) Under the constraint of convergence accuracy a = 0.001, b = 0.0001, ε 1 = 0.001, ε 2 = 0.001, both the proposed algorithm and the centralized algorithm can get a convergent optimal solution. Compared with the results of centralized optimization, the error of the calculation result of the proposed algorithm is less than 10 −3 . Optimization results verify the accuracy and effectiveness of the algorithm; (2) The ADMM-FE calculation time is longer than the centralized calculation time, but the information privacy of each subsystem is guaranteed. From the scale of formulating the daily dispatch plan, it is acceptable to sacrifice a small amount of calculation speed in exchange for the decentralized autonomy of each subsystem.

Convergence Comparison between ADMM-FE and Directly Generalized ADMM
Taking the iterative process of obtaining the initial point as an example, the convergence of the ADMM-FE algorithm and the directly generalized ADMM algorithm is compared and analyzed. According to the algorithm parameter settings, under the constraints of a = 0.001 and b = 0.0001, the maximum imbalance amount ε k in the optimization process is shown in Figure 2. The ratio of the target difference between the two iterations is shown in Figure 3.
As can be seen from Figures 2 and 3, in Example System 1, both ADMM and ADMM-FE can obtain convergence solutions, but in Example System 2, when the 66th iteration of ADMM-FE is completed, the boundary variable imbalance has reached the convergence condition, but the boundary variables cannot converge in ADMM. It can be seen that in solving the 3-Block problem, the directly popularized ADMM does not necessarily converge, but the ADMM-FE algorithm proposed in this paper has good convergence characteristics and convergence speed.

Convergence Comparison between ADMM-FE and Directly Generalized ADMM
Taking the iterative process of obtaining the initial point as an example, the convergence of the ADMM-FE algorithm and the directly generalized ADMM algorithm is compared and analyzed. According to the algorithm parameter settings, under the constraints of a = 0.001 and b = 0.0001, the maximum imbalance amount ε k in the optimization process is shown in Figure 2. The ratio of the target difference between the two iterations is shown in Figure 3.   As can be seen from Figures 2 and 3, in Example System 1, both ADMM and ADMM-FE can obtain convergence solutions, but in Example System 2, when the 66th iteration of ADMM-FE is completed, the boundary variable imbalance has reached the convergence condition, but the boundary variables cannot converge in ADMM. It can be seen that in solving the 3-Block problem, the directly popularized ADMM does not necessarily converge, but the ADMM-FE algorithm proposed in this paper has good convergence characteristics and convergence speed.

Model Optimization Run Results
We take Example System 1 as an example and analyzed the model optimization results. PCCP-ADMM-FE was used to optimize Example System 1. The calculation example sets the scheduling period T = 24 h, and the simulation step size was 1 h. Figures 4-6 respectively show the power network conditions, gas network conditions, and heating network conditions in the optimization results.
Appl. Sci. 2020, 10, 1214 13 of 18   It can be seen from Figure 4 that in the power network, in the early morning hours, the power load is low, and the wind turbine output is high. At this time, the power of the P2G device can convert the wind power into natural gas storage to absorb the wind power. After 7 o'clock, the daytime work starts, the electric load increases, the grid can absorb all the wind power, so the power of the P2G device drops to 0. In the evening, the power load reaches a peak, and the output of the coal-fired unit    It can be seen from Figure 4 that in the power network, in the early morning hours, the power load is low, and the wind turbine output is high. At this time, the power of the P2G device can convert the wind power into natural gas storage to absorb the wind power. After 7 o'clock, the daytime work starts, the electric load increases, the grid can absorb all the wind power, so the power of the P2G device drops to 0. In the evening, the power load reaches a peak, and the output of the coal-fired unit reaches the upper limit. At this time, the gas turbine output with higher output cost is called to    It can be seen from Figure 4 that in the power network, in the early morning hours, the power load is low, and the wind turbine output is high. At this time, the power of the P2G device can convert the wind power into natural gas storage to absorb the wind power. After 7 o'clock, the daytime work starts, the electric load increases, the grid can absorb all the wind power, so the power of the P2G device drops to 0. In the evening, the power load reaches a peak, and the output of the coal-fired unit reaches the upper limit. At this time, the gas turbine output with higher output cost is called to It can be seen from Figure 4 that in the power network, in the early morning hours, the power load is low, and the wind turbine output is high. At this time, the power of the P2G device can convert the wind power into natural gas storage to absorb the wind power. After 7 o'clock, the daytime work starts, the electric load increases, the grid can absorb all the wind power, so the power of the P2G device drops to 0. In the evening, the power load reaches a peak, and the output of the coal-fired unit reaches the upper limit. At this time, the gas turbine output with higher output cost is called to maintain the power balance. At other times, the coal-fired unit's own output can meet the load demand, so the gas turbine output has been maintained at the lowest level. Therefore, the gas turbine output will maintain the lowest value during the period except for fluctuations from 17 to 19 h. The electrical output of the CHP device is determined by the thermal output. Therefore, in addition to the CHP heating for the demand of heating networks during 6~12 h and 16~22 h, the CHP output grid power is zero.
It can be seen from Figure 5 that in the gas network, in the early morning hours, the gas load and the gas consumption of the gas turbine are low, and the grid inputs high natural gas for the wind power converted by P2G, so the gas tank input is high with 0 output, and the surplus natural gas is stored in the gas storage tank used as a slow storage energy. During the daytime working hours, as the air load increases and the heating-network-demand CHP increases, the gas tank input decreases and the output increases. At 9 h, the gas load peaks and the gas tank output peaks.
It can be seen from Figure 6 that in the heating network, in the early morning hours, the external heat source energy supply cost is lower than the CHP and is at a low price. Therefore, the heat load is satisfied by the external heat source, and the energy supply is greater than the heat load, so the excess heat is used as the heat storage energy. During the daytime working hours, the external heat source energy supply cost is higher than the CHP, so the heat load is transferred to be satisfied by the CHP energy supply. The CHP energy supply is greater than the heat load in some periods, so the excess CHP energy supply is used as the heat network energy storage to eliminate the storage from the gas network during nighttime. From 16 o'clock to 21 o'clock, the heat storage and external heat sources meet the heat load requirements.

The Effect of Heating Network Dynamic Properties on Optimal Operation
In order to explore the influence of the dynamic characteristics of the heating network on the overall optimization operation of the system, this paper draws the correlation between CHP output, external heat source energy supply and heat load under the condition of neglecting the dynamic characteristics of the heating network and considering the dynamic characteristics of the heating network. Ignore the dynamic characteristics of the heating network, that is, ignore Equations (11), (12), (13), and replace them by the following Equation (45) instead: As shown in Figure 7, under the condition of neglecting the dynamic characteristics of the heating network, the CHP output and the external heat source output are constrained by the heat load, and the output and load are coincident. After considering the dynamic characteristics of the heating network, the CHP output and the heat source output are no longer subject to the heat load. During 5~6 h, 16~22 h, the heat load is higher than the heat output, and the heat load is satisfied by CHP, external heat source and heat storage. In the rest of the time, the heat output is higher than the heat load, so the excess heat energy is stored in the heat network as the heat storage energy.
It can be seen that after considering the dynamic characteristics of the heating network, it can be used as energy storage to decouple the heat source output and heat load to a certain extent. When the heating cost is low or the renewable energy needs to be consumed, the heating network performs energy storage, and the energy storage is called during the peak of heating, which improves the flexibility of system scheduling. load, and the output and load are coincident. After considering the dynamic characteristics of the heating network, the CHP output and the heat source output are no longer subject to the heat load. During 5~6 h, 16~22 h, the heat load is higher than the heat output, and the heat load is satisfied by CHP, external heat source and heat storage. In the rest of the time, the heat output is higher than the heat load, so the excess heat energy is stored in the heat network as the heat storage energy. It can be seen that after considering the dynamic characteristics of the heating network, it can be used as energy storage to decouple the heat source output and heat load to a certain extent. When the heating cost is low or the renewable energy needs to be consumed, the heating network performs energy storage, and the energy storage is called during the peak of heating, which improves the flexibility of system scheduling.

Conclusions
In this paper, a day-ahead scheduling model of EGHS considering the dynamic characteristics of the networks is established, and a convex distributed optimization algorithm for inner and outer layer co-ordination based on the ADMM-FE algorithm is proposed. Simulation analysis verifies that the algorithm in this paper can effectively obtain a convergent optimal solution in the three-region distributed optimization problem of the EGHS under the condition of protecting the privacy of the information of each subsystem.
In reality, each subsystem may have different optimization objectives besides cost requirements, and there may be conflicts between different optimization targets. Therefore, in the follow-up work, the game between the power system, natural gas system, and heating system in an open market environment and the multi-objective dispersion optimization problem of EGHS will be the main research direction. In addition, the calculation speed of the ADMM-FE algorithm is slower than the centralized algorithm. In the next work, the author will further study how to improve the calculation speed of the ADMM-FE algorithm.

Acknowledgments:
The authors gratefully acknowledge the support of the National Natural Science Foundation of China.

Conflicts of Interest:
The authors declare no conflict of interest.  In reality, each subsystem may have different optimization objectives besides cost requirements, and there may be conflicts between different optimization targets. Therefore, in the follow-up work, the game between the power system, natural gas system, and heating system in an open market environment and the multi-objective dispersion optimization problem of EGHS will be the main research direction. In addition, the calculation speed of the ADMM-FE algorithm is slower than the centralized algorithm. In the next work, the author will further study how to improve the calculation speed of the ADMM-FE algorithm.

Appendix B
In Equations (2) and (3), M1 and M2 are constants depending on the parameters of the pipelines. The specific formulations are as follows (SI Units): In the above, where f and D denote the friction factor and diameter of the gas pipeline. T means the average temperature of natural gas in the pipeline. Ps and Ts are the pressure and temperature under standard conditions (101325 Pa, 298 K). Rg is the gas constant for natural gas which equals the standard gas constant divided by the molar mass of natural gas.