Interactive Equilibrium of Electricity-Gas Energy Distribution System and Integrated Load Aggregators Considering Energy Pricings: A Master-Slave Approach

Integration of electricity and gas distribution networks improves energy utilization and alleviates environmental pollution. In order to obtain a reasonable energy prices and optimize the day-ahead scheduling scheme under the premise of considering the interests of both supply and demand, an interactive equilibrium model considering electricity-gas energy distribution system and integrated load aggregators is established in this paper. In the upper level, the integrated energy distribution network (IEDN) maximizes the profit via solving the day-ahead optimal energy flow considering energy prices changing in the next day. In the lower level, integrated load aggregators (ILA) minimize charging costs, thus optimizing energy distribution based on integrated demand responses. By using the Karush-Kuhn-Tucker (KKT) optimality conditions and the duality theorem of linear programming, the master-slave game model is transformed into a single layer optimization problem. To solve this problem, second-order cone programming (SOCP) and penalty convex-concave procedure (PCCP) are used to deal with the non-convex functions of electricity-gas distribution system. Therefore, the proposed model is approximated by a mixed-integer second-order cone programming problem. Finally, case studies verify the effectiveness of the proposed model and method.

gas price of the mth integrated load aggregator P L,t net loss of power distribution network η L net loss cost coefficient in power system P T,t active power injected into power distribution gateway from the transmission network f j,t gas output of gas distribution station γ j cost coefficient of the jth gas distribution station P G,i,t active power output of the ith generator VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Q G,i,t reactive power output of the ith generator P D,i,t active power load of the ith node Q D,i,t reactive power load of the ith node P ij,t active transmitted power flow through the branch i − j Q ij,t reactive transmitted power flow through the branch i − j V i,t voltage amplitude of the ith node in power system θ ij,t voltage phase between the ith and jth node in power system L ij,t line pack of the pipe i − j f in ij,t gas flow input of the pipe i − j f out ij,t gas flow output of the pipe i − j f ij,t average of the pipe i − j gas flow k ij gas transmission coefficient of the pipe i − j b i,t gas pressure of the ith nodē b ij,t average of the pipe i − j gas pressure f g,t gas distribution station output vector f D,t gas load vector B g correlation matrix of node-gas distribution station in natural gas system B D correlation matrix of node-gas load B pipe correlation matrix of pipe-node f in average of the pipe i − j gas pressure f g,t gas distribution station output vector f D In recent years, with the increasing problem of environmental pollution, the resource allocation capacity of traditional energy networks has been unable to meet the needs of future development [1]- [3]. In 2019, new natural gas reserves reached 785.645 billion cubic meters in China, but natural gas still accounted for a significantly small proportion of primary energy consumption. With the development of gasfired power generation, combined heat and power (CHP), and power-to-gas (P2G) technologies, the importance of natural gas in energy systems has become increasingly significant. Because natural gas has the advantages of high calorific value, low carbon emissions, and easy storage in large quantities, the electric-gas integrated system that couples power network and natural gas network has gradually become a basic form in energy Internet, and the energy carrier also makes the transition from power flow to gas-power energy flow. The natural gas network is widely deployed in large and medium-sized cities in China, therefore, natural gas network can be coupled with the power network at key nodes. Integrating multiple networks can improve energy efficiency and enhance the reliability of power supply [4]- [6].
There are two main ways of a coupling power system and natural gas system, one is to establish a coupling on the network side through gas turbines and P2G conversion equipment, and the other one is to form an interconversion between two energy sources through an energy hub on the load side. The second way is used to model the supplydemand interaction problem of the integrated energy system. At present, the interaction between supply side and demand side of integrated energy systems has been studied in academia [7]- [9]. References [10]- [12] explored interactive optimization methods for users and operators of integrated energy systems in smart buildings, residential communities, industrial parks, etc. Besides, literatures [13]- [16] regarded the energy hub as an integrated load aggregators, and they proposed a two-level optimization model of the interaction between the integrated energy demand side and the power generators and gas sellers.
As an extension of electricity demand response, integrated demand response provides an entry point for the interaction of both sides of integrated energy supply and demand. Work done in [17] established a demand response load model from the perspective of flexible loads involved in grid operation regulation. By analyzing the relationship between the electricity price difference and the load transfer rate in each period, reference [18] established the time-sharing price optimization model of electricity price and user response, and provided the correct method for the model parameters.
As an effective method to promote the intelligent decisionmaking of multiple agents, game theory has excellent significance for the establishment of supply-demand interaction models [19]- [24]. Work done in [25] developed a masterslave game model based on the grid electricity price and the aggregator pricing mechanism. References [13], [26]- [28] used energy hubs to link energy consumers and distribution networks.
However, the current research on electricity-gas energy system has two problems, as follow: 1) The current research is to treat regional-level integrated energy systems as a single research object to solve problems in planning, construction, and operation. The higherlevel energy supply network is usually simplified and only maintains its supply-demand balance relationship to loads.
2) The flexible scheduling of integrated demand response is largely affected by prices. Most of researches are focused on solving the problem of resource allocation under single or fixed energy prices. But coupling relationship of multi energy conversions in the integrated energy system is complex, and optimal regulation cannot be achieved by relying only on time-sharing electricity price. Therefore, we focus on studying the integrated energy system interactive equilibrium model. The main contributions of this paper are as follows.
1) An interactive equilibrium model between an energy distribution system [29], [30] and energy hubs considering system-level multi-energy flow transmission n a day-ahead scheduling framework is proposed • .
2) The optimization of dynamic energy prices are considered in this paper. The electricity and natural gas prices have complementary characteristics, so making full use of timesharing energy prices can achieve indirect regulation of the comprehensive energy load.
3) A penalty convex-concave procedure (PCCP) are adopted to transform the non-convex dynamic gas flow into a second-order constraint programming 4) KKT conditions and strong duality theory are used to transform the master-slave game model into a mixed integer second-order cone programming (MISOCP). The IEEE 33node distribution network and 20-node gas distribution network are employed to demonstrate the validity and feasibility of the proposed model.
The remaining of this paper is organized as follows. Section II establishes the master-slave game modeling of IES. Section III introduces an improved solution method. Case studies are carried out in Section IV. Finally, Section V briefly summarizes the paper.

II. MASTER-SLAVE GAME MODELING OF IES A. PHYSICAL AND MARKET STRUCTURE OF IES INTERACTIVE EQUILIBRIUM MODEL
The integrated energy supply and demand interaction framework constructed in this paper is shown in Fig. 1. The physical structure can be divided into four parts, including multienergy transmission and supply, multi-energy flow coupling, multi-energy storages and multi-type demands. Multi-energy transmission and supply are provided by the power distribution system, gas distribution system, and wind turbines in small and medium-sized areas. Multi-energy flow coupling refers to the transformation of energy forms through gasfired power generation, CHP, and P2G. Multi-energy storages can include equipment for power storage, heat storage, and gas storage. Multi-type demands refers to the collection of multiple energy requirements of terminal users, which can include power loads, gas loads, and heat loads.
The market structure of IES interactive equilibrium model has three main participants: energy supplier, integrated load aggregators (ILA), and terminal loads. The integrated energy distribution network (IEDN), which consists of energy supplier distribution system and gas distribution system, is the maker of electricity and natural gas prices. ILA is built in this paper to combine many small and medium-sized electric-gasheat loads and take centralized control over them. Because energy hub (EH) can represent energy injection, conversion, transmission, and consumption in of the regional integrated energy system, it is a vital characterization form of an ILA. It can use the collaborative optimization of electricity-gas energy on the supply side and the demand side.
The trading scenarios of IEDN and multiple ILAs characterized by energy hubs are studied in this paper. Energy supplier (IEDN) and ILAs are independent stakeholders. During day-ahead scheduling periods, IEDN maximizes its total profit by optimizing operating conditions and energy pricing strategies for ILAs. Moreover, ILAs determine the optimal energy requirements based on the energy prices set by IEDN. The fluctuations in energy prices will inevitably lead to a shift in the timing of ILA's energy demands and the replacement of energy types, so IEDN and ILAs constitute a master-slave game relationship. Hence, in the upper layer of the masterslave game, IEDN is a leader to develop the energy prices, and in the lower layer, ILAs are followers respond to the change of IES energy prices. Therefore, an interactive equilibrium model of IEDN and multiple ILAs is constructed for dayahead integrated energy market clearing.

B. DAY-AHEAD CLEARING OF ENERGY SUPPLIER 1) OBJECTIVE FUNCTION
The day-ahead clearing model of IEDN is solved by centralized optimization in this paper. The objective of IEDN is to maximize the total profit, which is the difference between sales revenue and operating costs. Besides, the operating costs are the sum of purchased costs from the transmission network and the loss cost of power distribution network. The revenue of rigid load sales from other non-ILA access nodes is calculated as a constant.
where F 1 is the total profit of IEDN; T is the set of total scheduling periods; Hub is the set of ILAs; P m,t and f m,t indicate the electricity demand and the natural gas demand of the mth ILA, respectively; C p m,t and C g m,t define the electricity price and natural gas price of the mth ILA; E non is the sales revenue from rigid loads of non-ILA access nodes; P L,t is the net loss of power distribution network; η L is the network loss cost coefficient; P T ,t is the active power injected into the power distribution gateway from the transmission network; a and b are the cost coefficients; g is the set of gas distribution station; f j,t and γ j are respectively the gas distribution station output and the cost coefficient of the jth gas distribution station.

2) CONSTRAINTS OF IEDN
Constraints of IEDN include limitations in power distribution network, natural gas distribution network, and energy pricing constraints of IEDN.
Constraints of the power distribution network include AC power flow equation and inequality constraints, which can be summarized as follows: where N (i) is the set of nodes connected to the ith node; P G,i,t and Q G,i,t are respectively the active injection and reactive injection of the ith node; P D,i,t and Q D,i,t are respectively the active loads and reactive loads of the ith node; P ij,t and Q ij,t are respectively the active and reactive transmission power; g ij , b ij and θ ij,t are respectively the conductance, susceptance and voltage phase between the ith node and jth node; V i,t is the voltage of the ith node; node is the set of nodes in power distribution network.

b: NODAL VOLTAGE
The node voltage amplitude and phase angle should be limited within the upper and lower bound, as: where V max  The power between the distribution network and the transmission network should be within a specific fluctuation range, which can be described as: where P max T and P min T are respectively the maximum and minimum active power injection at distribution network gateway; Q max T and Q min T are respectively the maximum and minimum reactive power injection at the distribution network gateway.
Compared with the gas transmission network, the gas distribution network has no components such as compressors VOLUME 8, 2020 and valves, and the pipe diameter is relatively small. The high latency effect of the natural gas network allows some natural gas to be stored in the pipeline for a short time, so it is necessary to consider the impact of the line pack.

d: PIPES
The gas flow of a pipe is closely related to the transmission characteristics of the pipe, and the pipes are a non-cyclic network structure. According to the practical gas system operation, the intraday gas flow direction rarely change [31], [32], hence the natural gas flow direction can be determined, wheref ij,t is the average gas flow of the i − j pipe; k ij is the pipeline transmission coefficient; b i,t is the gas pressure of the ith node; f in ij,t and f out ij,t are gas inflow and outflow, respectively; GB is the set of gas nodes.

e: LINE PACK
According to the dynamic characteristics of natural gas network, the input gas inflow is not equal to the gas outflow of a pipe. Line pack equivalent to the reserve in the power system describes the gas stored in a pipe, which can be expressed as follows: where L ij,t indicates the line pack; M ij is a constant related to the physical properties of a pipe;b ij,t indicates the average gas pressure.
where f g,t and f D,t are respectively gas distribution station output vector and natural gas demand vector; B g , B D and B pipe are the correlation matrixes of node-gas station, nodegas load and pipe-node, respectively; f in l,t and f out l,t are the vector s of gas inflow and outflow.

g: VARIABLES CONSTRAINTS
The constraints of the gas pressure and gas station outputs can be described as follow: where b max I and b min I are respectively the maximum and minimum gas pressure of the ith node; f max j and f min j are the maximum and minimum gas station output, respectively.
The power pricing and natural gas pricing constraints of the IEDS include the upper and lower bounds and daily average price constraints, which can be expressed as follow: (13) where λ p t and λ g t are the electricity price and natural gas price of day-ahead integrated markets at time t, respectively; C p, m,t and C g m,t indicate the electricity price and natural gas price of the mth ILA, respectively; α 1 , α 2 , β 1 and β 2 are respectively the electricity price and the natural gas price adjustment coefficient; |T | is the number of elements in the set of total dispatch periods.
Based on the electricity prices and natural gas prices of day-ahead integrated markets, the energy pricing strategy of IEDN and ILAs is established. The fixed daily average prices can avoid IEDN raising energy prices to the upper limits in order to increase its total profit. Moreover, the energy prices rise in a certain period will lead to the prices drop in other periods, so ILAs can optimize its energy requirements by responding to prices.

C. ILA OPTIMIZATION MODEL
In the lower layer of the game, the energy hubs used as the form of ILAs are the follower reacts to the change of system energy prices. The internal energy conversions of ILAs include combined heat and power (CHP) units, transformers, gas boilers, power to gas (P2G) and various energy storage devices. ILAs realize the transfer of natural gas and electricity and achieve multi-energy shifting through the flexible configuration of energy conversion devices. Hence, the ILA m optimization model can be constructed as follows: L e m,t = P m,t η e Tran − P P2G,m,t + f CHP,m,t H GV η e CHP +P out S,m,t − P in S,m,t As (14) shows, the objective of an ILA is to minimize energy purchase cost which is expressed by F hub,m . Equations (15)- (17) imply that the imported electricity, natural gas and heat from the IEDN and conversion devices supply the power load, gas load, and heat load of each ILA. And the output power of conversion devices is equal to the product of input power and conversion efficiency. In (15)- (17), L e m,t , L h m,t and L g m,t are respectively the power load, gas load, and 70532 VOLUME 8, 2020 heat load of the mth ILA; P P2G,m,t and η g P2G indicate the input electricity and conversion efficiency of power to gas; f CHP,m,t denotes the natural gas input of CHP, η e CHP and η h CHP are respectively the power generation efficiency and heat production efficiency of CHP; f Fired,m,t is the natural gas input of the boiler, η h Fired denotes the heat production efficiency of the boiler; P in S,m,t and P out S,m,t indicate the electricity storage and release of the power storage equipment, respectively; f in S,m,t and f out S,m,t are the gas storage and release of the gas storage equipment. In addition, various energy storage equipment capacity equation can be uniformly expressed by (18), S m,t is energy storage capacity, Z in S,m,t and Z out S,m,t are energy storage and release of equipment, µ in m and µ out m denote respectively the storage and discharge efficiency, S cap m indicates the rated energy storage capacity. As (19) shows, Q m,t limited to the maximum and minimum value is the set of input power of conversion devices in the ILA and denoted by Q m,t = [f CHP,m,t f Fired,m,t P P2G,m,t Z in S,m,t Z out S,m,t S m,t ].

D. NASH EQUILIBRIUM EXISTENCE OF THE MASTER-SLAVE GAME
The existence of the master-slave game equilibrium solution for the supply and demand interaction of IEDN and ILAs can be proved by the definition and theorem as follow: Definition: Suppose N is the leader, X and F 1 (x, y) are its' strategy set and benefit function respectively. Let M is the follower, Y and F 2 (x, y) are its' strategy set and benefit function respectively. Leader N chooses the strategy x ∈ X, and follower M will choose y ∈ Y from its best response K (x) = {w ∈ Y : F 2 (x, w) = min F 2 (x, u)} after knowing the leader's strategy, and the leader N will choose x * ∈ X to make max y∈K (x * )  Theorem: Let F i : X × Y → R(i = 1, 2) be continuous on X × Y , and ∀x ∈ X, F 2 (x, y) is strictly convex, then ∃x * ∈ X and y ∈ K (x * ) makes To prove the existence of the master-slave game equilibrium solution, the following two problems must be explained: 1) The benefit functions of the game participants are continuous functions of the decision variables; 2) For the leader's decision x * ∈ X, the follower benefit function is strictly convex.
In this paper, the leader is IEDN and the benefit function is (1). The follower is ILAs in the lower layer, and the benefit function is expressed by formula (14). Due to the continuity of two decision variables obviously, 1) is proved. The benefit function of the ILAs are linear functions, so 2) is proved. As proved above, the Nash equilibrium of the game exist.

III. SOLUTION METHOD
In the objective of IEDN, C p, mt , C g m,t , P m,t and f m,t are decision variables. Moreover, the AC power flow function (2) and the gas flow function (5) are non-convex equality constraints, so the optimization model is still non-convex. Hence, how to eliminate the non-convex constraints of this mode to ensure the quality of the master-slave game equilibrium solution is also the research focus in this paper.

A. RELAXATION AND CONVEXITY
Due to non-convex equations (2) and (5) in the upper layer optimization model, the quality of master-slave game equilibrium solution is difficult to insure. Therefore, second-order cone programming (SOCP) is used to relax the AC power flow equation (2) of power distribution network, and penalty convex-concave procedure (PCCP) is proposed to convex the gas flow equation (5) for gas distribution network pipelines.

1) POWER DISTRIBUTION NETWORK CONVEXITY
The non-convexity of AC power flow (2) is derived from the trigonometric functions and negative quadratic terms in the equation. Since SOCP is convex optimization, the solution optimality and calculation efficiency can be guaranteed. First, replace the variables in Eq. (2), as: Then the branch active power and reactive power equations can be converted as follow: Only the last equation in Eq. (21) is a nonlinear constraint, so it is necessary to relax this equation as follows: Hence, Eq. (22) is a second-order cone constraint so that the non-convex source in the power distribution network is eliminated.

2) GAS DISTRIBUTION NETWORK CONVEXITY
Two inequality can be employed to replace Eq. (5), as: Through Eq. (23) is a standard second-order cone constraint, Eq. (24) is a non-convex restriction that need to be convexity [34], [35]. Using the first-order Taylor expansion to linearize the subtracted term in Eq. (24): where n represents the number of iterations; b n−1 j,t is the pressure value of jth node after the n − 1th iteration;f n−1 ij,t VOLUME 8, 2020 is the average natural gas flow of pipe i − j after the (n − 1)th iteration; δ ij ,t is the relaxation variable of the pipe i − j. The penalty term of relaxation variables is added to the objective of IEDN, and the PCCP is used to transform the problem into an iterative process which tights relaxation domains gradually.
B. MASTER-SLAVE GAME EQUILIBRIUM SOLUTION 1) SOCP EQUIVALENT CONDITION Due to energy prices known in the under layer, the ILAs optimization model is a linear model. Therefore, it can be converted into the equivalent KKT (Karush Kuhn Tucker) conditions and used as the constraints of the upper optimization model. Thus, the master-slave game model is transformed into a single-layer optimization model. Firstly, based on Eq. (14), the Lagrangian function L m can be defined as follows: (26) where h is the set of equality constraints; h j is the j th equality constraint; µ j is the dual variable of the jth equality constraint; d is the set of inequality constraints; g k is the kth inequality constraint; λ k is the dual variable of the kth inequality constraint.
Then, the ILAs optimization model Eq. (14)- (19) can be transformed into the KKT conditions as follows: The second term in Eq. (27) describes the complementary relaxation condition, which is nonlinear. So the Boolean variable θ k is introduced to convert (27) into a linear inequality (28), (29): where θ k is a binary number, M is a sufficiently large integer.

2) OBJECTIVE LINEARIZATION
The transformed single-layer optimization model is still nonconvex, and the non-convex source is the IEDN objective Eq. (1). According to the strong duality theory of linear programming, the original ILAs objective function corresponding to the optimal solution is equivalent to the objective function value of the dual problem. Therefore, the two variables in Eq. (1) can be linearized by the following form: where µ e , µ h and µ g are the dual variables corresponding to the equality constraint (15), (16) and (17), respectively; P max P2G,m,t , P min P2G,m,t are the input power limit of P2G; f max CHP,m,t , f min CHP,m,t are the CHP unit intake natural gas limit; f max Fired,m,t , f min Fired,m,t are the gas boiler intake natural gas limit; P in,max S,m,t , P in,min S,m,t , P out,max S,m,t , P out,min S,m,t indicate the storage capacity and discharge limit of the storage equipment; E max m,t and E min m,t are respectively the capacity limit of the storage device; f in,max S,m,t , f in,min S,m,t , f out,max S,m,t , f out,min S,m,t are the gas storage capacity and discharge limit of the gas storage equipment; G max m,t , G min m,t are the capacity limit of the gas storage equipment; λ is the dual variables that limit the constraints.

C. PCCP ALGORITHM
The process of PCCP for the above single layer optimization model is presented as follows.
In conclusion, the master-slave game model solving flow chart is shown in Fig. 2.

IV. CASE STUDIES
In this section, case studies are carried on a testing integrated energy system combined with an IEEE 33-bus distribution system and a modified 20-bus natural gas distribution system to establish a master-slave game model and investigate the application performance of the proposed convexity methods. The case studies are executed with Matlab R2017a and GAMS.

A. SIMULATION MODEL
As shown in Fig. 3, the IEEE 33-bus distribution network is connected to the transmission network at node 1. The cost coefficients of the upper power grid and network loss cost coefficient of the distribution network are a = 0.827$/(MWh) 2 ., b = 23.341$/(MWh) 2 , η L = 25$/(MWh) 2 . Also, wind farms are linked to nodes 13 and 29. For the IEEE 33-bus distribution system, the maximum demand is reduced to 75% of the original. Based on the network topology of the Belgian 20-bus gas system, a 20-node natural gas distribution network was constructed. This network has two gas distribution stations and 24 natural gas pipelines. The supply cost of gas distribution stations is 0.25$/m 3 . Compared with the original Belgian system, the modified gas distribution network has no compressors, valves and other components, and the pipe diameter is relatively small.
In Fig. 3, IEDN is connected to eight ILAs (represented by EH 1 to EH 8 ) in the simulation model, and the internal structures are different shown in Fig. 4. The day-ahead market energy prices are shown in Fig. 5 and Fig. 6 [36]- [39], and typical daily load data are presented in Table 1. Moreover, the loads of ILAs are proportional to the data in Table 1.

VOLUME 8, 2020
Algorithm 1 PCCP Algorithm 1) Ignore the constraint (25), solve a relaxed SCOP to get an initial point (x 0 , W 0 ): Initialize n = 1, initialize the penalty factor ρ 1 , and set the convergence tolerance ε 1 and ε 2 of PCCP. Set the dynamic adjustment ratio v c of the penalty factor, and v c > 1.
2) Take (25) into account, and solve the convex SOCP problem to get the nth point (x n , W n ): where ρ n is the penalty factor of nth iteration.
3) Convergence criteria check: If the above convergence conditions are satisfied, the calculation is stopped and the optimal solution is found. otherwise, the penalty factor is updated according to the following formula: where v c is the dynamic adjustment coefficient. Then set n = n + 1, and repeat steps 2, 3 until the convergence criteria are met.
Also, the internal equipment parameters of ILAs are shown are presented in Table 2.

B. PCCP PARAMETERS SETTING ANALYSIS
The initial value of the penalty factor ρ 0 is set to 1 × 10 −3 , and the dynamic adjustment coefficient v c is set to 2.5. As shown as Fig. 7, the relaxation variables of the natural gas pipeline after convergence are generally concentrated at 1.1 × 10 −4 to 1.25 × 10 −4 , so it can be approximated that Eq. (24) and Eq. (25) are equal. It can be seen that the PCCP algorithm is effective for solving this model.

C. ENERGY PRICING SENSITIVITY ANALYSIS
To research the impact of energy price optimization on scheduling operation of IEDN, the following four cases are studied and compared. Case 1: Mandatory scheduling case for fixed energy prices. Case 2: Natural gas prices are fixed, only to optimize electricity prices in one day.
Case 3: Electricity prices are fixed, only to optimize natural gas prices in one day.
Case 4: Optimize electricity prices and natural gas prices at the same time in one day.
According to the simulation results in Table 3, for the case of optimizing single energy price, the total profit of IEDN is improved compared with the mandatory scheduling case. Moreover, the case of optimizing two energy prices at the same time is better than Case 2 and Case 3, which can reduce ILAs total purchase energy cost and the system operating cost. This is because ILAs achieve a rational optimization   of energy demand through timing shift of multiple loads and different energy sources replacement under energy prices guidance. Hence the system energy prices obtained based on the interactive equilibrium model proposed in this paper have better flexibility than the fixed system energy prices.
Based on Case 4, the influence of energy pricing constraints on the equilibrium solution is explored in this section. Change the lower limits of electricity price and natural gas price α 1 , β 1 from 60% to 90%, and remain other parameters unchanged. Table. IV shows the impact of different pricing limits on total profits of IEDN and the energy purchase costs of ILAs.
As the lower limit of energy pricing increases, the total profit of IEDN and the energy purchase cost of ILAs show a downward trend. When the lower price limit is lower, due to the daily average prices constraints, the electricity price is lower during the period from 02:00 to 07:00, which is the periods of low load. And it is higher between 12:00 and 16:00, which is the periods of high load. However, the change in gas prices is the opposite. In this case, the energy storage can't increase the energy input continuously because of being limited by the maximum capacity, and CHP units are affected by the gas prices to maintain the lower limit of output. Therefore, the change in electricity prices has a greater impact on the total profit of IEDN and the energy purchase cost of ILAs, making them both higher. When the lower price limit is higher, the electricity price is higher in the low periods of power loads and is lower in the peak load periods, and the change in gas prices is the opposite. The output of CHP units are dispatched more, so the total profit of IEDN decreases and the energy purchase cost of the ILAs decreases. It can be seen that setting appropriate pricing constraints can balance the interests of both IEDN and ILAs. Table 5 shows the energy purchase costs of the eight ILAs. Comparing optimal objectives of four ILAs with consistent load levels, EH 1 only contains one power storage, as shown in Fig. 4(a), it can respond to the change of electricity prices VOLUME 8, 2020  and realize the timing shift of power, and transfer the peak electric load to the trough periods. EH 3 contains gas turbines, as shown in Fig. 4(b), which can replace natural gas to power during peak hours. As shown in Fig. 4(c), EH 5 includes power storage equipment and gas turbines. It can choose to convert high-efficiency power or natural gas with sufficient supply to meet its energy demand. The EH 7 , as shown in Fig. 4(d), has wind power access, and the cooperation of P2G and CHP units can convert power and natural gas, and realize the complementary advantages of various energy sources. It can be seen that the introduction of integrated demand response promotes the optimization of load energy structure and the improvement of energy efficiency.

2) OPTIMIZATION RESULTS AND DISCUSSIONS
In the simulation model, energy conversion devices and load levels included in ILAs are distinct from each other. Therefore, there is a difference in the energy prices and energy demand allocation structures corresponding to each ILAs  in the master-slave game equilibrium. Taking EH 5 as an example, the optimization results of energy prices and energy conversion devices are shown in Fig. 8.
It can be seen from Fig. 8 that the electricity price of EH 5 tends to the upper limit during the load trough, while the gas price is lower at night and higher in the rest of the periods. Besides, as shown in Fig. 9, the output power of the CHP unit, the power storage, and the boiler in EH 5 are closely related to the energy prices. With the effect of energy prices, the CHP unit works at the peak power load, while it burning natural gas to generate electricity can reduce the energy purchase cost of EH 5 . Because burning natural gas for the rest of periods will increase the energy purchase cost, so it will only be maintained at the lower limit output, and the heat output will also be less at the same time. Also, the power storage is equivalent to the price type demand response, which realizes the transfer of the power load in time according to the change of prices during the scheduling periods. The capacity change curve of the power storage shows that the storage device stores power at night when the power load and the electricity price are low. Then, the power load reaches the first peak in a day from 12:00 to 14:00, and the stored power energy is released. The supply and CHP are dispatched to maintain the electric demand balance. Finally, the power load reaches the peak again at 19:00, and the power storage discharges for decreasing purchase cost. EH 7 and EH 8 contain P2G, power storage, gas storage, boilers, and CHP units, and the access nodes connect to wind farms. Taking EH 7 as an example, energy prices and various types of electric power and gas power are shown in Fig. 9 to Fig.9. As shown in Fig. 9, the P2G in EH 7 works from 01:00 to 7:00 and 23:00 to24:00, at which time the electric price is lower, and the wind power output is larger due to P2G conversion efficiency not high, surplus wind power is first stored through power storage, and the rest of the output that can't be absorbed is converted to natural gas through P2G. Besides, it can be seen from Fig. 10 that due to the large proportion of thermal load on a typical winter day, the ratio of gas boiler and CHP to total natural gas load is higher than the actual natural gas load, and CHP mainly bears the role of peak shaving, and its output depends on power load. The gas storage stores the P2G converted natural gas from 04:00 to 07:00, and the power load and the gas load are highest from 12:00 to 18:00. At this time, the gas storage deflation can maintain the gas load balance and reduce energy purchase costs.
Without changing the charge and discharge limits of power storage equipment, the maximum capacity E max of EH 7 is altered from 30 kWh to 300 kWh. The energy purchase cost changes are shown in Fig. 11. When E max is less than 120 kWh, the energy purchase cost of EH 7 is significantly reduced with the increase of the maximum storage capacity. When E max is bigger than 120 kWh, the growth of the maximum storage capacity does not effect on the energy purchase cost. This is due to the increased capacity of the power storage, so that EH 7 can store more power at a lower price during the lower period of the electric price, and reduce the amount of electricity purchased during the higher electric price period. However, when E max increases to a certain extent, it is limited by the charge and discharge power. Although there is still more storage capacity, it cannot be effectively utilized, so the energy purchase cost doesn't change significantly. Besides, the capacity of the gas storage device is changed in the same way as shown in Fig. 12.

V. CONCLUSION
This paper proposes an interactive equilibrium model of integrated energy distribution system and multiple integrated load aggregators based on a master-slave game. This model guides the decision of the integrated load aggregator by optimizing the energy price to achieve the benefits balance of maximizing the total profit of the energy distribution system and minimizing the energy purchase cost of the comprehensive load aggregator. The case studies with the IEEE 33-node distribution network and 20-node gas distribution network show that: • PCCP algorithm can converge in a short time by continuously punishing the duct gas flow relaxation domain. Therefore, for the day-ahead scheduling problem, the calculation efficiency of the algorithm is completely acceptable.
• The energy price determined by the interactive equilibrium model proposed in this paper is closely related to the energy demand of the integrated load aggregators, therefore the energy pricing mechanism proposed has better flexibility than the energy price of fixed systems.
• The integrated demand response strategies considered in this paper are conducive to improving the multi energy utilization rates. Moreover, increasing the capacity of multiple types of energy storage devices within a certain range can reduce the energy purchase cost of integrated load aggregators. The next stage will consider independent scheduling of distribution networks. That is a master-slave game model including a non-cooperative game between the upper-layer network, and optimization of the energy purchase cost of the lower-level integrated load aggregator. BO YANG received the B.Eng. degree in electrical power system and the M.Eng. degree in electrical power system from the South China University of Technology, Guangzhou, China, in 2010 and 2011, respectively, and the Ph.D. degree in electrical engineering from the University of Liverpool, Liverpool, U.K., in 2015. He is currently an Associate Professor of power system with the Kunming University of Science and Technology, Kunming, China. His main research directions are the optimization and control of new energy power generation systems, and the application of artificial intelligence in smart grids.