Bi-Level Decentralized Optimal Economic Dispatch for Urban Regional Integrated Energy System Under Carbon Emission Constraints

Based on the real-time environmental constraints in urban regional construction, this paper constructed a bi-level decentralized low-carbon optimal dispatching model of the urban regional integrated energy system (RIES), including the park integrated energy systems (IESs). In this model, a bi-level optimal allocation model of carbon emission constraints between the urban and the park is proposed for the first time. The upper urban will formulate the real-time carbon emission constraints based on real-time environmental monitoring, decomposing the historical carbon emissions to the lower park IESs; the lower park will meet the real-time carbon emission constraints during optimization. We through the upper urban with the lower park between the bi-level decentralized optimization to ensure that the objective function’s upper urban power, natural gas, and heat distribution network system is minimum total network loss. In addition, it is necessary to ensure the minimum operating cost of each park IESs and focus on how to meet the requirements of the overall environment of urban RIES. Furthermore, we study the influence of optimal allocation strategy of carbon emission constraints on network loss, and operating cost of urban RIES under different scenarios. Then, an improved analytical target cascading (ATC) method is applied to solve the bi-level decentralized optimal dispatching model of urban RIES. Finally, an example under three different scenarios is given to verify the superiority and effectiveness of the proposed model and the improved method.


I. INTRODUCTION A. BACKGROUND AND MOTIVATION
At present, China is still in the middle and late stage of rapid urbanization development. The energy consumption and carbon emission generated by new urban buildings will continue to increase [1]. However, faced with the continuous expansion of energy demand and the worsening of ecological and environmental problems, the urban Regional Integrated Energy System (RIES) shows a new development trend [2], [3]. Traditional integrated energy system (IES) planning and operation are limited to single energy forms such as electricity, natural gas, cold and thermal energy, which cannot fully play the complementary advantages and synergistic benefits of IES in urban areas. Compared with a single energy system, urban RIES is mainly composed of integrated energy supply and use systems, such as regional electricity, natural gas, and heating distribution systems through IESs interconnection at the terminal, which can realize the transmission, conversion, distribution and other key roles of multi-types of energy [4]- [6]. It can also effectively realize the decentralized coordination and optimization among subsystems in the urban area quickly in urban areas rapidly and fully exploit the potential advantages of different subsystems. Meanwhile, overall consideration of multiple types of energy sources can also significantly reduce the integrated energy consumption of urban RIES [7], [8]. Therefore, it is of great significance to study urban RIES' integrated energy-saving and carbon emissions reduction.

B. LITERATURE REVIEW AND RESEARCH GAP
It was considering that different research subjects showed a diversified distribution of RIES at all levels. At present, many studies have focused on the influence of models, strategies, demand response (DR) or integrated DR mechanism design, and methods on the optimal operation of RIES. Soheil et al. [9] proposed a general energy hub economic dispatch model. A new self-adaptive learning with a timevarying acceleration coefficient gravitational search algorithm was employed to solve the model. In [10] and [11], an optimal operating model for electricity-thermal-natural gas network systems was presented that considered economy and environmental pollution. Li et al. [12] proposed a linear electricity and district heating networks model to coordinate the short-term operation of electric power and district heating systems. And then, in [13], a combined heat and power dispatch model for electric power systems and district heating systems was presented. In [14] and [15], a novel mixedinteger linear programming (MILP) optimization model was introduced for multi-energy networks, in which electricity, heating, and cooling loads and sources were considered in the IES. The authors of [16] proposed an optimized dispatching model for city-scale IES that realized the flexibilities of energy utilization. Pan et al. [17] proposed a reasonable power to heat and hydrogen model with start-up or shut-down constraints and a novel seasonal hydrogen storage model for an electricity-hydrogen IES. Taylor et al. [18] presented a new mixed-integer quadratic, quadratically constrained, and second-order cone programming model for distribution system reconfiguration. The above studies have contributed significantly to IES modeling in their respective ways.
Some operation factors are also considered in the optimal operation of IES, such as the strategies, DR, or Integrated DR programs combined with different application scenarios. For example, He et al. [19] presented an environmental economic dispatch model for IES. The carbon trading scheme and integrated DR strategies were introduced to reduce pollutant gas emissions. Guo et al. [20] presented a RIES optimization model considering DR, in which different DR strategies were applied to improve the energy utilization efficiency. In [21] and [22], an integrated DR uncertainty or certainty model for the community networks was proposed, and a distributed algorithm was employed to determine their optimal strategies. Tek et al. [23] developed an integrated water-foodenergy nexus optimization model for an irrigated agricultural production IES, in which carbon and water strategies were applied to reduce carbon footprints. He et al. [24] presented a low-carbon economic dispatch model for electricity and natural gas systems, in which power to gas (P2G) equipment and carbon emission capture strategies were applied to solve it. The authors of [25] proposed a bi-level optimization model for an IES network, in which a hybrid algorithm with an adaptive particle swarm optimization algorithm was applied to increase the proportion of renewable energy. Xu et al. [26] presented a novel RIES structure in which an energy management strategy was applied to improve energy efficiency. Pan et al. [27] also offered a two-stage MES planning method and novel energy pricing strategy that considered the integrated DR program.
Moreover, some researchers have studied the operation factors (such as different economic dispatching models and fast solving methods or algorithms.) influencing carbon emissions in the optimal operation of IES. For example, Chen et al. [28] proposed a price-based low-carbon hierarchical model for IES, in which a carbon emission flow model was employed to trace the carbon flow and calculate the carbon emission, and an enhanced particle swarm optimization algorithm was applied to solve the model. The authors of [29] proposed a complex non-convex, non-smooth and non-linear multi-objective dynamic low-carbon dispatch model for the combined heat and power units, in which a novel selfadaptive probabilistic mutation method was introduced to improve the performance of the algorithm. The authors of [30] proposed a low-carbon distributed energy resource optimization model, in which a MILP method was employed to minimize the energy cost and reduce carbon emissions. In [31], a long-term co-optimization planning model for an integrated electricity and district heating system was proposed that considered transmission lines and heat pipelines, in which a parallel Benders decomposition combined with the sequential bound-tightening method was applied to solve it. Tsao et al. [32] presented a sustainable advanced distribution management system considering the DR program and carbon emission strategy. In [33] and [34], a bi-level joint economic dispatch model was constructed that considered carbon emission constraints for new urban with multi-energy parks. Sequential quadratic programming combined with the path tracking interior point method was applied to solve it. Furthermore, other studies have focused on the mechanism design of multi-energy prices. For example, Gu et al. [35] proposed a double iterative optimal low-carbon dispatch model for an industrial park with multi-energy, in which a bi-level primedual path following the interior point method was applied to solve it.
The above studies have focused on the two-layer operation of the urban RIES with park IESs considering optimal carbon emission constraints. First, few studies have constructed a centralized or distributed optimal economic dispatch framework for urban RIES with a park IES. Second, the integrated DR program affected the operation of urban RIES by regulating flexible loads, and corresponding constraints were considered. The above economic dispatch model was designed only to ensure that urban carbon emissions are roughly limited. In recent years, China has put forward a series of strict requirements for urban carbon emissions reduction. Hence, some scholars have considered carbon trading scheme in the objective function of the economic dispatch model, introduced low-carbon equipment operations such as power-to-gas devices and ground source heat pumps, and proposed low-carbon emission reduction technologies such as a given total carbon emissions quota index, to realize low-carbon economic dispatch for urban. In addition, previous research on the IDR program corresponding to the multi-energy pricing incentive strategies affected the operation of park IESs by regulating flexible loads. Although few papers have considered the design of urban carbon emission constrained optimal allocation mechanism with park IESs, when urban areas get a given total carbon emission constraint target value, we cannot be sure that the given target value is the best. Therefore, we need a new carbon emission constraint optimization allocation mechanism to achieve the above goal.

C. CONTRIBUTIONS
Motivated by the abovementioned research gap, this paper proposes a bi-level decentralized optimal low-carbon economic dispatch model for urban RIES considering park IESs. This model is optimized through decentralized coordination between the upper and lower levels of the urban RIES and park IESs, and takes the optimal network loss of the urban and the operation cost of park IESs as the optimization objective, respectively. The upper-level optimization involves the electricity, natural gas, and district heating distribution networks, considering the dispatch of two coal-fired and one diesel generator unit and battery devices at the power distribution network side, and considering the urban historical realtime unit power carbon emissions constraints. Overall, the upper level involves real-time low-carbon constraints aimed at the optimal network loss of urban electricity-gas-heating distribution networks. The lower level contains real-time lowcarbon constraints aimed at the optimal economical operation cost of park IESs. The IDR program model of the temperature-controlled loads is constructed. The contributions of this paper are as follows: 1) A bi-level decentralized optimal low-carbon economic dispatch model for urban RIES with consideration of park IESs is proposed, which is beneficial to reduce the total carbon emissions of urban RIES and the economic operation costs of park IESs. 2) An improved analytical target cascading method (ATC) is employed, which can effectively improve the convergence characteristics of the algorithm. 3) A carbon emission constrained optimal allocation mechanism is introduced, which can help develop a more accurate carbon emissions reduction strategy for urban RIES.

D. PAPER ORGANIZATION
The rest of this paper is organized as follows. Part 2 introduces the bi-level decentralized optimization model of the urban RIES. Part 3 presents the solution algorithm for the proposed model in detail. In part 4, case studies and simulation results are presented. Part 5 presents the summary and conclusions. In part 6, future work is presented. Parts 7,8,9, and 10 represent the appendix, acknowledgment, references, and team respectively.

II. BI-LEVEL DECENTRALIZED OPTIMIZATION FRAMEWORK AND MODEL OF URBAN RIEs WITH PARK IESs
The bi-level decentralized optimal scheduling framework between the urban RIES and park IESs is shown in figure 1. The whole structure is divided into two levels: the urban RIES and the park IESs. The urban is responsible for setting and allocating the total carbon emission target value. New carbon emission constraints are reset and decomposed based on the historical carbon emissions generated by the upperlevel urban RIES and the lower-level park IESs, respectively. Under the constraints of real-time carbon emissions, the lower park IESs will transfer the accumulated carbon emissions generated by their actual operation to the upper urban, respectively, and reset the total carbon emission target value of the urban. Therefore, the total carbon emissions generated by the units with carbon emission sources in the urban will be updated and iterated in real-time during the optimized operation. At the upper level, considering the real-time carbon emission and historical data constraints of the unit power supply, the urban is responsible for the decentralized management of the park, reasonably scheduling the coal-fired and diesel generator, and battery units of the urban, integrated power supply, and the optimized operation of energy conversion equipment of the park, aiming to minimize the loss of urban's electricity-gas-heating distribution networks. At the lower level, considering the real-time unit power of the park takes into account carbon emissions constraints, the park is responsible for reasonably scheduling the optimized operation of the multi-type energy equipment in each IES, aiming to minimize parks' comprehensive economic operation cost. The integrated DR power can be housed in the park IESs as the urban RIES terminal energy hub, along with the integrated energy supply, energy conversion, power, gas, cold and heat, energy storage units, and the flexible load combinations by reasonably scheduling the integrated coordination of park IESs. Through the decentralized optimal equipment operation and carbon emission allocation, output energies can complement each other to integrate the urban RIES and park IESs. Figure 2 shows a bi-level structure of the urban RIES with park IESs.

A. UPPER AND LOWER OPTIMIZATION MODEL
At the upper level, the urban is responsible for the unified allocation of the whole total carbon emission target value within the urban itself and the park area. The carbon emission allocation mechanism of urban RIES can be referred to Fig.1 and the corresponding description section. In addition, the upper urban RIES involves the electricity, gas, and heating distribution networks, and integrated two coal-fired generator units and one diesel generator unit, and battery units (DN-Battery) within the power distribution network. The lower level is responsible for the unified allocation of the lower total carbon emission target value to IES within the park. The lower park IESs are responsible for the unified management of the integrated power supply, energy conversion, and the optimal operation of electric power, gas, heat, and cold energy and energy storage equipment within the park. The power supply equipment in the park IES includes micro-gas turbine (MT), wind turbine (WT), photovoltaic (PV), heat recovery (H-REC), and gas boiler (GB). The Energy conversion equipment includes electrical chiller (EC), electrical heater (EH), heat exchanger (H-EX), H-REC, and absorption chiller (AC). The energy storage devices include batteries (BT), cold (CD), heat (HT), and natural gas (NG) storage tanks. The terminal load side of the park with electricity, gas, cold and heat energy, and transfers electricity, cold and gas to users through rigid and shiftable loads (R-EL and S-EL, R-CL and S-CL, R-GL and S-GL), and heat to terminal users through flexible air heating and hot water supply temperature loads (A-HL and W-HL).

B. UPPER OBJECTIVE FUNCTION
The upper level aims to optimize the operation of the coal, diesel generators and energy storage equipment that can be dispatched by the upper layer reasonably on the premise of meeting the load demand of the urban regional distribution network system, and to minimize the network loss of the power distribution, gas, and heat networks respectively. Therefore, the upper optimization objective functions FDN, FNG, and FHWP are respectively shown as follows:

1) OBJECTIVE FUNCTION OF POWER DISTRIBUTION NETWORK
Considering the power transmission loss of the distribution network, this section selects the minimum loss of the distribution network within the day-ahead optimization scheduling period as the objective function. The details are as follows: Here, F DN is the optimization objective function of the distribution network. I t n,n represents the current of the line nn at time t, A. R n,n means the resistance of branch nn , ohms. N /N -1 indicates the number of nodes and branches of the distribution network, respectively. t ∈ T means hourly periods, running from 1 to T hours.

2) OBJECTIVE FUNCTION OF NATURAL GAS NETWORK
Considering that the pressurizer in the upper urban regional natural gas network is mainly used to compensate for the pressure or energy loss in the process of natural gas transmission, this section selects the minimum energy loss in the day-ahead dispatching cycle as the objective function, the details are as follows: where F NG is the optimization objective function of the natural gas distribution network. Q t com,κ,l means the airflow consumed by the pressurizer and the natural gas flow into the pressurizer at time t, m3/h. − Q t κl indicates the average flow rate of the natural gas pipeline of the gas network at time t, m3/h. f κ,l represents the damping coefficient of the natural gas pipeline in branch KL, which is usually 0.06-0.10. K /L indicates the number of nodes and branches where the pressurizer is installed in the branch of the gas network.

3) OBJECTIVE FUNCTION OF HEAT NETWORK
Considering the energy loss in the hot water transmission process of the upper urban area heat network, this paper selects the minimum energy loss in the day-ahead scheduling cycle as the objective function, the details are as follows: where, F HWP is the optimization objective function of the heat network.
Here, P t n,n /Q t n,n represents the active/reactive power of the line nn , (kW). P t n /Q t n means the active/reactive power injected by node n at time t, (kW). Pg t gm µ is the output power of coal-fired or diesel generator set dispatched at time t, kW, where µ ∈ {gm}; Pl t n,u=dl denotes the active power of each node at time t, (kW). Pg t n,c=es /Pg t n,d=es means the charging/discharging power of the battery device at the distribution network side at time t, (kW). P t GD_IESi represents the active power transmitted by the distribution network from node n to each IES in the park at time t, (kW). X n,n denotes the reactance of branch nn , (ohms). V min /V max means the lower/upper limit value of each node voltage, (kV).S n,n represents the upper limit of the apparent capacity of branch nn , (kVA). e ∈ {1, 2, 3, 4} means the variation range of the distribution network's perceived capacity constraint, as shown in figure. 3 where, W t es,dn is the storage capacity of the battery at the distribution network side at time t, (kWh). Cap t dn,bt is the rated capacity of the battery storage device, (kWh). W t es,dn,min and W t es,dn,max are the allowable minimum and maximum rated stored capacity of the battery device, kWh, respectively.
where, R gm µ ,down /R gm µ ,upper is the constraint coefficient of upper/lower limit climbing output for coal-fired and diesel generators, (kW/min), respectively. When the subscript µ is 1 or 3, it represents coal-fired units, and when µ is 2, it represents diesel units.

d: INEQUALITY CONSTRAINTS ON TRANSMISSION POWER BETWEEN POWER DISTRIBUTION NETWORK AND EACH PARK IES
where,P t GD_IESi,max /P t GD_IESi,min is the transmission power between the power distribution network and each IES in the park at time t, kW, respectively.

e: CARBON EMISSION CONSTRAINTS OF UPPER URBAN RIES ITSELF BASED ON REAL-TIME ENVIRONMENTAL MONITORING
where, ρ e 1 and ρ e 2 are the carbon emission intensity coefficients of coal-fired and diesel generating set per unit time, (t/kWh), respectively. D δ carbon,e 1 and D δ carbon,e 2 are the upper limits of carbon emission constraints that can be relaxed for coal and diesel generator units, (t), respectively. E t,δ carbon,e 1 and E t,δ carbon,e 2 represent real-time optimized carbon emissions of coal and diesel generator units, (t/h), respectively. Formula (8) indicates that the real-time carbon emissions of the coal/diesel units that can be allocated by the power distribution network in the upper urban area after optimization, and the carbon emissions are respectively calculated and transmitted to the upper urban. The superscript δ represents the number of iterations.

2) EQUALITY AND INEQUALITY CONSTRAINTS OF THE GAS NETWORK a: POWER BALANCE EQUATIONS AND INEQUALITY CONSTRAINTS FOR GAS NETWORK NODES
Here, Q t Air,s is the gas flow output from the natural gas source point of the gas network at time t, (m 3 /h). Q t κ,load is the natural gas load flow of the pipeline node κ at time t, m3/h. Q t κ,l,in and Q t κ,l,out represents the gas flow at the inlet and outlet of the pipeline κl at time t, (m 3 /h). P t NG_IESi refers to the transmission flow of natural gas pipeline for the gas network from node κ to each IES at time t, (m 3 /h). Q t Air,s,max /Q t Air,s,min means the upper/lower limits of the output gas flow at the air source point of the gas network at time t, (m 3 /h), respectively.

b: INEQUALITY CONSTRAINTS FOR NATURAL GAS PRESSURIZER DEVICES
Friction resistance exists in the natural gas transmission process of the medium and low-pressure gas network, resulting in pressure loss of the gas network, so a certain number of pressurizer devices are needed. This paper assumes that a variable ratio fixed pressurizer device is adopted [14], The details are as follows: where p t l is the pressure of node l of the natural gas pipeline at time t, (bar). β com,c is the compression coefficient of the pressurizer, which is the pressure boost ratio of the pressurizer, and the value is set to 1.2. p t κ is the pressure of node κ of natural gas pipelines at time t, (bar).

c: EQUATION AND INEQUALITY CONSTRAINTS OF NATURAL GAS PIPELINE FLOW
Considering that there is no general formula to describe the steady flow of natural gas in pipelines, we adopt the simplified Polyflo's medium-pressure pipe network equation to describe the branch gas flow of the natural gas pipe network and the pressure at both ends of the pipe nodes [14], the details are as follows: , and where, − Q t κl is the average flow rate of natural gas pipelines at time t, (m3/h). K κl represents the pipe constant of the pipe κl at time t, which is related to the length L (m) and diameter D (mm) of the pipe segment.
κl,max and − Q t κl,min mean the upper/lower limit constraint value of natural gas pipeline κl flow at time t, (m3/h), respectively. S k,l represents the gas flow direction of the natural gas pipeline κl flow at time t, (m3/h). In this paper, the flow direction of the gas in the pipeline is fixed by fixing the pressure difference of the nodes in the optimization process. v l is the flow index of natural gas, approximately value is 2 in the low-pressure pipe network. p t κ,max and p t κ,min are the upper and lower pressure of node κ, (bar), respectively.

d: EQUATION AND INEQUALITY CONSTRAINTS FOR GAS NETWORK DYNAMIC CHARACTERISTICS
The dynamic characteristics of the gas network refer to that due to the slow transmission speed and compressibility of natural gas, the natural gas flow at the first end of the pipeline is often different from the output flow at the end. This part of natural gas flow will be temporarily stored in the natural gas pipeline, which is called ''pipe storage'' [14]. In addition, to make reasonable use of pipe storage, the stock of natural gas pipes after one week of operation is usually restored to the initial value and a certain adjustment margin is reserved for the next dispatching cycle, which can buffer the fluctuation of natural gas load and ensure the reliability of natural gas supply. The details are as follows: where M t κl is the pipe storage of natural gas pipeline κl at time t, (m 3 /h). M t κl,max andM t κl,min are the upper and lower limit constraint values of pipe storage for pipeline κl at time t, (m 3 /h), respectively. M 0 κl and M T κl represent the initial value of pipe storage and that after one cycle operation for pipeline κl at time t, (m 3 /h), respectively. κl means the collection of natural gas pipeline κl.

e: INEQUALITY CONSTRAINTS ON THE TRANSMISSION POWER BETWEEN THE GAS NETWORK AND EACH IES
GAS_IESi,max (13) where, P t GAS_IESi,max and P t GAS_IESi,min are the upper and lower limit of the transmission flow between the gas network and each IES in the park at time t, (m 3 /h), respectively.

3) EQUALITY AND INEQUALITY CONSTRAINTS OF HEAT NETWORK
In this paper, a simplified 12-node heat network system is considered [13], and the energy loss caused by frictional resistance during heat distribution pipeline transmission is also considered. The details are as follows: where, Q t m is the thermal power injected into node m of the heat network at time t, (kW). Ql t m,u=hl denotes the thermal power load of node m at time t, (kW). Qg t m,j means the heating power of node m at time t, (kW). If j ∈ HNS, the heating power is the power of the upper heat network injected into the heat distribution network. P t HEAT_IESi represents the transmission power of hot water pipelines that the heat distribution network transmits to each IES in the park at time t, (kW).

b: INEQUALITY CONSTRAINTS OF SOURCE POINT OUTPUT POWER
where, Qg t m,j,max and Qg t m,j,min are the upper and lower limit value of thermal power output of heat source point at time t, (kW), respectively.
where, Q t m,m ,min andQ t m,m ,max represent the lower and upper limit value of thermal power flow of line mm at time t, (kW), respectively. where, P t HEAT_IESi,max andP t HEAT_IESi,min are the upper and lower limit value of thermal power transmitted by hot water pipes between the heat network and each IES in the park at time t, (kW), respectively.

D. LOWER OBJECTIVE FUNCTION
The lower optimization goal is to minimize the economic operation cost of IESs in the park. The operating cost of IESs mainly includes the gas acquisition cost and the operating cost of each piece of equipment. Therefore, the lower level of the optimization objective function F IES is shown as follows: Here, C t i,IES is the economic operating cost of i-th IES, (¥). C t i,FU represents the gas purchase cost of MT in i-th IES, (¥). C t i,OM means the total operating cost of integrated energy supply and storage equipment in i-th IES at time t, (¥). is the number set of IESs in the park.
where, R FU is the price of natural gas, (¥/m 3 ). H FU is the low calorific value of natural gas, (kWh/m 3 ). η MT and η GB represent the electricity generation efficiency of MT consumption of natural gas and the heat generation efficiency of boiler consumption of natural gas, respectively. indicates the number of MT equipment.
Here, ω is the waste heat distribution coefficient. P t i,R_DL and P t i,S_DL in (21) are the rigid and shiftable electrical loads, (kW), respectively. Q t i,R_GL and Q t i,S_GL are the rigid and shiftable gas loads, (kW), respectively. Q t i,R_CL and Q t i,S_CL are the rigid and shiftable cold loads, (kW), respectively.
are the shiftable electrical/gas/cold energy participating in the virtual energy output increase/ decrease power of the integrated DR program, (kW), respectively. Q t i,A_HL is the temperature-controlled hot air power load that the park i supplies to end users at time t, (kW). Q t i,W_HL is flexible hot water power load that the park i supplies to end users at time t, (kW).

b: ENERGY COUPLING CONSTRAINTS BETWEEN ENERGY SUPPLY DEVICES
where, P t i,MT and Q t i,G2MT are the output power of MT and GB in each IES at time t, (kW), respectively. Q t i,H-REC /Q t i,H-REC,ex represents the waste heat power input/output of waste heat boiler at time t, (kW), respectively. Q t i,G2GB /Q t i,GB,ex means the thermal power input/output of GB at time t, (kW), respectively. COP MT indicates the efficiency of converting natural gas to electricity for MT. COP GB and COP REC are the thermal conversion efficiency of GB and waste heat boiler, respectively.
where, P t i,χ ,MT,max and P t i,χ ,MT,min are the upper and lower limit of the output power of MT at time t, (kW), respectively. R i,χ ,MT,dn and R i,χ ,MT,up represent the upper and lower limit of climbing force constraint coefficient where, Q t EC,ex,max /Q t EC,ex,min , Q t EH,ex,max /Q t EH,ex,min , Q t AC,ex,max /Q t AC,ex,min , and Q t H-EX,ex,max /Q t H-EX,ex,min represent the upper/ lower limits of the output power of EC, EH, AC, and H-EX at time t, (kW), respectively.

g: INEQUALITY CONSTRAINTS ON THE TRANSMISSION POWER OF LINK LINE BETWEEN EACH IES AND URBAN DISTRIBUTION NETWORK
where, P t grid_IESi,max , P t gas_IESi,max , P t heat_IESi,max , and P t cold_IESi,max are the upper limit of the transmission power of the connection line between each IES and urban distribution network at time t, (kW), respectively.

h: CARBON EMISSION CONSTRAINTS OF EACH PARK IES BASED ON REAL-TIME ENVIRONMENTAL MONITORING
where, τ cchp e is the carbon emission intensity coefficient of integrated energy supply of each IES in unit time, (t/kWh). D δ carbon,park,i is the upper limit of IES's total carbon emission quota constraint value, (t). Each historical total carbon emission quota constraint value (D δ carbon,park,1 , D δ carbon,park,2 , and D δ carbon,park,3 ) can be obtained based on real-time air environment monitoring. When the historical accumulated carbon emission data of each park is obtained, the historical method based on historical carbon emission intensity can be used to determine the initial carbon emission constraint target value [35]. E t,δ carbon,park,i is the real-time carbon emission of each IES optimized operation, (t/h). The third to fifth lines in (28) indicate that the unified calculation of real-time carbon emission E t,δ carbon,park,1 /E t,δ carbon,park,2 /E t,δ carbon,park,3 of each IES in the lower park after optimization and the upper limit of the total carbon emission quota constraint value of each IES will be transmitted to the upper urban.

i: ELECTRICITY, NATURAL GAS AND COLD LOAD SHIFTING
The electrical, gas and cold loads of the park IESs comprise the basic and shiftable loads, among which the total shiftable electrical, gas and cold load increase or decrease should equal the sum of the load change. The specific constraints of shiftable electricity, gas and cold loads are as follows:  (29) represent the allowable maximum change power of the shiftable electrical, gas and cold load, (kW). ξ t in,dl , ξ t de,dl , ξ t in,gl , ξ t de,gl , ξ t in,cl and ξ t de,cl are binary 0-1 variables. We adopt the equivalent thermal parameter model [29] based on circuit simulation for the study of building air temperature control loads in the park. Therefore, based on this equivalent model, the constraints of the power equivalent model of the flexible cooling or heating supply for end-users in the park can be obtained as follows: where Q t i,air,fore is the predicted cold and heat load demand power at time t, (kW). The second constraint in (30) means that the sum of the heating and cooling virtual energy supply is equal to the sum of the power demand. T t i,in,opt is the comfortable indoor temperature in the office buildings of the park at time t, ( • C). T t i,in and T t i,out in (30) are the indoor and outdoor air temperatures of buildings, ( • C). R is the equivalent thermal resistance of buildings, ( • C/kW). T t i,in,min and T t i,in,max are the allowable minimum and maximum indoor temperature fluctuations change, • C, respectively. The third and fourth constraints in (30) mean that the cold and heat power in the park can be flexibly regulated by the air temperature values. T t i,air is the difference between the comfortable indoor temperature and the outdoor temperature at time t ( • C). ξ t air,− and ξ t air,+ are the fluctuation parameters of cold and heat loads.
Considering the widespread application of indoor hot water supplies in buildings, we establish the constraints of the indoor flexible heating water supply load power equivalent model [29] and the specific constraints as follows:   where the first constraint in (31) is the forecasted indoor flexible hot water demand that can maintain an optimal water storage temperature for users. C water is the hot water parameter (kWh/(L · • C)), and V t i,ws,cold is the storage volume of cold water instead of hot water at time t (L). T t i,ws,cold is the temperature at time t when cold water replaces hot water at time t ( • C). T t i,ws is the storage temperature of hot water at time t ( • C), T t i,ws,max and T t i,ws,min are the allowable minimum and maximum hot water storage temperatures ( • C), respectively. T t i,ws,opt is the comfortable indoor hot water temperature value at time t ( • C). The third constraint in (31) ensures that the sum of the hot virtual power supply is equal to the sum of the  users' demand. Q t i,ws,fore represents the predicted hot water load demand power (kW). ξ t ws,− and ξ t ws,+ are the fluctuation coefficient of hot water load. T t i,ws,min and T t i,ws,max are the allowable minimum and maximum hot water storage temperature at time t ( • C).

III. BI-LEVEL DECENTRALIZED OPTIMIZATION SOLUTION STRATEGY
The bi-level optimization model constructed in this paper aims to optimize the network loss of the urban at the upper level and minimize the comprehensive economic operation costs of park IESs at the lower level.

A. LINEARIZATION OF THE MODEL
Due to the existence of nonlinear terms in formula (1) and (4) in the model, to reduce the difficulty of optimization solution, we consider that the urban RIES nonlinear optimization problem with park IESs can be transformed into a mixed integer optimization problem to solve it. Then, Second-order cone programming (SOCP) is used to perform equivalent transformation for the quadratic equation model in LinDistFlow model of the distribution network [12]. However, for the quadratic equation model in the distribution network model, the piecewise linearization method can be adopted to perform the equivalent transformation in the original model In addition, considering the complexity of the first line of formula (11) in the upper optimization constraints of urban RIES, we refer to the voltage amplitude processing method of each node in the distribution network constraints for equivalent processing, that is, we consider replacing the quadratic pressure difference of natural gas nodes in the distribution network with a quadratic pressure difference p t k − p t l .

B. APPLICATION OF ANALYTICAL TARGET CASCADING
Analytical target cascading (ATC) is an effective method to quickly solve decentralized [37], hierarchical coordination problems. It allows each subject in the hierarchy to make independent decisions, coordinates and optimizes the decision between each subject and each sub-subject, and obtains the overall optimal solution of the system. ATC transfers the optimized coupling variables of each principal system to the objective optimization function of each sub-principal system. Then the optimized coupling variables of each sub-principal system are transferred to the objective optimization function of each principal system. Therefore, the coupling variables of each link line are decoupled. The first and second terms of the Lagrange penalty function are introduced into the objective function of each agent, respectively, to ensure the optimization of each agent. Compared with other optimization methods, the objective cascade method has the advantages of parallel optimization, unlimited series, and strict proof of convergence. Therefore, we adopt the improved ATC algorithm to solve the above bi-level decentralized coordinated optimal scheduling model, and the optimization problem is described as follows:    Comparison of the total output power of devices with carbon emission sources in the upper urban regional from scenarios one to three.  Here, P t λ,GD_EHi , P t γ ,GAS_EHi , and P t φHEAT_EHi represent the coupling variable of transmission power between the upper power, gas, and heat distribution network and each IES in the lower park at time t, (kW), respectively. ω t λ , ω t γ , and ω t φ mean the first term multiplier of the Lagrange penalty function for the objective function of power, gas, and heat distribution network at time t; ς t λ , ς t γ , and ς t φ are the quadratic multiplier of Lagrange penalty function for the objective function of power, gas, and heat distribution network at time t, respectively.

E. THE CONVERGENCE CRITERION AND LAGRANGE MULTIPLIER RENEWAL PRINCIPLE
The optimized coupling variables of each principal system in upper/lower regions can be iteratively transferred to each other by ATC to ensure the optimal goal of each principal system. Therefore, whether the difference of the optimized coupling variable and objective function between the upper/lower-level system meets the accuracy requirement as the convergence condition is considered in this paper, the details are as follows: Here, ε 1 , ε 2 , and ε 3 represent the convergence accuracy of the difference of coupling variable and the difference of objective function between the upper and lower systems in the urban area, respectively.
If the above convergence criterion cannot satisfy the formula at the same time, the Lagrange penalty function multiplier should be updated according to the following formula (37), the details are as follows: To accelerate the convergence speed of ATC, the value of β in (37) is set as 2.5, and the value of the Lagrange penalty function multiplier in (37) is uniformly set as 1.5.

F. THE SOLUTION FLOW OF BI-LEVEL DECENTRALIZED OPTIMIZATION ALGORITHM
The process of urban RIES bi-level decentralized coordination optimization algorithm based on ATC with IESs in the park is shown as follows: Step 1: Input carbon emission constraint value distribution parameters and system equipment parameters of urban RIES. Set the initial value of the system coupling variable, the initial value of the first and second multiplier of the Lagrange penalty function, and the initial value of the iteration k = 1. The primary variable zero matrix is reserved for data recording observation.
Step 2: According to formula (35), (18), and (31), the optimization problems of IES in the lower park are solved respectively. Then, the solved coupling variables (P t λ,grid_IESi /P t γ ,gas_IESi /P t φ,heat_IESi ) are transmitted to the upper power/gas/heat distribution network system for optimization, and the real-time carbon emissions (E t,δ carbon,park,1 /E t,δ carbon,park,2 /E t,δ carbon,park,3 ) of IESs in the lower park after optimization are transmitted to the upper. And then update the carbon emission constraint value (D δ carbon,park,1 /D δ carbon,park,2 /D δ carbon,park,3 ) of IESs decomposed from the upper urban area to the lower park.
Step 3: When the upper power distribution system receives all the data transmitted by each IES in the lower park. According to the formula (32), (4) ∼  (14) ∼ (17). The optimization problem of the upper urban power/gas/heat distribution network system is solved, and the real-time carbon emission (E t,δ carbon,e 1 /E t,δ carbon,e 2 ) after the optimization of the power distribution network is transmitted to the upper layer, to update the carbon emission constraint value (D δ carbon,e 1 /D δ carbon,e 2 ) decomposed by the urban itself.
Step 4: Refer to Formula (36) to check whether the convergence conditions of the optimization algorithm meet the convergence criteria. If convergence, the iteration is terminated and output the optimal scheduling results. If not, set k = k + 1, update the multiplier of Langerin penalty function by Formula (37), then return to Step 2, and then continue iterative solution.
The flow chart of the proposed bi-level iterative algorithm is shown in figure 4.

IV. CASE STUDIES A. CASE DESCRIPTION
This section establishes a bi-level decentralized optimal lowcarbon economic dispatch model for urban RIES with park IESs considering the optimal allocation mechanism of carbon emissions. The model and algorithm in this paper are written in Matlab2017a and run on a computer with an Intel Core i5 5257U CPU, 3.00 GHz main frequency and 8 GB memory. The correctness and validity of the proposed method are verified by three cases of different scenarios. In this case, the urban area as a whole is divided into a bi-level decentralized coordination structure of upper RIES and lower IESs. Among them, the upper urban considers three network structures: power, gas, and heat distribution network. Each distribution network has its conventional loads, the node 24 and 29 of the network are connected to two dispatchable coal-fired generating units (gm 1 and gm 3 ), node 7 is connected to a diesel generating unit (gm 2 ), and node 32 is connected to a set of battery storage (DN-battery) equipment.
Three different electrical/gas/heat/cold rigid loads and flexible loads are also considered in the lower park area. It is assumed that the total power of electricity/gas/ heat/cooling predicted by the lower level IES is equal, which is 4095, 3024, 5364 and 4335KW respectively. The flexible load of electricity/gas/cold in the lower layer of the park is shiftable load, and the flexible heat load includes flexible air heating and hot water supply temperature loads, which account for 60% and 30%, respectively. In addition, the wind and photovoltaic power generation equipment (WT and PV) are also contained   in the park IESs, which are typically represented by the predicted curves of load and rated output power of WT and PV on a certain day, as shown in figure 5. Refer to tables 1 and 2 for the prediction curves of outdoor temperature and storage volume of cold water instead of hot water in different urban areas on a certain day. The equipment parameters of the calculation example are shown in table 3, and the initial carbon emission constraint value of urban RIES (t/day), as shown in table 4 [35]. The symbols and descriptions of the main equipment types of IES in the lower park area are shown in section 2.1, and the Lagrange function multiplier parameter of ATC algorithm is shown in formula (37).
In fact, the carbon emissions generated by the whole urban depend on the operation of equipment with carbon emission sources in the upper and lower areas, respectively, which depends on the load distribution of the upper urban distribution network itself and the comprehensive energy consumption of IES in the lower park.

B. SCENARIOS DESCRIPTION
In this paper, the following 3 different scenarios are set for the analysis of the network loss of upper urban regional power, gas, and heat distribution network, the economic operating cost and corresponding comprehensive energy consumption of park IESs, the carbon emissions of each system in the urban RIES, and the convergence characteristics of the algorithm. By analyzing the differences between the data of multiple types of indicators in the above different scenarios, it is helpful for urban to build a reasonable regional low-carbon RIES. The details are shown in table 5.
In table 5, scenario 1 does not consider the overall carbon emission constraints of urban RIES, while scenario 2 considers the overall carbon emission constraints of the urban based on real-time environmental monitoring of urban RIES. Based on scenario 2, scenario 3 further considers the end users of park IESs to participate in integrated demand response (IDR).  Table 6 shows the comparison and analysis results of network losses of the upper urban power distribution/gas/heat network system in scenarios 1, 2, and 3. Table 7 shows the comparison and analysis results of system operating cost of IES1∼3 in scenarios 1, 2, and 3. Table 8 shows the comparison and analysis results of carbon emissions of the upper urban distribution system and the lower park IESs in scenarios 1,2 and 3. Results a and b in tables 6, 7, and 8 represent the comparison results of scenarios 1 and 2, and 1 and 3, respectively.
In table 6, compared with scenarios 2 and 1 (see the result a), the network loss of the upper urban regional distribution network system is reduced by about 2.3%, and the network loss of the gas distribution network system is almost unchanged. However, the network loss of the heat distribution network system is increased by about 14.8%, which also indirectly led to the total network loss of the upper urban regional distribution network system increased by about 0.5%. Compared with scenario 1 (see result b), the network loss of power and gas distribution network system is reduced by 15.9% and 1.9%, respectively. Although the network loss of heat distribution network system is still increased by 9.1%, however, compared with scenario 2, scenario 3 reduces the network loss of the upper urban power, gas, and heat distribution network system by about 13.9%, 1.9%, and 5.1% through introducing IDR to end users in the park, respectively. Through data analysis, even based on scenario 1, the network loss of the heat distribution network system in scenario 3 still increases. However, the total network loss of the urban RIES in scenario 3 is effectively reduced. VOLUME 10, 2022 Therefore, the strategy in scenario 3 is better than that in scenario 2.
In table 7, compared with scenario 1, the economic operating costs of park IES1 in scenario 2 are reduced by about 4.6%, those of park IES2 are increased by about 2.6%, while those of IES3 are virtually unchanged. Although the total operating costs of park IESs are reduced by about 1.2%, the effect is not obvious (see result a). Moreover, the operating costs of each IES in the lower park area are greatly reduced (see result b), and the total operating costs of the lower park IESs are reduced by about 89.6%. Therefore, the operating costs of park IESs in scenario 3 are optimal.
In table 8, scenario 1 shows that the upper urban is responsible for the initial carbon emission constraint upper limit values of park IESs and the urban itself. When only the carbon emission constraints are considered, the carbon emissions generated by equipment with carbon emission sources in the upper urban and lower park in scenario 2 are reduced by about 1.2 t/day and 0.03 t/day, respectively. However, according to table 8, it is found that the carbon emission reduction effect in the middle and lower levels of the park in scenario 2 is not ideal (see result a). Therefore, the carbon emissions of the upper urban RIES and the lower park IESs are further optimized in scenario 3. And, the carbon emissions generated by IESs in the lower park area are almost zero (see result b). Therefore, scenario 3 has the best carbon emission reduction effect.
Overall, the bi-level carbon-constrained optimization allocation mechanism is an effective IDR strategy, which can greatly reduce total network loss of the upper urban RIES, decrease the operating cost of park IESs, and further reduce the overall carbon emissions of urban. Even, to a certain extent, zero carbon emissions can be achieved in the park, which can ultimately ensure the overall environment of the urban.
2) THE TOTAL OUTPUT POWER OF GM1, GM2, AND GM3 NUMERICAL ANALYSIS Table 9 shows the total output power of main devices (gm 1 , gm 2 , and gm 3 ) with carbon emission sources in the upper urban RIES, which are described as follows.
In table 9, in scenarios 2 and 3, the output power of the dispatchable diesel generator unit (gm 2 ) decreases by about 70.4% and 76.1%, respectively. Then, the output power of the coal-fired generator unit (gm 3 ) decreases by about 2% and 11.2%. Compared with scenario 2, the total output power of generator sets (gm 1 , gm 2 , and gm 3 ) with carbon emission sources in the power distribution system of upper urban in scenario 3 is further reduced by about 6.8%. Table 10 shows the total output power of devices (MT, H-REC, and GB) with carbon emission sources in the lower park area of IES1, IES2, and IES3.

3) THE TOTAL OUTPUT POWER OF MT, H-REC, AND GB NUMERICAL ANALYSIS
Combining the analysis in tables 10, 8, and 9, in scenario 2, the total output power of each device with carbon emission source in each park may decrease, increase or remain unchanged. However, in scenario 3, the total output power of equipment with carbon emission sources in the park will reduced by about 100%, which can theoretically achieve zero carbon emission in the park.
Combining the analysis in table 10 and figure 8, the IDR strategy can smooth the energy consumption curve of each end user relatively, and achieve the effect of peak clipping and valley filling, and greatly improves the user's energy consumption comfort (refer to figure 8). Figure 6 shows the comparison of real-time charging/ discharging power of energy storage devices on the distribution network side, the real-time active power injected by power nodes in the power distribution network, the real-time natural gas flow output by natural gas source points in the gas distribution network, and the real-time thermal power output by heat source points in the heat distribution network under different scenarios. Figure 8 shows the optimization dispatching results of real-time power balance output power curves of each IES in the lower park under different scenarios. Figure 9 shows the real-time charging/discharging state of battery, cold/ heat/gas storage devices of each IES in the lower park under different scenes.

4) THE SYSTEM OPERATING RESULTS COMPREHENSIVE ANALYSIS
Combining the analysis in table 9 and figure 7, the realtime status of scenarios 2 and 1 are basically the same. In scenario 3, the charging/discharging status of the battery storage device on the side of the upper urban RIES. However, the real-time active power injected by the power supply and heat source nodes has a certain influence on the output of the devices with carbon emission sources, while the real-time gas flow injected by the natural gas source node is almost unchanged.
According to figure 9, the battery device of IES3, cold and gas storage tank devices of IES1∼3, and heat storage tank device of IES1∼2 in the lower level of the park are basically not in operation. However, the battery device of IES1∼2 has a more frequent charge/discharge state, and only the heat storage tank device of park IES3 stored and released heat for part of the time.
Therefore, we can basically confirm that the regulation of energy storage devices in the upper urban RIES and the lower park IESs is an indirect factor to reduce the overall urban regional carbon emissions. However, the regulation of electricity and heat power injected by power nodes and heat nodes of the distribution network is the main factor to reduce the urban regional overall carbon emissions. In addition, we also need to specifically analyze how the upper urban power, gas, and heat distribution network system optimally discretely coordinates the operation of each IES in the park, as shown in figure 6.
Combining the analysis of figure 6 and figure 8, the operating state of scenario 3 is the best. Compared with scenario 2, the energy load of each IES in the lower park area tends to be stable. Among them, the gas/cold power transmitted by the upper urban RIES to the lower park IESs can mostly meet the use of gas/cold energy users. The heat power used by thermal users is satisfied by the heat transmission power of the heat distribution network system and EH equipment, while, the power of electricity users is met by the transmission power of the distribution network system, WT, and PV devices, and so on. Through comparison, it is found that all types of energy storage devices in the lower park are barely used, and most energy conversion devices (such as EC, AC, HX, etc.) are not used. On the one hand, it is because we have introduced the bi-level optimal control strategy of carbon emission in urban areas. On the other hand, it is because we have introduced IDR program for electricity, gas, cold, and heat for all end users in the park area. Therefore, we can better realize the minimum network loss of the upper urban RIES, the minimum operation cost of the lower park IESs, and maintain the overall environment of the urban region to maintain the best.

D. CONVERGENCE OF BI-LEVEL DECENTRALIZED OPTIMAL METHOD
In section three, the improved ATC algorithm is proposed to solve the bi-level iterative decentralized optimization problem, and the parameters of the above method are given. In this paper, figure 10 shows the convergence property of coupling variables between the upper urban RIES and the lower park IESs under different scenarios. Figure 8 shows the convergence characteristic curves of objective function difference iteration of each park IES under different scenarios.
In figures 10 and 11, the improved ATC algorithm proposed by us can converge in different scenarios. Among them, scenarios 1, 2, and 3 completed iterative convergence at the 118th, 10th and 14th times, respectively. Therefore, scenario1 has the slowest convergence speed, followed by scenarios 3 and 2. This is because we have introduced the bi-level optimization allocation model of carbon emission constraints in scenarios 2 and 3, which can reduce the constraint space of the original system optimization operation. Scenario 3 introduces the IDR strategy on the basis of scenario 2, which will inevitably lead to the increase of a large number of interaction information of coupling variables in the iterative process, especially affecting the convergence rate of optimization of complex urban RIES. Therefore, the convergence speed of scenario 3 is lower than that of scenario 2. Even if the urban RIES has a large number of initial parameters and iterative coupling variables, scenarios 2 and 3 can quickly realize iterative convergence without sacrificing accuracy. However, through the comprehensive comparison of the network loss of urban power, gas, and heat distribution network, the operating cost, and integrated energy conversion efficiency of park IESs, and the overall carbon emissions of the urban area, even if the convergence speed of scenario 3 is slightly lower than that of scenario 2, the effectiveness of the proposed urban RIES bi-level decentralized optimization scheduling model and the improved ATC algorithm considering environmental constraints can be fully demonstrated.

V. SUMMARY AND CONCLUSION
This paper proposes a bi-level decentralized optimal lowcarbon economic dispatch model for urban RIES with consideration of park IESs under the background of China's urban energy reform. The improved ATC method is adopted to solve the bi-level decentralized iterative optimization model, and three different types of scenarios in cases are used to verify the effectiveness of this method.
The conclusions of this paper are summarized as follows: 1) Real-time air environment monitoring is adopted in the urban RIES. According to historical data, the carbon emissions of the unit integrated energy supply with carbon sources of urban RIES and park IESs are constrained in real-time. In addition, a bi-level carbon emission constrained optimal allocation mechanism is constructed, which can determine the actual carbon emissions and help the urban meet the overall environmental requirements.
2) By decentralizing and coordinating the optimized operation of urban RIES and each park IES, the network loss of the upper urban regional power and gas distribution network system is reduced by about 15.9% and 1.9% respectively. Although the network loss of the heat distribution network is increased by about 9.1%, the total urban regional network loss and carbon emissions are further reduced. In addition, the total operating cost of the lower park system is reduced by about 1.2%, which can almost achieve zero carbon emission of park IESs.
3) The IDR program includes electricity, cold, gas, and heat load shifting, a flexible cooling or heating power supply, and a flexible indoor hot water supply, which can not only further reduce the operating cost and carbon emissions of each park IES but also improve the comfort level of end users and energy conversion efficiency of the park IESs. 4) An improved ATC algorithm is proposed, which can not only arbitrarily set the initial value parameters of coupling variables between the upper and lower systems but also guarantee the privacy of urban RIES. Moreover, the feasible region constraints can be continuously optimized and compressed to achieve convergence quickly without sacrificing convergence accuracy.

VI. FUTURE WORK
Currently, there is no unified definition of integrated energy efficiency of urban RIES. Therefore, how to build an effectively integrated energy efficiency theory for the urban or park RIES will be the focus of this study in the future.